From: Ralf Hartmann Date: Tue, 9 Mar 2004 16:07:28 +0000 (+0000) Subject: Implementation of a FE class based on the PolynomialsP polynomial space. X-Git-Tag: v8.0.0~15646 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3a2ccdc4f2b540db126918aaf1af1069897eedc4;p=dealii.git Implementation of a FE class based on the PolynomialsP polynomial space. git-svn-id: https://svn.dealii.org/trunk@8701 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_dgp_monomial.h b/deal.II/deal.II/include/fe/fe_dgp_monomial.h new file mode 100644 index 0000000000..8e19bf24dd --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_dgp_monomial.h @@ -0,0 +1,456 @@ +//---------------------------------------------------------------- +// $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 +#include +#include +#include + +template class MappingQ; + + +/*!@addtogroup fe */ +/*@{*/ + +/** + * Discontinuous finite elements based on monomials of degree up to + * k + * + * 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 +class FE_DGPMonomial : public FiniteElement +{ + 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(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 &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 &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 &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 &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 &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 &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::ExcInterpolationNotImplemented} + * is thrown. + */ + virtual void + get_interpolation_matrix (const FiniteElementBase &source, + FullMatrix &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 & + 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 *clone() const; + + /** + * Prepare internal data + * structures and fill in values + * independent of the cell. + */ + virtual + typename Mapping::InternalDataBase * + get_data (const UpdateFlags, + const Mapping& mapping, + const Quadrature& quadrature) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_face_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_subface_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& 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 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 P_k. + */ + const PolynomialsP 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::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; +}; + +/*@}*/ + +template +inline unsigned int +FE_DGPMonomial::get_degree () const +{ + return polynomial_space.degree(); +} + + +#endif diff --git a/deal.II/deal.II/source/fe/fe_dgp_monomial.cc b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc new file mode 100644 index 0000000000..4bac467c9b --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc @@ -0,0 +1,713 @@ +//---------------------------------------------------------------- +// $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 +#include +#include +#include +#include +#include +#include +#include + +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + + +// TEST +#include + + + +// 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 + void generate_unit_points (const unsigned int, + vector > &); + + template <> + void generate_unit_points (const unsigned int, + vector > &) + { + Assert(false, ExcNotImplemented()); + } + + + template <> + void generate_unit_points (const unsigned int k, + vector > &p) + { + Assert(k<=4, ExcNotImplemented()); + Assert(p.size()==start_index2d[k+1]-start_index2d[k], ExcInternalError()); + for (unsigned int i=0; i + void generate_unit_points (const unsigned int k, + vector > &p) + { + Assert(k<=2, ExcNotImplemented()); + Assert(p.size()==start_index3d[k+1]-start_index3d[k], ExcInternalError()); + for (unsigned int i=0; i +FE_DGPMonomial::FE_DGPMonomial (const unsigned int degree) + : + FiniteElement (FiniteElementData(get_dpo_vector(degree), 1, degree), + std::vector(FiniteElementData(get_dpo_vector(degree), 1, degree).dofs_per_cell,true), + std::vector >(FiniteElementData(get_dpo_vector(degree), 1, degree).dofs_per_cell, + std::vector(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 +std::string +FE_DGPMonomial::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 +FiniteElement * +FE_DGPMonomial::clone() const +{ + return new FE_DGPMonomial(get_degree()); +} + + + +template +double +FE_DGPMonomial::shape_value (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + return polynomial_space.compute_value(i, p); +} + + + +template +double +FE_DGPMonomial::shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + Assert (component == 0, ExcIndexRange (component, 0, 1)); + return polynomial_space.compute_value(i, p); +} + + + +template +Tensor<1,dim> +FE_DGPMonomial::shape_grad (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + return polynomial_space.compute_grad(i, p); +} + + +template +Tensor<1,dim> +FE_DGPMonomial::shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + Assert (component == 0, ExcIndexRange (component, 0, 1)); + return polynomial_space.compute_grad(i, p); +} + + + +template +Tensor<2,dim> +FE_DGPMonomial::shape_grad_grad (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); + return polynomial_space.compute_grad_grad(i, p); +} + + + +template +Tensor<2,dim> +FE_DGPMonomial::shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (idofs_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 +void +FE_DGPMonomial:: +get_interpolation_matrix (const FiniteElementBase &x_source_fe, + FullMatrix &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*>(&x_source_fe) != 0), + typename FiniteElementBase:: + ExcInterpolationNotImplemented()); + + // ok, source is a Q element, so + // we will be able to do the work + const FE_DGPMonomial &source_fe + = dynamic_cast&>(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() +void +FE_DGPMonomial::initialize_embedding () +{ + vector > unit_points(this->dofs_per_cell); + generate_unit_points(get_degree(), unit_points); + + FullMatrix cell_interpolation (this->dofs_per_cell, + this->dofs_per_cell); + FullMatrix subcell_interpolation (this->dofs_per_cell, + this->dofs_per_cell); + for (unsigned int child=0; child::children_per_cell; ++child) + this->prolongation[child].reinit (this->dofs_per_cell, + this->dofs_per_cell); + for (unsigned int child=0; child::children_per_cell; ++child) + { + for (unsigned int j=0; jdofs_per_cell; ++j) + { + const Point &p_subcell=unit_points[j]; + + const Point p_cell = + GeometryInfo::child_to_cell_coordinates (p_subcell, child); + + for (unsigned int i=0; idofs_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; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + if (std::fabs(this->prolongation[child](i,j)) < 2e-14*get_degree()*dim) + this->prolongation[child](i,j) = 0.; + } +} + + +template +void +FE_DGPMonomial::initialize_restriction () +{ + Assert(false, ExcNotImplemented()); +} + + +//---------------------------------------------------------------------- +// Auxiliary functions +//---------------------------------------------------------------------- + + +template +std::vector +FE_DGPMonomial::get_dpo_vector(unsigned int deg) +{ + std::vector dpo(dim+1, 0U); + dpo[dim] = ++deg; + for (unsigned int i=1;i +UpdateFlags +FE_DGPMonomial::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 +UpdateFlags +FE_DGPMonomial::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 +typename Mapping::InternalDataBase * +FE_DGPMonomial::get_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &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 values(0); + std::vector > grads(0); + std::vector > 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; idofs_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 +void +FE_DGPMonomial::fill_fe_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &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 (fedata); + + const UpdateFlags flags(fe_data.current_update_flags()); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + if (flags & update_values) + for (unsigned int i=0;icompute_2nd (mapping, cell, + QProjector::DataSetDescriptor::cell(), + mapping_data, fe_data, data); +} + + + +template +void +FE_DGPMonomial::fill_fe_face_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &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 (fedata); + + // offset determines which data set + // to take (all data sets for all + // faces are stored contiguously) + const typename QProjector::DataSetDescriptor offset + = (QProjector::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; kdofs_per_cell; ++k) + { + for (unsigned int i=0;icompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +template +void +FE_DGPMonomial::fill_fe_subface_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const unsigned int subface, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &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 (fedata); + + // offset determines which data set + // to take (all data sets for all + // sub-faces are stored contiguously) + const typename QProjector::DataSetDescriptor offset + = (QProjector::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; kdofs_per_cell; ++k) + { + for (unsigned int i=0;icompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +template +unsigned int +FE_DGPMonomial::n_base_elements () const +{ + return 1; +} + + + +template +const FiniteElement & +FE_DGPMonomial::base_element (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return *this; +} + + + +template +unsigned int +FE_DGPMonomial::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 +unsigned int +FE_DGPMonomial::memory_consumption () const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + + +template class FE_DGPMonomial;