From d2600fc1382130d4e907a104a46779eccba32f33 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Wed, 25 Nov 2009 05:22:16 +0000 Subject: [PATCH] add FE_FaceQ and FE_PolyFace for testing git-svn-id: https://svn.dealii.org/trunk@20164 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_dgq.h | 5 +- deal.II/deal.II/include/fe/fe_face.h | 86 ++++ deal.II/deal.II/include/fe/fe_poly.h | 9 +- deal.II/deal.II/include/fe/fe_poly_face.h | 387 ++++++++++++++++++ .../include/fe/fe_poly_face.templates.h | 281 +++++++++++++ deal.II/deal.II/source/fe/fe_face.cc | 78 ++++ 6 files changed, 840 insertions(+), 6 deletions(-) create mode 100644 deal.II/deal.II/include/fe/fe_face.h create mode 100644 deal.II/deal.II/include/fe/fe_poly_face.h create mode 100644 deal.II/deal.II/include/fe/fe_poly_face.templates.h create mode 100644 deal.II/deal.II/source/fe/fe_face.cc diff --git a/deal.II/deal.II/include/fe/fe_dgq.h b/deal.II/deal.II/include/fe/fe_dgq.h index a50375e7af..e0b44726b1 100644 --- a/deal.II/deal.II/include/fe/fe_dgq.h +++ b/deal.II/deal.II/include/fe/fe_dgq.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -13,9 +13,6 @@ #ifndef __deal2__fe_dgq_h #define __deal2__fe_dgq_h - -//TODO: we should probably turn the derivation around, i.e. make FE_DGQ a special case of FE_DGQ_Generic/Arbitrary, and then have other classes that only have the constructor and clone() that use other sets of quadrature nodes. - #include #include #include diff --git a/deal.II/deal.II/include/fe/fe_face.h b/deal.II/deal.II/include/fe/fe_face.h new file mode 100644 index 0000000000..97cab1f999 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_face.h @@ -0,0 +1,86 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 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_face_h +#define __deal2__fe_face_h + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +/** + * @warning This class has not been tested + * + * A finite element, which is a tensor product polynomial on each face + * and undefined in the interior of the cells. + * + * This finite element is the trace space of FE_RavirtThomas on the + * faces and serves in hybridized methods. + * + * @author Guido Kanschat, 2009 + */ +template +class FE_FaceQ : FE_PolyFace, dim, spacedim> +{ + public: + /** + * Constructor for tensor product + * polynomials of degree + * p. The shape + * functions created using this + * constructor correspond to + * Legendre polynomials in each + * coordinate direction. + */ + FE_FaceQ(unsigned int p); + + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * FE_DGQ(degree) , with + * dim and degree + * replaced by appropriate + * values. + */ + virtual std::string get_name () 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 + * FiniteElement + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const; + + private: + /** + * Return vector with dofs per + * vertex, line, quad, hex. + */ + static std::vector get_dpo_vector (const unsigned int deg); +}; + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/fe/fe_poly.h b/deal.II/deal.II/include/fe/fe_poly.h index ea3d772b96..cc0875e653 100644 --- a/deal.II/deal.II/include/fe/fe_poly.h +++ b/deal.II/deal.II/include/fe/fe_poly.h @@ -26,10 +26,12 @@ DEAL_II_NAMESPACE_OPEN * FiniteElement classes based on a polynomial spaces like the * TensorProductPolynomials or a PolynomialSpace classes. * - * Every class that implements following functions can be used as + * Every class conforming to the following interface can be used as * template parameter POLY. * * @code + * static const unsigned int dimension; + * * double compute_value (const unsigned int i, * const Point &p) const; * @@ -55,7 +57,10 @@ DEAL_II_NAMESPACE_OPEN * functions depend on the cells in real space, the update_once() and * update_each() functions must be overloaded. * - * @author Ralf Hartmann 2004 + * @todo Since nearly all functions for spacedim != dim are + * specialized, this class needs cleaning up. + * + * @author Ralf Hartmann 2004, Guido Kanschat, 2009 **/ template class FE_Poly : public FiniteElement diff --git a/deal.II/deal.II/include/fe/fe_poly_face.h b/deal.II/deal.II/include/fe/fe_poly_face.h new file mode 100644 index 0000000000..5f5e96e5fe --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_poly_face.h @@ -0,0 +1,387 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 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_face_h +#define __deal2__fe_poly_face_h + + +#include + +DEAL_II_NAMESPACE_OPEN + +/*!@addtogroup febase */ +/*@{*/ + +/** + * @warning This class is not sufficiently tested yet! + * + * This class gives a unified framework for the implementation of + * FiniteElement classes only located on faces of the mesh. Theu are + * based on polynomial spaces like the TensorProductPolynomials or a + * PolynomialSpace classes. + * + * Every class that implements following functions can be used as + * template parameter POLY. + * + * @code + * double compute_value (const unsigned int i, + * const Point &p) const; + * @endcode + * Example classes are TensorProductPolynomials, PolynomialSpace or + * PolynomialsP. + * + * This class is not a fully implemented FiniteElement class. Instead + * there are several pure virtual functions declared in the + * FiniteElement and FiniteElement classes which cannot + * implemented by this class but are left for implementation in + * derived classes. + * + * Furthermore, this class assumes that shape functions of the + * FiniteElement under consideration do not depend on the + * actual shape of the cells in real space, i.e. update_once() + * includes update_values. For FiniteElements whose shape + * functions depend on the cells in real space, the update_once() and + * update_each() functions must be overloaded. + * + * @author Guido Kanschat, 2009 + **/ +template +class FE_PolyFace : public FiniteElement +{ + public: + /** + * Constructor. + */ + FE_PolyFace (const POLY& poly_space, + const FiniteElementData &fe_data, + const std::vector &restriction_is_additive_flags); + + /** + * Return the polynomial degree + * of this finite element, + * i.e. the value passed to the + * constructor. + */ + unsigned int get_degree () const; + + /** + * Return the value of the + * ith shape function at + * the point p. See the + * FiniteElement 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 + * componentth vector + * component of the ith + * shape function at the point + * p. See the + * FiniteElement 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 + * _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 + * ith shape function at + * the point p. See the + * FiniteElement 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 + * componentth vector + * component of the ith + * shape function at the point + * p. See the + * FiniteElement 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 + * _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 + * ith shape function at + * point p on the unit + * cell. See the + * FiniteElement 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 componentth + * vector component of the + * ith shape function at + * the point p. See the + * FiniteElement 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 + * _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, + * base_element(0) is + * this, and all other + * indices throw an error. + */ + virtual const FiniteElement & + base_element (const unsigned int index) const; + + /** + * Multiplicity of base element + * index. Since this is + * a scalar element, + * element_multiplicity(0) + * returns one, and all other + * indices will throw an error. + */ + virtual unsigned int element_multiplicity (const unsigned int index) const; + + + protected: + + virtual + typename Mapping::InternalDataBase * + get_data (const UpdateFlags, + const Mapping& mapping, + const Quadrature& quadrature) const ; + + typename Mapping::InternalDataBase * + get_face_data (const UpdateFlags, + const Mapping& mapping, + const Quadrature& quadrature) const ; + + typename Mapping::InternalDataBase * + get_subface_data (const UpdateFlags, + const Mapping& mapping, + const Quadrature& quadrature) const ; + + virtual void + fill_fe_values (const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData &data, + CellSimilarity::Similarity &cell_similarity) const; + + virtual void + fill_fe_face_values (const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const ; + + virtual void + fill_fe_subface_values (const Mapping &mapping, + const typename Triangulation::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 + * flags. + * + * For the purpuse of this + * function, refer to the + * documentation in + * FiniteElement. + * + * This class assumes that shape + * functions of this + * FiniteElement do not + * depend on the actual shape of + * the cells in real + * space. Therefore, the effect + * in this element is as follows: + * if update_values is + * set in 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 + * flags. + * + * For the purpuse of this + * function, refer to the + * documentation in + * FiniteElement. + * + * This class assumes that shape + * functions of this + * FiniteElement do not + * depend on the actual shape of + * the cells in real + * space. + * + * The effect in this element is + * as follows: + *
    + + *
  • if + * update_gradients is + * set, the result will contain + * update_gradients and + * 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 + * update_covariant_transformation + * is actually performed by the + * Mapping object used in + * conjunction with this finite + * element.
  • if + * update_hessians + * is set, the result will + * contain + * update_hessians + * and + * 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. + * + *
+ */ + 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 FiniteElement::InternalDataBase + { + public: + /** + * Array with shape function + * values in quadrature + * points on one face. 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 one face of 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. + * + * In particular, we can + * simply copy the same set + * of values to each of the + * faces. + */ + std::vector > shape_values; + }; + + /** + * The polynomial space. Its type + * is given by the template + * parameter POLY. + */ + POLY poly_space; +}; + +/*@}*/ + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/fe/fe_poly_face.templates.h b/deal.II/deal.II/include/fe/fe_poly_face.templates.h new file mode 100644 index 0000000000..16aaee25cb --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_poly_face.templates.h @@ -0,0 +1,281 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 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 + + +DEAL_II_NAMESPACE_OPEN + +template +FE_PolyFace::FE_PolyFace ( + const POLY& poly_space, + const FiniteElementData &fe_data, + const std::vector &restriction_is_additive_flags): + FiniteElement (fe_data, + restriction_is_additive_flags, + std::vector > (1, std::vector(1,true))), + poly_space(poly_space) +{ + AssertDimension(dim, POLY::dimension+1); +} + + +template +unsigned int +FE_PolyFace::get_degree () const +{ + return this->degree; +} + + +//--------------------------------------------------------------------------- +// Auxiliary functions +//--------------------------------------------------------------------------- + + + + +template +UpdateFlags +FE_PolyFace::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_PolyFace::update_each (const UpdateFlags flags) const +{ + UpdateFlags out = update_default; + + if (flags & update_gradients) + out |= update_gradients | update_covariant_transformation; + if (flags & update_hessians) + out |= update_hessians | update_covariant_transformation; + if (flags & update_cell_normal_vectors) + out |= update_cell_normal_vectors | update_JxW_values; + + return out; +} + + + +//--------------------------------------------------------------------------- +// Data field initialization +//--------------------------------------------------------------------------- + +template +typename Mapping::InternalDataBase * +FE_PolyFace::get_data ( + const UpdateFlags, + const Mapping&, + const Quadrature&) const +{ + Assert(false, ExcNotImplemented()); + return 0; +} + + +template +typename Mapping::InternalDataBase * +FE_PolyFace::get_face_data ( + const UpdateFlags update_flags, + const 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.size(); + + // 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.resize (this->dofs_per_cell, + std::vector (n_q_points)); + for (unsigned int i=0; idofs_per_cell; ++k) + data->shape_values[k][i] = values[k]; + } + } + // No derivatives of this element + // are implemented. + if (flags & update_gradients || flags & update_hessians) + { + Assert(false, ExcNotImplemented()); + } + + return data; +} + + + +template +typename Mapping::InternalDataBase * +FE_PolyFace::get_subface_data ( + const UpdateFlags flags, + const Mapping& mapping, + const Quadrature& quadrature) const +{ + return get_face_data (flags, mapping, + QProjector::project_to_all_children(quadrature)); +} + + + + + +//--------------------------------------------------------------------------- +// Fill data of FEValues +//--------------------------------------------------------------------------- +template +void +FE_PolyFace::fill_fe_values + (const Mapping&, + const typename Triangulation::cell_iterator&, + const Quadrature&, + typename Mapping::InternalDataBase&, + typename Mapping::InternalDataBase&, + FEValuesData&, + CellSimilarity::Similarity&) const +{ + Assert(false, ExcNotImplemented()); +} + + + +template +void +FE_PolyFace::fill_fe_face_values ( + const Mapping&, + const typename Triangulation::cell_iterator&, + const unsigned int, + const Quadrature& quadrature, + typename Mapping::InternalDataBase&, + 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 + Assert (dynamic_cast (&fedata) != 0, ExcInternalError()); + InternalData &fe_data = static_cast (fedata); + + 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; i +void +FE_PolyFace::fill_fe_subface_values ( + const Mapping&, + const typename Triangulation::cell_iterator&, + const unsigned int, + const unsigned int subface, + const Quadrature& quadrature, + typename Mapping::InternalDataBase&, + 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 + Assert (dynamic_cast (&fedata) != 0, ExcInternalError()); + InternalData &fe_data = static_cast (fedata); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + unsigned int offset = subface*quadrature.size(); + for (unsigned int k=0; kdofs_per_cell; ++k) + { + if (flags & update_values) + for (unsigned int i=0; i +unsigned int +FE_PolyFace::n_base_elements () const +{ + return 1; +} + + + +template +const FiniteElement & +FE_PolyFace::base_element (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return *this; +} + + + +template +unsigned int +FE_PolyFace::element_multiplicity (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return 1; +} + + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/fe/fe_face.cc b/deal.II/deal.II/source/fe/fe_face.cc new file mode 100644 index 0000000000..43d966b5a9 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_face.cc @@ -0,0 +1,78 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002, 2003, 2004, 2005, 2006, 2007, 2008 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 + +DEAL_II_NAMESPACE_OPEN + +template +FE_FaceQ::FE_FaceQ (const unsigned int degree) + : + FE_PolyFace, dim, spacedim> ( + TensorProductPolynomials(Polynomials::Legendre::generate_complete_basis(degree)), + FiniteElementData(get_dpo_vector(degree), 1, degree, FiniteElementData::L2), + std::vector(1,true)) +{} + + +template +std::string +FE_FaceQ::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 + + std::ostringstream namebuf; + namebuf << "FE_FaceQ<" << dim << ">(" << this->degree << ")"; + + return namebuf.str(); +} + +template +bool +FE_FaceQ::has_support_on_face ( + const unsigned int shape_index, + const unsigned int face_index) const +{ + return (face_index == (shape_index/this->dofs_per_face)); +} + + +template +std::vector +FE_FaceQ::get_dpo_vector (const unsigned int deg) +{ + std::vector dpo(dim+1, 0U); + dpo[dim-1] = deg+1; + for (unsigned int i=1;i 1 + template class FE_PolyFace >; +//template class FE_PolyFace, deal_II_dimension>; +//template class FE_PolyFace, deal_II_dimension>; + template class FE_FaceQ; +#endif + +DEAL_II_NAMESPACE_CLOSE -- 2.39.5