From d1ce06e6fb856cc842b7b2bb3d371258028d0dc8 Mon Sep 17 00:00:00 2001 From: guido Date: Thu, 19 Sep 2002 22:21:54 +0000 Subject: [PATCH] new class FE_DGPNonparametric git-svn-id: https://svn.dealii.org/trunk@6488 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/fe/fe_dgp_nonparametric.h | 457 ++++++++++++++++++ .../deal.II/source/fe/fe_dgp_nonparametric.cc | 439 +++++++++++++++++ deal.II/doc/news/2002/c-3-4.html | 10 + 3 files changed, 906 insertions(+) create mode 100644 deal.II/deal.II/include/fe/fe_dgp_nonparametric.h create mode 100644 deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc diff --git a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h new file mode 100644 index 0000000000..b755739d1b --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h @@ -0,0 +1,457 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002 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_nonparametric_h +#define __deal2__fe_dgp_nonparametric_h + +#include +#include +#include +#include + +template class PolynomialSpace; +template class MappingQ; + + +/** + * Discontinuous finite elements evaluated at the mapped quadrature points. + * + * This finite element implements complete polynomial spaces, that is, + * $d$-dimensional polynomials of order $k$. + * + * The polynomials are not mapped. Therefore, they are constant, + * linear, quadratic, etc. on any grid cell. + * + * Since the polynomials are evaluated at the quadrature points of the + * actual grid cell, no grid transfer and interpolation matrices are + * available. + * + * The purpose of this class is experimental, therefore the + * implementation will remain incomplete. + * + * @author Guido Kanschat, 2002 + */ +template +class FE_DGPNonparametric : public FiniteElement +{ + public: + /** + * Constructor for tensor product + * polynomials of degree @p{k}. + */ + FE_DGPNonparametric (const unsigned int k); + + /** + * 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; + + /** + * 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; + + /** + * 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; + + + /** + * Declare a nested class which + * will hold static definitions of + * various matrices such as + * constraint and embedding + * matrices. The definition of + * the various static fields are + * in the files @p{fe_dgp_[123]d.cc} + * in the source directory. + */ + struct Matrices + { + /** + * Pointers to the embedding + * matrices, one for each + * polynomial degree starting + * from constant elements + */ + static const double * const embedding[][GeometryInfo::children_per_cell]; + + /** + * Number of elements (first + * index) the above field + * has. Equals the highest + * polynomial degree plus one + * for which the embedding + * matrices have been + * computed. + */ + static const unsigned int n_embedding_matrices; + + /** + * As @p{embedding} but for + * projection matrices. + */ + static const double * const projection_matrices[][GeometryInfo::children_per_cell]; + + /** + * As + * @p{n_embedding_matrices} + * but for projection + * matrices. + */ + static const unsigned int n_projection_matrices; + }; + + 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); + + /** + * 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: + // have some scratch arrays + std::vector values; + std::vector > grads; + std::vector > grad_grads; + }; + + /** + * Allow access from other dimensions. + */ + template friend class FE_DGPNonparametric; + + /** + * Allows @p{MappingQ} class to + * access to build_renumbering + * function. + */ + friend class MappingQ; +}; + + +// declaration of explicit specializations of member variables, if the +// compiler allows us to do that (the standard says we must) +#ifndef DEAL_II_MEMBER_VAR_SPECIALIZATION_BUG +template <> +const double * const FE_DGPNonparametric<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell]; + +template <> +const unsigned int FE_DGPNonparametric<1>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_DGPNonparametric<1>::Matrices::projection_matrices[][GeometryInfo<1>::children_per_cell]; + +template <> +const unsigned int FE_DGPNonparametric<1>::Matrices::n_projection_matrices; + +template <> +const double * const FE_DGPNonparametric<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell]; + +template <> +const unsigned int FE_DGPNonparametric<2>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_DGPNonparametric<2>::Matrices::projection_matrices[][GeometryInfo<2>::children_per_cell]; + +template <> +const unsigned int FE_DGPNonparametric<2>::Matrices::n_projection_matrices; + +template <> +const double * const FE_DGPNonparametric<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell]; + +template <> +const unsigned int FE_DGPNonparametric<3>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_DGPNonparametric<3>::Matrices::projection_matrices[][GeometryInfo<3>::children_per_cell]; + +template <> +const unsigned int FE_DGPNonparametric<3>::Matrices::n_projection_matrices; +#endif + +#endif diff --git a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc new file mode 100644 index 0000000000..0d18b01271 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc @@ -0,0 +1,439 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002 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 + + + +template +FE_DGPNonparametric::FE_DGPNonparametric (const unsigned int degree) + : + FiniteElement (FiniteElementData(get_dpo_vector(degree),1), + std::vector(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell,true), + std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + std::vector(1,true))), + degree(degree), + polynomial_space (Legendre::generate_complete_basis(degree)) +{ + // if defined, copy over matrices + // from precomputed arrays +// if ((degree < Matrices::n_embedding_matrices) && +// (Matrices::embedding[degree][0] != 0)) +// for (unsigned int c=0; c::children_per_cell; ++c) +// this->prolongation[c].fill (Matrices::embedding[degree][c]); +// else +// for (unsigned int i=0; i::children_per_cell;++i) +// this->prolongation[i].reinit(0,0); + + // restriction can be defined + // through projection for + // discontinuous elements, but is + // presently not implemented for DGPNonparametric + // elements. + // + // if it were, then the following + // snippet would be the right code +// if ((degree < Matrices::n_projection_matrices) && +// (Matrices::projection_matrices[degree] != 0)) +// { +// restriction[0].fill (Matrices::projection_matrices[degree]); +// } +// else +// // matrix undefined, set size to zero +// for (unsigned int i=0;i::children_per_cell;++i) +// restriction[i].reinit(0, 0); + // since not implemented, set to + // "empty" + for (unsigned int i=0;i::children_per_cell;++i) + restriction[i].reinit(0, 0); + + // note further, that these + // elements have neither support + // nor face-support points, so + // leave these fields empty +}; + + + +template +FiniteElement * +FE_DGPNonparametric::clone() const +{ + return new FE_DGPNonparametric(degree); +} + + + +template +double +FE_DGPNonparametric::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_DGPNonparametric::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_DGPNonparametric::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_DGPNonparametric::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_DGPNonparametric::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_DGPNonparametric::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 +//---------------------------------------------------------------------- + + +template +std::vector +FE_DGPNonparametric::get_dpo_vector(unsigned int deg) +{ + std::vector dpo(dim+1, static_cast(0)); + dpo[dim] = ++deg; + for (unsigned int i=1;i +UpdateFlags +FE_DGPNonparametric::update_once (const UpdateFlags) 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; +} + + +template +UpdateFlags +FE_DGPNonparametric::update_each (const UpdateFlags flags) const +{ + UpdateFlags out = flags & update_values; + + 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_DGPNonparametric::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); + + // initialize fields only if really + // necessary. otherwise, don't + // allocate memory + if (flags & update_values) + { + data->values.resize (this->dofs_per_cell); + } + + if (flags & update_gradients) + { + data->grads.resize (this->dofs_per_cell); + } + + // 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); + return data; +} + + + +//---------------------------------------------------------------------- +// Fill data of FEValues +//---------------------------------------------------------------------- + +template +void +FE_DGPNonparametric::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()); + + const unsigned int n_q_points = data.quadrature_points.size(); + + if (flags & (update_values | update_gradients)) + for (unsigned int i=0; idofs_per_cell; ++k) + { + if (flags & update_values) + data.shape_values[k][i] = fe_data.values[k]; + if (flags & update_gradients) + data.shape_gradients[k][i] = fe_data.grads[k]; + } + } + + if (flags & update_second_derivatives) + compute_2nd (mapping, cell, 0, mapping_data, fe_data, data); + + fe_data.first_cell = false; +} + + + +template +void +FE_DGPNonparametric::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); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + const unsigned int n_q_points = data.quadrature_points.size(); + + if (flags & (update_values | update_gradients)) + for (unsigned int i=0; idofs_per_cell; ++k) + { + if (flags & update_values) + data.shape_values[k][i] = fe_data.values[k]; + if (flags & update_gradients) + data.shape_gradients[k][i] = fe_data.grads[k]; + } + } + + // offset determines which data set + // to take (all data sets for all + // faces are stored contiguously) + const unsigned int offset = face * quadrature.n_quadrature_points; + + if (flags & update_second_derivatives) + compute_2nd (mapping, cell, offset, mapping_data, fe_data, data); + + fe_data.first_cell = false; +} + + + +template +void +FE_DGPNonparametric::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); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + const unsigned int n_q_points = data.quadrature_points.size(); + + if (flags & (update_values | update_gradients)) + for (unsigned int i=0; idofs_per_cell; ++k) + { + if (flags & update_values) + data.shape_values[k][i] = fe_data.values[k]; + if (flags & update_gradients) + data.shape_gradients[k][i] = fe_data.grads[k]; + } + } + + // offset determines which data set + // to take (all data sets for all + // sub-faces are stored contiguously) + const unsigned int offset = (face * GeometryInfo::subfaces_per_face + subface) + * quadrature.n_quadrature_points; + + if (flags & update_second_derivatives) + compute_2nd (mapping, cell, offset, mapping_data, fe_data, data); + + fe_data.first_cell = false; +} + + + +template +unsigned int +FE_DGPNonparametric::n_base_elements () const +{ + return 1; +}; + + + +template +const FiniteElement & +FE_DGPNonparametric::base_element (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return *this; +}; + + + +template +bool +FE_DGPNonparametric::has_support_on_face (const unsigned int, + const unsigned int) const +{ + return true; +} + + + +template +unsigned int +FE_DGPNonparametric::memory_consumption () const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + + +template +unsigned int +FE_DGPNonparametric::get_degree () const +{ + return degree; +}; + + + + +template class FE_DGPNonparametric; diff --git a/deal.II/doc/news/2002/c-3-4.html b/deal.II/doc/news/2002/c-3-4.html index 4c0927a9e0..b6b2595ab4 100644 --- a/deal.II/doc/news/2002/c-3-4.html +++ b/deal.II/doc/news/2002/c-3-4.html @@ -299,6 +299,16 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

deal.II

    +
  1. New: the class FE_DGPNonparametric implements finite elements + where shape functions are polynomials of order k on the + actual grid cell. This is achieved by evaluating the polynomials at + the mapped quadrature points. No grid transfer matrices are + available for this class. +
    + (GK 09/19/2002) +

    +
  2. Fixed: Some of the various instances of the VectorTools::interpolate_boundary_values functions -- 2.39.5