--- /dev/null
+//----------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// 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
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------
+#ifndef __deal2__fe_dgp_monomial_h
+#define __deal2__fe_dgp_monomial_h
+
+#include <base/config.h>
+#include <base/polynomials_p.h>
+#include <grid/geometry_info.h>
+#include <fe/fe.h>
+
+template <int dim> class MappingQ;
+
+
+/*!@addtogroup fe */
+/*@{*/
+
+/**
+ * Discontinuous finite elements based on monomials of degree up to
+ * <tt>k</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.
+ *
+ * @author Ralf Hartmann, 2004
+ */
+template <int dim>
+class FE_DGPMonomial : public FiniteElement<dim>
+{
+ public:
+ /**
+ * Constructor for the polynomial
+ * space of degree @p{k}.
+ */
+ FE_DGPMonomial (const unsigned int k);
+
+ /**
+ * 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.
+ */
+ 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,
+ * 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 @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,
+ * 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.
+ *
+ * This function returns
+ * @p{true}, if the shape
+ * function @p{shape_index} has
+ * non-zero values on the face
+ * @p{face_index}.
+ *
+ * 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;
+
+ 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:
+
+ /**
+ * 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(unsigned int degree);
+
+ /**
+ * Initialize the embedding
+ * matrices. Called from the
+ * constructor.
+ */
+ void initialize_embedding ();
+
+ /**
+ * Initialize the restriction
+ * matrices. Called from the
+ * 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
+ * access to build_renumbering
+ * function.
+ */
+ friend class MappingQ<dim>;
+};
+
+/*@}*/
+
+template <int dim>
+inline unsigned int
+FE_DGPMonomial<dim>::get_degree () const
+{
+ return polynomial_space.degree();
+}
+
+
+#endif
--- /dev/null
+//----------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// 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
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------
+
+#include <base/quadrature.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_dgp_monomial.h>
+#include <fe/fe_values.h>
+
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#endif
+
+
+// TEST
+#include <base/polynomial.h>
+
+
+
+// namespace for some functions that are used in this file.
+namespace
+{
+ // storage of hand-chosen support
+ // points
+ //
+ // For dim=2, dofs_per_cell of
+ // FE_DGPMonomial(k) is given by
+ // 0.5(k+1)(k+2), i.e.
+ //
+ // k 0 1 2 3 4 5 6 7
+ // dofs 1 3 6 10 15 21 28 36
+ //
+ // indirect access of unit points:
+ // the points for degree k are
+ // located at
+ //
+ // points[start_index[k]..start_index[k+1]-1]
+ const unsigned int start_index2d[6]={0,1,4,10,20,35};
+ const double points2d[35][2]=
+ {{0,0},
+ {0,0},{1,0},{0,1},
+ {0,0},{1,0},{0,1},{1,1},{0.5,0},{0,0.5},
+ {0,0},{1,0},{0,1},{1,1},{1./3.,0},{2./3.,0},{0,1./3.},{0,2./3.},{0.5,1},{1,0.5},
+ {0,0},{1,0},{0,1},{1,1},{0.25,0},{0.5,0},{0.75,0},{0,0.25},{0,0.5},{0,0.75},{1./3.,1},{2./3.,1},{1,1./3.},{1,2./3.},{0.5,0.5}
+ };
+
+ //
+ // For dim=3, dofs_per_cell of
+ // FE_DGPMonomial(k) is given by
+ // 1./6.(k+1)(k+2)(k+3), i.e.
+ //
+ // k 0 1 2 3 4 5 6 7
+ // dofs 1 4 10 20 35 56 84 120
+ const unsigned int start_index3d[6]={0,1,5,15/*,35*/};
+ const double points3d[35][3]=
+ {{0,0,0},
+ {0,0,0},{1,0,0},{0,1,0},{0,0,1},
+ {0,0,0},{1,0,0},{0,1,0},{0,0,1},{0.5,0,0},{0,0.5,0},{0,0,0.5},{1,1,0},{1,0,1},{0,1,1}
+ };
+
+
+ template<int dim>
+ void generate_unit_points (const unsigned int,
+ vector<Point<dim> > &);
+
+ template <>
+ void generate_unit_points (const unsigned int,
+ vector<Point<1> > &)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+
+ template <>
+ void generate_unit_points (const unsigned int k,
+ vector<Point<2> > &p)
+ {
+ Assert(k<=4, ExcNotImplemented());
+ Assert(p.size()==start_index2d[k+1]-start_index2d[k], ExcInternalError());
+ for (unsigned int i=0; i<p.size(); ++i)
+ {
+ p[i](0)=points2d[start_index2d[k]+i][0];
+ p[i](1)=points2d[start_index2d[k]+i][1];
+ }
+ }
+
+ template <>
+ void generate_unit_points (const unsigned int k,
+ vector<Point<3> > &p)
+ {
+ Assert(k<=2, ExcNotImplemented());
+ Assert(p.size()==start_index3d[k+1]-start_index3d[k], ExcInternalError());
+ for (unsigned int i=0; i<p.size(); ++i)
+ {
+ p[i](0)=points3d[start_index3d[k]+i][0];
+ p[i](1)=points3d[start_index3d[k]+i][1];
+ p[i](2)=points3d[start_index3d[k]+i][2];
+ }
+ }
+}
+
+
+
+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)
+{
+ Assert(polynomial_space.n()==dofs_per_cell, ExcInternalError());
+
+ // DG doesn't have constraints, so
+ // leave them empty
+
+ // initialize the interpolation
+ // matrices
+ initialize_embedding ();
+// initialize_restriction ();
+
+ // note, that these elements have
+ // neither support nor face-support
+ // points, so leave these fields
+ // empty
+}
+
+
+
+template <int dim>
+std::string
+FE_DGPMonomial<dim>::get_name () const
+{
+ // note that the
+ // FETools::get_fe_from_name
+ // function depends on the
+ // particular format of the string
+ // this function returns, so they
+ // have to be kept in synch
+
+#ifdef HAVE_STD_STRINGSTREAM
+ std::ostringstream namebuf;
+#else
+ std::ostrstream namebuf;
+#endif
+
+ namebuf << "FE_DGPMonomial<" << dim << ">(" << get_degree() << ")";
+
+#ifndef HAVE_STD_STRINGSTREAM
+ namebuf << std::ends;
+#endif
+ return namebuf.str();
+}
+
+
+
+template <int dim>
+FiniteElement<dim> *
+FE_DGPMonomial<dim>::clone() const
+{
+ return new FE_DGPMonomial<dim>(get_degree());
+}
+
+
+
+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>::
+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
+ // DGPMonomial element
+ AssertThrow ((x_source_fe.get_name().find ("FE_DGPMonomial<") == 0)
+ ||
+ (dynamic_cast<const FE_DGPMonomial<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_DGPMonomial<dim> &source_fe
+ = dynamic_cast<const FE_DGPMonomial<dim>&>(x_source_fe);
+
+ const unsigned int m=interpolation_matrix.m();
+ const unsigned int n=interpolation_matrix.n();
+ Assert (m == this->dofs_per_cell, ExcDimensionMismatch (m, this->dofs_per_cell));
+ Assert (n == source_fe.dofs_per_cell, ExcDimensionMismatch (n, source_fe.dofs_per_cell));
+
+ const unsigned int min_mn=
+ interpolation_matrix.m()<interpolation_matrix.n() ?
+ interpolation_matrix.m() : interpolation_matrix.n();
+
+ for (unsigned int i=0; i<min_mn; ++i)
+ interpolation_matrix(i,i)=1.;
+}
+
+
+
+template <int dim>
+void
+FE_DGPMonomial<dim>::initialize_embedding ()
+{
+ vector<Point<dim> > unit_points(this->dofs_per_cell);
+ generate_unit_points(get_degree(), unit_points);
+
+ 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)
+ {
+ const Point<dim> &p_subcell=unit_points[j];
+
+ 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)
+ {
+ cell_interpolation(j,i) = polynomial_space.compute_value (i, p_cell);
+ subcell_interpolation(j,i) = polynomial_space.compute_value (i, p_subcell);
+ }
+ }
+
+ // then compute the embedding
+ // matrix for this child and
+ // this coordinate direction
+ subcell_interpolation.gauss_jordan ();
+ subcell_interpolation.mmult (this->prolongation[child], cell_interpolation);
+
+ // cut off very small values
+ 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*get_degree()*dim)
+ this->prolongation[child](i,j) = 0.;
+ }
+}
+
+
+template <int dim>
+void
+FE_DGPMonomial<dim>::initialize_restriction ()
+{
+ Assert(false, ExcNotImplemented());
+}
+
+
+//----------------------------------------------------------------------
+// Auxiliary functions
+//----------------------------------------------------------------------
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_DGPMonomial<dim>::get_dpo_vector(unsigned int deg)
+{
+ std::vector<unsigned int> dpo(dim+1, 0U);
+ dpo[dim] = ++deg;
+ for (unsigned int i=1;i<dim;++i)
+ {
+ dpo[dim] *= deg+i;
+ dpo[dim] /= i+1;
+ }
+ return dpo;
+}
+
+
+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 <>
+bool
+FE_DGPMonomial<1>::has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ return face_index==1 || (face_index==0 && get_degree()==0);
+}
+
+#endif
+
+#if deal_II_dimension == 2
+
+template <>
+bool
+FE_DGPMonomial<2>::has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ bool support_on_face=false;
+ if (face_index==1 || face_index==2)
+ support_on_face=true;
+ else
+ {
+ unsigned int degrees[2];
+ polynomial_space.directional_degrees(shape_index, degrees);
+ if ((face_index==0 && degrees[1]==0) ||
+ (face_index==3 && degrees[0]==0))
+ support_on_face=true;
+ }
+ return support_on_face;
+}
+
+#endif
+
+#if deal_II_dimension == 3
+
+template <>
+bool
+FE_DGPMonomial<3>::has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ bool support_on_face=false;
+ if (face_index==1 || face_index==3 || face_index==4)
+ support_on_face=true;
+ else
+ {
+ unsigned int degrees[3];
+ polynomial_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))
+ support_on_face=true;
+ }
+ return support_on_face;
+}
+
+#endif
+
+
+template <int dim>
+unsigned int
+FE_DGPMonomial<dim>::memory_consumption () const
+{
+ Assert (false, ExcNotImplemented ());
+ return 0;
+}
+
+
+
+template class FE_DGPMonomial<deal_II_dimension>;