From 8a04dacc4fd2fbaf8f78bc88ec0c34fb597c0e4b Mon Sep 17 00:00:00 2001 From: guido Date: Mon, 11 Feb 2002 12:58:57 +0000 Subject: [PATCH] new FE_DGP git-svn-id: https://svn.dealii.org/trunk@5501 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_dgp.h | 364 ++++++++++++++++++++++++ deal.II/deal.II/source/fe/fe_dgp.cc | 411 ++++++++++++++++++++++++++++ 2 files changed, 775 insertions(+) create mode 100644 deal.II/deal.II/include/fe/fe_dgp.h create mode 100644 deal.II/deal.II/source/fe/fe_dgp.cc diff --git a/deal.II/deal.II/include/fe/fe_dgp.h b/deal.II/deal.II/include/fe/fe_dgp.h new file mode 100644 index 0000000000..787ea5fb30 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_dgp.h @@ -0,0 +1,364 @@ +//--------------------------------------------------------------- +// $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_h +#define __deal2__fe_dgp_h + +#include +#include +#include + +template class PolynomialSpace; +template class MappingQ; + + +/** + * Discontinuous tensor product elements based on equidistant support points. + * + * This is a discontinuous finite element using interpolating tensor + * product polynomials. The shape functions are Lagrangian + * interpolants of an equidistant grid of points on the unit cell. The + * points are numbered in lexicographical order, @p{x} running fastest. + * + * @author Guido Kanschat, Ralf Hartmann, 2001 + */ +template +class FE_DGP : public FiniteElement +{ + public: + /** + * Constructor for tensor product + * polynomials of degree @p{k}. + */ + FE_DGP (unsigned int k); + /** + * Destructor. + */ + ~FE_DGP (); + + /** + * Return the value of the + * @p{i}th shape function at the + * point @p{p}. @p{p} is a point + * on the reference element. + */ + virtual double shape_value (const unsigned int i, + const Point &p) const; + + /** + * Return the gradient of the + * @p{i}th shape function at the + * point @p{p}. @p{p} is a point + * on the reference element, and + * likewise the gradient is the + * gradient on the unit cell with + * respect to unit cell + * coordinates. + */ + virtual Tensor<1,dim> shape_grad (const unsigned int i, + const Point &p) const; + + /** + * Return the tensor of second + * derivatives of the @p{i}th + * shape function at point @p{p} + * on the unit cell. The + * derivatives are derivatives on + * the unit cell with respect to + * unit cell coordinates. + */ + virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, + const Point &p) 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; + + 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: + + /** + * Declare a nested class which + * will has static definitions of + * various matrices such as + * constraint and embedding + * matrices. The definition of + * the various static fields are + * in the files @p{fe_q_[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[]; + + /** + * 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[]; + + /** + * As + * @p{n_embedding_matrices} + * but for projection + * matrices. + */ + static const unsigned int n_projection_matrices; + }; + + + /** + * 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); + + /** + * Compute flags for initial update only. + */ + virtual UpdateFlags update_once (const UpdateFlags flags) const; + + /** + * Compute flags for update on each cell. + */ + virtual UpdateFlags update_each (const UpdateFlags flags) const; + + /** + * Degree of the polynomials. + */ + const unsigned int degree; + + /** + * Pointer to the tensor + * product polynomials. + */ + PolynomialSpace* poly; + + /** + * Fields of cell-independent data. + */ + class InternalData : public FiniteElementBase::InternalDataBase + { + public: + /** + * Array with shape function + * values in quadrature + * points. There is one + * vector for each shape + * function, containing + * values for each quadrature + * point. + */ + std::vector > shape_values; + + /** + * Array with shape function + * gradients in quadrature + * points. There is one + * vector for each shape + * function, containing + * values for each quadrature + * point. + */ + typename std::vector > > shape_gradients; + }; + + /** + * Allow access from other dimensions. + */ + template friend class FE_DGP; + + /** + * Allows @p{MappingQ} class to + * access to build_renumbering + * function. + */ + friend class MappingQ; +}; + + +// declaration of explicit specializations + +template <> +const double * const FE_DGP<1>::Matrices::embedding[]; + +template <> +const unsigned int FE_DGP<1>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_DGP<1>::Matrices::projection_matrices[]; + +template <> +const unsigned int FE_DGP<1>::Matrices::n_projection_matrices; + +template <> +const double * const FE_DGP<2>::Matrices::embedding[]; + +template <> +const unsigned int FE_DGP<2>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_DGP<2>::Matrices::projection_matrices[]; + +template <> +const unsigned int FE_DGP<2>::Matrices::n_projection_matrices; + +template <> +const double * const FE_DGP<3>::Matrices::embedding[]; + +template <> +const unsigned int FE_DGP<3>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_DGP<3>::Matrices::projection_matrices[]; + +template <> +const unsigned int FE_DGP<3>::Matrices::n_projection_matrices; + + +#endif diff --git a/deal.II/deal.II/source/fe/fe_dgp.cc b/deal.II/deal.II/source/fe/fe_dgp.cc new file mode 100644 index 0000000000..80372c2c5c --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_dgp.cc @@ -0,0 +1,411 @@ +//---------------------------------------------------------------- +// $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 +#include +#include +#include + + + +template +FE_DGP::FE_DGP (unsigned int degree) + : + FiniteElement (FiniteElementData(get_dpo_vector(degree),1), + std::vector(1,true)), + degree(degree), + poly(0) +{ + // create array of Legendre polynomials + std::vector > v; + for (unsigned int i=0;i<=degree;++i) + v.push_back(Legendre(i)); + poly = new PolynomialSpace (v); + + // if defined, copy over matrices + // from precomputed arrays +// if ((degree < Matrices::n_embedding_matrices) && +// (Matrices::embedding[degree] != 0)) +// { +// prolongation[0].fill (Matrices::embedding[degree]); +// } +// else +// // matrix undefined, set size to zero +// for (unsigned int i=0;i::children_per_cell;++i) +// prolongation[i].reinit(0); + +// // same as above: copy over matrix +// // from predefined values and +// // generate all others by rotation +// 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); +}; + + + +template +FE_DGP::~FE_DGP () +{ + // delete poly member and set it to + // zero to prevent accidental use + delete poly; + poly = 0; +} + + + +template +FiniteElement * +FE_DGP::clone() const +{ + return new FE_DGP(degree); +} + + + +template +double +FE_DGP::shape_value (const unsigned int i, + const Point &p) const +{ + return poly->compute_value(i, p); +} + + + +template +Tensor<1,dim> +FE_DGP::shape_grad (const unsigned int i, + const Point &p) const +{ + return poly->compute_grad(i, p); +} + + + +template +Tensor<2,dim> +FE_DGP::shape_grad_grad (const unsigned int i, + const Point &p) const +{ + return poly->compute_grad_grad(i, p); +} + + +//---------------------------------------------------------------------- +// Auxiliary functions +//---------------------------------------------------------------------- + + +template +std::vector +FE_DGP::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_DGP::update_once (const UpdateFlags flags) const +{ + UpdateFlags out = update_default; + + if (flags & update_values) + out |= update_values; + + return out; +} + + +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 +{ + InternalData* data = new InternalData; + std::vector values(0); + std::vector > grads(0); + std::vector > grad_grads(0); + + // 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; + + if (flags & update_values) + { + values.resize (dofs_per_cell); + data->shape_values.resize(dofs_per_cell, + std::vector(n_q_points)); + } + + if (flags & update_gradients) + { + grads.resize (dofs_per_cell); + data->shape_gradients.resize(dofs_per_cell, + std::vector >(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); + + + if (flags & (update_values | update_gradients)) + for (unsigned int i=0; icompute(quadrature.point(i), values, grads, grad_grads); + for (unsigned int k=0; kshape_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; k +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 unsigned int offset = face * quadrature.n_quadrature_points; + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + for (unsigned int k=0; k +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 unsigned int offset = (face * GeometryInfo::subfaces_per_face + subface) + * quadrature.n_quadrature_points; + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + for (unsigned int k=0; k +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 +bool +FE_DGP::has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const +{ + return true; +} + + + +template +unsigned int +FE_DGP::memory_consumption () const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + + +template +unsigned int +FE_DGP::get_degree () const +{ + return degree; +}; + + + + +template class FE_DGP; -- 2.39.5