From 675bc74fc1baf770f9b44d8f09236e9c82c1d3b8 Mon Sep 17 00:00:00 2001 From: guido Date: Tue, 24 May 2005 06:51:17 +0000 Subject: [PATCH] checking in preliminary version git-svn-id: https://svn.dealii.org/trunk@10717 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_poly_tensor.h | 312 ++++++++++++++ deal.II/deal.II/source/fe/fe_poly_tensor.cc | 454 ++++++++++++++++++++ 2 files changed, 766 insertions(+) create mode 100644 deal.II/deal.II/include/fe/fe_poly_tensor.h create mode 100644 deal.II/deal.II/source/fe/fe_poly_tensor.cc diff --git a/deal.II/deal.II/include/fe/fe_poly_tensor.h b/deal.II/deal.II/include/fe/fe_poly_tensor.h new file mode 100644 index 0000000000..f18a0061f5 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_poly_tensor.h @@ -0,0 +1,312 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005 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_tensor_h +#define __deal2__fe_poly_tensor_h + + +#include +#include + +/*!@addtogroup fe */ +/*@{*/ + +/** + * This class gives a unified framework for the implementation of + * FiniteElement classes based on Tensor valued polynomial spaces like + * PolynomialsBDM and PolynomialsRaviartThomas. + * + * Every class that implements following function can be used as + * template parameter POLY. + * + * @code + * void compute (const Point &unit_point, + * std::vector > &values, + * std::vector > &grads, + * std::vector > &grad_grads) const; + * @endcode + * + * The polynomial spaces are usually described as direct sum of + * simpler spaces. In such a case, the usual basis of node functionals + * is not dual to the basis of the polynomial space. Therefore, the + * matrix #inverse_node_matrix can be filled by the constructor of a + * derived class such that the usual interpolation condition + * Ni(vj) holds on the reference cell. + * + * In many cases, the node functionals depend on the shape of the mesh + * cell, since they evaluate normal or tangential components on the + * faces. In order to allow for a set of transformations, the variable + * #mapping_type has been introduced. It should also be set in the + * constructor of a derived class. + * + * This class is not a fully implemented FiniteElement class, but + * implements some common features of vector valued elements based on + * vector valued polynomial classes. What's missing here in particular + * is information on the topological location of the node values. + * + * @see PolynomialsBDM, PolynomialsRaviartThomas + * + * @author Guido Kanschat, 2005 + **/ +template +class FE_PolyTensor : public FiniteElement +{ + public: + /** + * Constructor. + */ + FE_PolyTensor (unsigned int degree, + const FiniteElementData &fe_data, + const std::vector &restriction_is_additive_flags, + const std::vector > &nonzero_components); + + /** + * Since these elements are + * vector valued, an exception is + * thrown. + */ + virtual double shape_value (const unsigned int i, + const Point &p) const; + + virtual double shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Since these elements are + * vector valued, an exception is + * thrown. + */ + virtual Tensor<1,dim> shape_grad (const unsigned int i, + const Point &p) const; + + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Since these elements are + * vector valued, an exception is + * thrown. + */ + virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, + const Point &p) const; + + 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 not a composed element, + * return one. + */ + virtual unsigned int n_base_elements () const; + + /** + * Access to base element + * objects. Since this element is + * not composed of several elements, + * 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 + * not a composed 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: + /** + * Different options for + * transforming the basis + * functions from the reference + * cell to the actual mesh cell. + * + * Most vector valued elements + * either transform shape + * functions to keep node values + * on edges meaningful. Still, in + * special cases, it may be + * possible to avoid the mapping. + */ + enum MappingType { + /// Shape functions do not depend on actual mesh cell + independent, + /// Shape functions are transformed covariant. + covariant, + /// Shape functions are transformed contravariant. + contravariant + }; + + /** + * The mapping type to be used to + * map shape functions from the + * reference cell to the msh + * cell. + */ + MappingType mapping_type; + + virtual + typename Mapping::InternalDataBase * + get_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) 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 ; + + /** + * Fields of cell-independent + * data for FE_PolyTensor. Stores + * the values of the shape + * functions and their + * derivatives on the reference + * cell for later use. + * + * All tables are organized in a + * way, that the value for shape + * function i at + * quadrature point k is + * accessed by indices + * (i,k). + */ + 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. + */ + std::vector > > shape_values; + + /** + * Array with shape function + * gradients in quadrature + * points. There is one + * row for each shape + * function, containing + * values for each quadrature + * point. + */ + std::vector > > shape_grads; + }; + + /** + * Degree of the polynomials. + */ + unsigned int degree; + + /** + * The polynomial space. Its type + * is given by the template + * parameter POLY. + */ + POLY poly_space; + /** + * The inverse of the matrix + * aij of node + * values Ni + * applied to polynomial + * pj. This + * matrix is used to convert + * polynomials in the "raw" basis + * provided in #poly_space to the + * basis dual to the node + * functionals on the reference cell. + * + * This object is not filled by + * FE_PolyTensor, but is a chance + * for a derived class to allow + * for reorganization of the + * basis functions. If it is left + * empty, the basis in + * #poly_space is used. + */ + FullMatrix inverse_node_matrix; + + /** + * If a shape function is + * computed at a single point, we + * must compute all of them to + * apply #inverse_node_matrix. In + * order to avoid too much + * overhead, we cache the point + * and the function values for + * the next evaluation. + */ + mutable Point cached_point; + + /** + * Cached shape function values after + * call to + * shape_value_component(). + */ + mutable std::vector > cached_values; + + /** + * Cached shape function gradients after + * call to + * shape_grad_component(). + */ + mutable std::vector > cached_grads; + + /** + * Cached second derivatives of + * shape functions after call to + * shape_grad_grad_component(). + */ + mutable std::vector > cached_grad_grads; +}; + +/*@}*/ + +#endif diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc new file mode 100644 index 0000000000..4c9cebdd16 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -0,0 +1,454 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005 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 + + +template +FE_PolyTensor::FE_PolyTensor ( + unsigned int degree, + 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), + degree(degree), + poly_space(POLY(degree)) +{ + cached_point(0) = -1; +} + + +template +double +FE_PolyTensor::shape_value ( + const unsigned int, const Point &) const +{ + Assert(false, typename FiniteElementBase::ExcFENotPrimitive()); + return 0.; +} + + +template +double +FE_PolyTensor::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 < dim, ExcIndexRange (component, 0, dim)); + + if (cached_point != p || cached_values.size() == 0) + { + cached_point = p; + cached_values.resize(poly_space.n()); + poly_space.compute(p, cached_values, cached_grads, cached_grad_grads); + } + + double s = 0; + if (inverse_node_matrix.n_cols() == 0) + return cached_values[i][component]; + else + for (unsigned int j=0;j +Tensor<1,dim> +FE_PolyTensor::shape_grad ( + const unsigned int, const Point &) const +{ + Assert(false, typename FiniteElementBase::ExcFENotPrimitive()); + return Tensor<1,dim>(); +} + + + +template +Tensor<1,dim> +FE_PolyTensor::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 < dim, ExcIndexRange (component, 0, dim)); + + if (cached_point != p || cached_grads.size() == 0) + { + cached_point = p; + cached_grads.resize(poly_space.n()); + poly_space.compute(p, cached_values, cached_grads, cached_grad_grads); + } + + Tensor<1,dim> s; + if (inverse_node_matrix.n_cols() == 0) + return cached_grads[i][component]; + else + for (unsigned int j=0;j +Tensor<2,dim> +FE_PolyTensor::shape_grad_grad ( + const unsigned int, const Point &) const +{ + Assert(false, typename FiniteElementBase::ExcFENotPrimitive()); + return Tensor<2,dim>(); +} + + + +template +Tensor<2,dim> +FE_PolyTensor::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 < dim, ExcIndexRange (component, 0, dim)); + + if (cached_point != p || cached_grad_grads.size() == 0) + { + cached_point = p; + cached_grad_grads.resize(poly_space.n()); + poly_space.compute(p, cached_values, cached_grads, cached_grad_grads); + } + + Tensor<2,dim> s; + if (inverse_node_matrix.n_cols() == 0) + return cached_grad_grads[i][component]; + else + for (unsigned int j=0;j +typename Mapping::InternalDataBase * +FE_PolyTensor::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.resize (this->dofs_per_cell); + for (unsigned int i=0;idofs_per_cell;++i) + data->shape_values[i].resize (n_q_points); + } + + if (flags & update_gradients) + { + grads.resize (this->dofs_per_cell); + data->shape_grads.resize (this->dofs_per_cell); + for (unsigned int i=0;idofs_per_cell;++i) + data->shape_grads[i].resize (n_q_points); + } + + // if second derivatives through + // finite differencing is required, + // then initialize some objects for + // that + if (flags & update_second_derivatives) + { +// grad_grads.resize (this->dofs_per_cell); + data->initialize_2nd (this, mapping, quadrature); + } + + // Compute shape function values + // and derivatives on the reference + // cell. Make sure, that for the + // node values N_i holds + // N_i(v_j)=\delta_ij for all basis + // functions v_j + if (flags & (update_values | update_gradients)) + for (unsigned int k=0; kdofs_per_cell; ++i) + data->shape_values[i][k] = values[i]; + else + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + data->shape_values[i][k] = inverse_node_matrix(i,j) * values[j]; + + if (flags & update_gradients) + if (inverse_node_matrix.n_cols() == 0) + for (unsigned int i=0; idofs_per_cell; ++i) + data->shape_grads[i][k] = grads[i]; + else + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + data->shape_grads[i][k] = inverse_node_matrix(i,j) * grads[j]; + } + return data; +} + + + + +//--------------------------------------------------------------------------- +// Fill data of FEValues +//--------------------------------------------------------------------------- + +template +void +FE_PolyTensor::fill_fe_values (const Mapping &mapping, + const typename Triangulation::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 unsigned int n_dofs = this->dofs_per_cell; + const unsigned int n_quad = quadrature.n_quadrature_points; + const UpdateFlags flags(fe_data.current_update_flags()); + + Assert(mapping_type == independent, ExcNotImplemented()); + + Assert(!(flags & update_values) || fe_data.shape_values.size() == n_dofs, + ExcDimensionMismatch(fe_data.shape_values.size(), n_dofs)); + Assert(!(flags & update_values) || fe_data.shape_values[0].size() == n_quad, + ExcDimensionMismatch(fe_data.shape_values[0].size(), n_quad)); + + for (unsigned int i=0; i > shape_grads1 (n_quad); + mapping.transform_covariant(fe_data.shape_grads[i], 0, + shape_grads1, + mapping_data); + for (unsigned int k=0; k::DataSetDescriptor dsd; + if (flags & update_second_derivatives) + this->compute_2nd (mapping, cell, dsd.cell(), + mapping_data, fe_data, data); +} + + + +template +void +FE_PolyTensor::fill_fe_face_values (const Mapping &mapping, + const typename Triangulation::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 unsigned int n_dofs = this->dofs_per_cell; + const unsigned int n_quad = quadrature.n_quadrature_points; + // 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), + n_quad); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + Assert(mapping_type == independent, ExcNotImplemented()); +//TODO: Size assertions + + for (unsigned int i=0; i > shape_grads1 (n_quad); + mapping.transform_covariant(fe_data.shape_grads[i], offset, + shape_grads1, mapping_data); + for (unsigned int k=0; kcompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +template +void +FE_PolyTensor::fill_fe_subface_values (const Mapping &mapping, + const typename Triangulation::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 unsigned int n_dofs = this->dofs_per_cell; + const unsigned int n_quad = quadrature.n_quadrature_points; + // 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), + n_quad); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + Assert(mapping_type == independent, ExcNotImplemented()); +//TODO: Size assertions + + for (unsigned int i=0; i > shape_grads1 (n_quad); + mapping.transform_covariant(fe_data.shape_grads[i], offset, + shape_grads1, mapping_data); + for (unsigned int k=0; kcompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +template +unsigned int +FE_PolyTensor::n_base_elements () const +{ + return 1; +} + + + +template +const FiniteElement & +FE_PolyTensor::base_element (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return *this; +} + + + +template +unsigned int +FE_PolyTensor::element_multiplicity (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return 1; +} + + + +template class FE_PolyTensor,deal_II_dimension>; -- 2.39.5