From 806909b7ba17af30b03827734f79f68524b364d6 Mon Sep 17 00:00:00 2001 From: hartmann Date: Tue, 16 Mar 2004 16:36:47 +0000 Subject: [PATCH] Make use of the new FE_Poly class. git-svn-id: https://svn.dealii.org/trunk@8776 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_dgp.h | 300 +-------------- deal.II/deal.II/include/fe/fe_dgp_monomial.h | 312 +--------------- deal.II/deal.II/source/fe/fe_dgp.cc | 368 +----------------- deal.II/deal.II/source/fe/fe_dgp_monomial.cc | 374 +------------------ 4 files changed, 39 insertions(+), 1315 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_dgp.h b/deal.II/deal.II/include/fe/fe_dgp.h index 4839bbbbce..9d1d0b9591 100644 --- a/deal.II/deal.II/include/fe/fe_dgp.h +++ b/deal.II/deal.II/include/fe/fe_dgp.h @@ -14,15 +14,11 @@ #define __deal2__fe_dgp_h #include -#include #include -#include -#include +#include -template class PolynomialSpace; template class MappingQ; - /*!@addtogroup fe */ /*@{*/ @@ -30,21 +26,21 @@ template class MappingQ; * Discontinuous finite elements based on Legendre polynomials. * * This finite element implements complete polynomial spaces, that is, - * $d$-dimensional polynomials of order $k$. The underlying + * dim-dimensional polynomials of degree p. The underlying * polynomials form a Legendre basis on the unit square. Thus, the * mass matrix is diagonal, if the grid cells are parallelograms. * - * @author Guido Kanschat, 2001, 2002 + * @author Guido Kanschat, 2001, 2002, Ralf Hartmann 2004 */ template -class FE_DGP : public FiniteElement +class FE_DGP : public FE_Poly,dim> { public: /** * Constructor for tensor product - * polynomials of degree @p{k}. + * polynomials of degree @p p. */ - FE_DGP (const unsigned int k); + FE_DGP (const unsigned int p); /** * Return a string that uniquely @@ -57,109 +53,6 @@ class FE_DGP : public FiniteElement */ 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, @@ -167,34 +60,6 @@ class FE_DGP : public FiniteElement * constructor. */ unsigned int get_degree () const; - - /** - * Number of base elements in a - * mixed discretization. Since - * this is a scalar element, - * return one. - */ - virtual unsigned int n_base_elements () const; - - /** - * Access to base element - * objects. Since this element is - * scalar, @p{base_element(0)} is - * @p{this}, and all other - * indices throw an error. - */ - virtual const FiniteElement & - 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. @@ -282,59 +147,6 @@ class FE_DGP : public FiniteElement * 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: @@ -349,111 +161,11 @@ class FE_DGP : public FiniteElement * @p{FiniteElementData}. */ static std::vector get_dpo_vector(unsigned int degree); - - /** - * Given a set of flags indicating - * what quantities are requested - * from a @p{FEValues} object, - * return which of these can be - * precomputed once and for - * all. Often, the values of - * shape function at quadrature - * points can be precomputed, for - * example, in which case the - * return value of this function - * would be the logical and of - * the input @p{flags} and - * @p{update_values}. - * - * For the present kind of finite - * element, this is exactly the - * case. - */ - virtual UpdateFlags update_once (const UpdateFlags flags) const; - - /** - * This is the opposite to the - * above function: given a set of - * flags indicating what we want - * to know, return which of these - * need to be computed each time - * we visit a new cell. - * - * If for the computation of one - * quantity something else is - * also required (for example, we - * often need the covariant - * transformation when gradients - * need to be computed), include - * this in the result as well. - */ - virtual UpdateFlags update_each (const UpdateFlags flags) const; /** * Degree of the polynomials. */ const unsigned int degree; - - /** - * Pointer to an object - * representing the polynomial - * space used here. - */ - const PolynomialSpace 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; - }; /** * Allow access from other dimensions. diff --git a/deal.II/deal.II/include/fe/fe_dgp_monomial.h b/deal.II/deal.II/include/fe/fe_dgp_monomial.h index 7de1f4ad18..e6865e48a6 100644 --- a/deal.II/deal.II/include/fe/fe_dgp_monomial.h +++ b/deal.II/deal.II/include/fe/fe_dgp_monomial.h @@ -15,8 +15,7 @@ #include #include -#include -#include +#include template class MappingQ; @@ -26,138 +25,36 @@ template class MappingQ; /** * Discontinuous finite elements based on monomials of degree up to - * k + * p. * - * This finite element makes use of the @ref{PolynomialsP} class which - * implements $d$-dimensional polynomials of degree $k$ based the - * @ref{Polynomials::Monomial} and the @ref{PolynomialSpace} classes. + * This finite element makes use of the PolynomialsP class which + * implements dim-dimensional polynomials of degree + * p based the Polynomials::Polynomial and the + * PolynomialSpace classes. * * @author Ralf Hartmann, 2004 */ template -class FE_DGPMonomial : public FiniteElement +class FE_DGPMonomial : public FE_Poly,dim> { public: /** * Constructor for the polynomial - * space of degree @p{k}. + * space of degree p. */ - FE_DGPMonomial (const unsigned int k); + FE_DGPMonomial (const unsigned int p); /** * Return a string that uniquely * identifies a finite * element. This class returns - * @p{FE_DGPMonomial(degree)}, with - * @p{dim} and @p{degree} - * replaced by appropriate - * values. + * FE_DGPMonomial(degree), + * with dim and + * p 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, @@ -185,34 +82,6 @@ class FE_DGPMonomial : public FiniteElement 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. @@ -253,59 +122,6 @@ class FE_DGPMonomial : public FiniteElement * 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: @@ -334,106 +150,6 @@ class FE_DGPMonomial : public FiniteElement * 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 @@ -449,7 +165,7 @@ template inline unsigned int FE_DGPMonomial::get_degree () const { - return polynomial_space.degree(); + return this->poly_space.degree(); } diff --git a/deal.II/deal.II/source/fe/fe_dgp.cc b/deal.II/deal.II/source/fe/fe_dgp.cc index 9230f796cf..7a42f22e45 100644 --- a/deal.II/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/deal.II/source/fe/fe_dgp.cc @@ -2,7 +2,7 @@ // $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 @@ -30,12 +30,13 @@ template FE_DGP::FE_DGP (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))), - degree(degree), - polynomial_space (Polynomials::Legendre::generate_complete_basis(degree)) + FE_Poly, dim> ( + PolynomialSpace(Polynomials::Legendre::generate_complete_basis(degree)), + 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))), + degree(degree) { // if defined, copy over matrices // from precomputed arrays @@ -114,76 +115,6 @@ FE_DGP::clone() const -template -double -FE_DGP::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_DGP::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_DGP::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_DGP::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_DGP::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_DGP::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); -} - - //---------------------------------------------------------------------- // Auxiliary functions //---------------------------------------------------------------------- @@ -204,289 +135,6 @@ FE_DGP::get_dpo_vector(unsigned int deg) } -template -UpdateFlags -FE_DGP::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_DGP::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_DGP::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_DGP::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_DGP::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_DGP::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_DGP::n_base_elements () const -{ - return 1; -} - - - -template -const FiniteElement & -FE_DGP::base_element (const unsigned int index) const -{ - Assert (index==0, ExcIndexRange(index, 0, 1)); - return *this; -} - - - -template -unsigned int -FE_DGP::element_multiplicity (const unsigned int index) const -{ - Assert (index==0, ExcIndexRange(index, 0, 1)); - return 1; -} - - - template bool FE_DGP::has_support_on_face (const unsigned int, diff --git a/deal.II/deal.II/source/fe/fe_dgp_monomial.cc b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc index e9773a7817..9382a11fe1 100644 --- a/deal.II/deal.II/source/fe/fe_dgp_monomial.cc +++ b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc @@ -122,13 +122,14 @@ namespace template 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) + FE_Poly, dim> ( + PolynomialsP(degree), + 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))) { - Assert(polynomial_space.n()==this->dofs_per_cell, ExcInternalError()); + Assert(this->poly_space.n()==this->dofs_per_cell, ExcInternalError()); // DG doesn't have constraints, so // leave them empty @@ -182,77 +183,6 @@ FE_DGPMonomial::clone() const -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:: @@ -313,8 +243,8 @@ FE_DGPMonomial::initialize_embedding () 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); + cell_interpolation(j,i) = this->poly_space.compute_value (i, p_cell); + subcell_interpolation(j,i) = this->poly_space.compute_value (i, p_subcell); } } @@ -361,288 +291,6 @@ FE_DGPMonomial::get_dpo_vector(unsigned int deg) } -template -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 <> @@ -668,7 +316,7 @@ FE_DGPMonomial<2>::has_support_on_face (const unsigned int shape_index, else { unsigned int degrees[2]; - polynomial_space.directional_degrees(shape_index, degrees); + this->poly_space.directional_degrees(shape_index, degrees); if ((face_index==0 && degrees[1]==0) || (face_index==3 && degrees[0]==0)) support_on_face=true; @@ -691,7 +339,7 @@ FE_DGPMonomial<3>::has_support_on_face (const unsigned int shape_index, else { unsigned int degrees[3]; - polynomial_space.directional_degrees(shape_index, degrees); + this->poly_space.directional_degrees(shape_index, degrees); if ((face_index==0 && degrees[1]==0) || (face_index==2 && degrees[2]==0) || (face_index==5 && degrees[0]==0)) -- 2.39.5