From: hartmann Date: Tue, 16 Mar 2004 16:35:30 +0000 (+0000) Subject: New (intermediate) FiniteElement class which allows implementation of FiniteElements... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=74bc0ff48bb41e058a85a4e4c84238ab3638ab3f;p=dealii-svn.git New (intermediate) FiniteElement class which allows implementation of FiniteElements based on PolynomialSpace and TensorProductPolynomials in a unified way. git-svn-id: https://svn.dealii.org/trunk@8775 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_poly.h b/deal.II/deal.II/include/fe/fe_poly.h new file mode 100644 index 0000000000..1bdb3860ee --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_poly.h @@ -0,0 +1,384 @@ +//--------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------- +#ifndef __deal2__fe_poly_h +#define __deal2__fe_poly_h + + +#include + + +/** + * This class gives a unified framework for implementing a + * FiniteElement class based on a polynomial space like a + * @p TensorProductSpace or a @p PolynomialSpace. + * + * Every class that implements following functions can be used as + * template parameter @p POLY. Example classes are @p + * TensorProductSpace and @p PolynomialSpace. + * + * @code + * double compute_value (const unsigned int i, + * const Point &p) const; + * + * Tensor<1,dim> compute_grad (const unsigned int i, + * const Point &p) const; + * + * Tensor<2,dim> compute_grad_grad (const unsigned int i, + * const Point &p) const; + * @endcode + * + * This class is not a fully implemented FiniteElement class. Instead + * there are several pure virtual functions declared in the + * FiniteElement and FiniteElementBase classes which cannot + * implemented by this class but are left for implementation in + * derived classes. + * + * Todos: + * - checke dim of POLY + * - templatisiere nur auf POLY, dim ergibt sich durch POLY::dim + * + * @author Ralf Hartmann 2004 + **/ +template +class FE_Poly : public FiniteElement +{ + public: + /** + * Constructor. + */ + FE_Poly (const POLY& poly_space, + const FiniteElementData &fe_data, + const std::vector &restriction_is_additive_flags, + const std::vector > &nonzero_components); + + /** + * Return the value of the + * @p{i}th shape function at the + * point @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + */ + virtual double shape_value (const unsigned int i, + const Point &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; + + /** + * 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; + + + protected: + + /** + * 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 ; + + + /** + * Determine the values that need + * to be computed on the unit + * cell to be able to compute all + * values required by @p{flags}. + * + * For the purpuse of this + * function, refer to the + * documentation in + * @p{FiniteElement}. + * + * The effect in this element is + * as follows: if + * @p{update_values} is set in + * @p{flags}, copy it to the + * result. All other flags of the + * result are cleared, since + * everything else must be + * computed for each cell. + */ + virtual UpdateFlags update_once (const UpdateFlags flags) const; + + /** + * Determine the values that need + * to be computed on every + * cell to be able to compute all + * values required by @p{flags}. + * + * For the purpuse of this + * function, refer to the + * documentation in + * @p{FiniteElement}. + * + * The effect in this element is + * as follows: + * @begin{itemize} + * @item if @p{update_gradients} + * is set, the result will + * contain @p{update_gradients} + * and + * @p{update_covariant_transformation}. + * The latter is required to + * transform the gradient on the + * unit cell to the real + * cell. Remark, that the action + * required by + * @p{update_covariant_transformation} + * is actually performed by the + * @p{Mapping} object used in + * conjunction with this finite + * element. + * @item if + * @p{update_second_derivatives} + * is set, the result will + * contain + * @p{update_second_derivatives} + * and + * @p{update_covariant_transformation}. + * The rationale is the same as + * above and no higher + * derivatives of the + * transformation are required, + * since we use difference + * quotients for the actual + * computation. + * @end{itemize} + */ + virtual UpdateFlags update_each (const UpdateFlags flags) const; + + + /** + * Fields of cell-independent data. + * + * For information about the + * general purpose of this class, + * see the documentation of the + * base class. + */ + class InternalData : public FiniteElementBase::InternalDataBase + { + public: + /** + * Array with shape function + * values in quadrature + * points. There is one + * row for each shape + * function, containing + * values for each quadrature + * point. + * + * In this array, we store + * the values of the shape + * function in the quadrature + * points on the unit + * cell. Since these values + * do not change under + * transformation to the real + * cell, we only need to copy + * them over when visiting a + * concrete cell. + */ + Table<2,double> shape_values; + + /** + * Array with shape function + * gradients in quadrature + * points. There is one + * row for each shape + * function, containing + * values for each quadrature + * point. + * + * We store the gradients in + * the quadrature points on + * the unit cell. We then + * only have to apply the + * transformation (which is a + * matrix-vector + * multiplication) when + * visiting an actual cell. + */ + Table<2,Tensor<1,dim> > shape_gradients; + }; + + /** + * The polynomial space. + */ + POLY poly_space; +}; + + + +#endif diff --git a/deal.II/deal.II/source/fe/fe_poly.cc b/deal.II/deal.II/source/fe/fe_poly.cc new file mode 100644 index 0000000000..8bde0740bb --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_poly.cc @@ -0,0 +1,419 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------- + +// #include +// #include +// #include +// #include +// #include +// #include +#include +#include +#include +#include +//#include +// #include +// #include +#include + +// #ifdef HAVE_STD_STRINGSTREAM +// # include +// #else +// # include +// #endif + + +template +FE_Poly::FE_Poly (const POLY& poly_space, + const FiniteElementData &fe_data, + const std::vector &restriction_is_additive_flags, + const std::vector > &nonzero_components): + FiniteElement (fe_data, + restriction_is_additive_flags, + nonzero_components), + poly_space(poly_space) +{} + + + +template +double +FE_Poly::shape_value (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + return poly_space.compute_value(i, p); +} + + +template +double +FE_Poly::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 poly_space.compute_value(i, p); +} + + + +template +Tensor<1,dim> +FE_Poly::shape_grad (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + return poly_space.compute_grad(i, p); +} + + + +template +Tensor<1,dim> +FE_Poly::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 poly_space.compute_grad(i, p); +} + + + +template +Tensor<2,dim> +FE_Poly::shape_grad_grad (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + return poly_space.compute_grad_grad(i, p); +} + + + +template +Tensor<2,dim> +FE_Poly::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 poly_space.compute_grad_grad(i, p); +} + + + +//---------------------------------------------------------------------- +// Auxiliary functions +//---------------------------------------------------------------------- + + + + +template +UpdateFlags +FE_Poly::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_Poly::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_Poly::get_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature) const +{ + // generate a new data object and + // initialize some fields + InternalData* data = new InternalData; + + // check what needs to be + // initialized only once and what + // on every cell/face/subface we + // visit + data->update_once = update_once(update_flags); + data->update_each = update_each(update_flags); + data->update_flags = data->update_once | data->update_each; + + const UpdateFlags flags(data->update_flags); + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // some scratch arrays + std::vector 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) + data->shape_values[k][i] = values[k]; + + if (flags & update_gradients) + for (unsigned int k=0; kdofs_per_cell; ++k) + data->shape_gradients[k][i] = grads[k]; + } + return data; +} + + + + +//---------------------------------------------------------------------- +// Fill data of FEValues +//---------------------------------------------------------------------- + +template +void +FE_Poly::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; i::DataSetDescriptor dsd; + if (flags & update_second_derivatives) + this->compute_2nd (mapping, cell, dsd.cell(), + mapping_data, fe_data, data); +} + + + +template +void +FE_Poly::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 dsd; + const typename QProjector::DataSetDescriptor offset + = dsd.face (face, cell->face_orientation(face), + quadrature.n_quadrature_points); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + if (flags & update_values) + for (unsigned int i=0; icompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +template +void +FE_Poly::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 dsd; + const typename QProjector::DataSetDescriptor offset + = dsd.sub_face (face, subface, cell->face_orientation(face), + quadrature.n_quadrature_points); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + if (flags & update_values) + for (unsigned int i=0; icompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +template +unsigned int +FE_Poly::n_base_elements () const +{ + return 1; +} + + + +template +const FiniteElement & +FE_Poly::base_element (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return *this; +} + + + +template +unsigned int +FE_Poly::element_multiplicity (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return 1; +} + + + +template class FE_Poly, deal_II_dimension>; +template class FE_Poly, deal_II_dimension>; +template class FE_Poly, deal_II_dimension>;