--- /dev/null
+//---------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004 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_poly_h
+#define __deal2__fe_poly_h
+
+
+#include <fe/fe.h>
+
+
+/**
+ * This class gives a unified framework for implementing a
+ * FiniteElement class based on a polynomial space like a
+ * @p TensorProductSpace or a @p PolynomialSpace.
+ *
+ * Every class that implements following functions can be used as
+ * template parameter @p POLY. Example classes are @p
+ * TensorProductSpace and @p PolynomialSpace.
+ *
+ * @code
+ * double compute_value (const unsigned int i,
+ * const Point<dim> &p) const;
+ *
+ * Tensor<1,dim> compute_grad (const unsigned int i,
+ * const Point<dim> &p) const;
+ *
+ * Tensor<2,dim> compute_grad_grad (const unsigned int i,
+ * const Point<dim> &p) const;
+ * @endcode
+ *
+ * This class is not a fully implemented FiniteElement class. Instead
+ * there are several pure virtual functions declared in the
+ * FiniteElement and FiniteElementBase classes which cannot
+ * implemented by this class but are left for implementation in
+ * derived classes.
+ *
+ * Todos:
+ * - checke dim of POLY
+ * - templatisiere nur auf POLY, dim ergibt sich durch POLY::dim
+ *
+ * @author Ralf Hartmann 2004
+ **/
+template <class POLY, int dim/*=POLY::dim*/>
+class FE_Poly : public FiniteElement<dim>
+{
+ public:
+ /**
+ * Constructor.
+ */
+ FE_Poly (const POLY& poly_space,
+ const FiniteElementData<dim> &fe_data,
+ const std::vector<bool> &restriction_is_additive_flags,
+ const std::vector<std::vector<bool> > &nonzero_components);
+
+ /**
+ * 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;
+
+ /**
+ * 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;
+
+
+ protected:
+
+ /**
+ * 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 ;
+
+
+ /**
+ * Determine the values that need
+ * to be computed on the unit
+ * cell to be able to compute all
+ * values required by @p{flags}.
+ *
+ * For the purpuse of this
+ * function, refer to the
+ * documentation in
+ * @p{FiniteElement}.
+ *
+ * The effect in this element is
+ * as follows: if
+ * @p{update_values} is set in
+ * @p{flags}, copy it to the
+ * result. All other flags of the
+ * result are cleared, since
+ * everything else must be
+ * computed for each cell.
+ */
+ virtual UpdateFlags update_once (const UpdateFlags flags) const;
+
+ /**
+ * Determine the values that need
+ * to be computed on every
+ * cell to be able to compute all
+ * values required by @p{flags}.
+ *
+ * For the purpuse of this
+ * function, refer to the
+ * documentation in
+ * @p{FiniteElement}.
+ *
+ * The effect in this element is
+ * as follows:
+ * @begin{itemize}
+ * @item if @p{update_gradients}
+ * is set, the result will
+ * contain @p{update_gradients}
+ * and
+ * @p{update_covariant_transformation}.
+ * The latter is required to
+ * transform the gradient on the
+ * unit cell to the real
+ * cell. Remark, that the action
+ * required by
+ * @p{update_covariant_transformation}
+ * is actually performed by the
+ * @p{Mapping} object used in
+ * conjunction with this finite
+ * element.
+ * @item if
+ * @p{update_second_derivatives}
+ * is set, the result will
+ * contain
+ * @p{update_second_derivatives}
+ * and
+ * @p{update_covariant_transformation}.
+ * The rationale is the same as
+ * above and no higher
+ * derivatives of the
+ * transformation are required,
+ * since we use difference
+ * quotients for the actual
+ * computation.
+ * @end{itemize}
+ */
+ virtual UpdateFlags update_each (const UpdateFlags flags) const;
+
+
+ /**
+ * 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;
+ };
+
+ /**
+ * The polynomial space.
+ */
+ POLY poly_space;
+};
+
+
+
+#endif
--- /dev/null
+//----------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004 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/polynomial.h>
+// #include <base/tensor_product_polynomials.h>
+// #include <grid/tria.h>
+// #include <grid/tria_iterator.h>
+// #include <dofs/dof_accessor.h>
+#include <base/quadrature.h>
+#include <base/tensor_product_polynomials.h>
+#include <base/polynomials_p.h>
+#include <fe/fe_poly.h>
+//#include <fe/fe.h>
+// #include <fe/mapping.h>
+// #include <fe/fe_q.h>
+#include <fe/fe_values.h>
+
+// #ifdef HAVE_STD_STRINGSTREAM
+// # include <sstream>
+// #else
+// # include <strstream>
+// #endif
+
+
+template <class POLY, int dim>
+FE_Poly<POLY,dim>::FE_Poly (const POLY& poly_space,
+ const FiniteElementData<dim> &fe_data,
+ const std::vector<bool> &restriction_is_additive_flags,
+ const std::vector<std::vector<bool> > &nonzero_components):
+ FiniteElement<dim> (fe_data,
+ restriction_is_additive_flags,
+ nonzero_components),
+ poly_space(poly_space)
+{}
+
+
+
+template <class POLY, int dim>
+double
+FE_Poly<POLY,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 poly_space.compute_value(i, p);
+}
+
+
+template <class POLY, int dim>
+double
+FE_Poly<POLY,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 poly_space.compute_value(i, p);
+}
+
+
+
+template <class POLY, int dim>
+Tensor<1,dim>
+FE_Poly<POLY,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 poly_space.compute_grad(i, p);
+}
+
+
+
+template <class POLY, int dim>
+Tensor<1,dim>
+FE_Poly<POLY,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 poly_space.compute_grad(i, p);
+}
+
+
+
+template <class POLY, int dim>
+Tensor<2,dim>
+FE_Poly<POLY,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 poly_space.compute_grad_grad(i, p);
+}
+
+
+
+template <class POLY, int dim>
+Tensor<2,dim>
+FE_Poly<POLY,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 poly_space.compute_grad_grad(i, p);
+}
+
+
+
+//----------------------------------------------------------------------
+// Auxiliary functions
+//----------------------------------------------------------------------
+
+
+
+
+template <class POLY, int dim>
+UpdateFlags
+FE_Poly<POLY,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 <class POLY, int dim>
+UpdateFlags
+FE_Poly<POLY,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 <class POLY, int dim>
+typename Mapping<dim>::InternalDataBase *
+FE_Poly<POLY,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;
+
+ // 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)
+ {
+ poly_space.compute(quadrature.point(i),
+ values, grads, grad_grads);
+
+ if (flags & update_values)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+ data->shape_values[k][i] = values[k];
+
+ if (flags & update_gradients)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+ data->shape_gradients[k][i] = grads[k];
+ }
+ return data;
+}
+
+
+
+
+//----------------------------------------------------------------------
+// Fill data of FEValues
+//----------------------------------------------------------------------
+
+template <class POLY, int dim>
+void
+FE_Poly<POLY,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);
+ };
+ }
+
+ const typename QProjector<dim>::DataSetDescriptor dsd;
+ if (flags & update_second_derivatives)
+ this->compute_2nd (mapping, cell, dsd.cell(),
+ mapping_data, fe_data, data);
+}
+
+
+
+template <class POLY, int dim>
+void
+FE_Poly<POLY,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 dsd;
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = dsd.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)
+ {
+ 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+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 <class POLY, int dim>
+void
+FE_Poly<POLY,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 dsd;
+ const typename QProjector<dim>::DataSetDescriptor offset
+ = dsd.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)
+ {
+ 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+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 <class POLY, int dim>
+unsigned int
+FE_Poly<POLY,dim>::n_base_elements () const
+{
+ return 1;
+}
+
+
+
+template <class POLY, int dim>
+const FiniteElement<dim> &
+FE_Poly<POLY,dim>::base_element (const unsigned int index) const
+{
+ Assert (index==0, ExcIndexRange(index, 0, 1));
+ return *this;
+}
+
+
+
+template <class POLY, int dim>
+unsigned int
+FE_Poly<POLY,dim>::element_multiplicity (const unsigned int index) const
+{
+ Assert (index==0, ExcIndexRange(index, 0, 1));
+ return 1;
+}
+
+
+
+template class FE_Poly<TensorProductPolynomials<deal_II_dimension>, deal_II_dimension>;
+template class FE_Poly<PolynomialSpace<deal_II_dimension>, deal_II_dimension>;
+template class FE_Poly<PolynomialsP<deal_II_dimension>, deal_II_dimension>;