From 9e9c63f8ae1765208158ea7db7154d9c20f908cb Mon Sep 17 00:00:00 2001 From: guido Date: Tue, 24 May 2005 10:20:27 +0000 Subject: [PATCH] new RT compiles and links git-svn-id: https://svn.dealii.org/trunk@10719 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_poly_tensor.h | 64 +++- .../deal.II/include/fe/fe_raviart_thomas.h | 142 +++++++-- deal.II/deal.II/source/fe/fe_poly_tensor.cc | 90 ++++-- .../source/fe/fe_raviart_thomas_nodal.cc | 294 ++++++++++++++++++ 4 files changed, 525 insertions(+), 65 deletions(-) create mode 100644 deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.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 index f18a0061f5..4c8559023b 100644 --- a/deal.II/deal.II/include/fe/fe_poly_tensor.h +++ b/deal.II/deal.II/include/fe/fe_poly_tensor.h @@ -63,6 +63,11 @@ class FE_PolyTensor : public FiniteElement public: /** * Constructor. + * + * @arg @c degree: constructor + * argument for poly. May be + * different from @p + * fe_data.degree. */ FE_PolyTensor (unsigned int degree, const FiniteElementData &fe_data, @@ -134,6 +139,28 @@ class FE_PolyTensor : public FiniteElement */ virtual unsigned int element_multiplicity (const unsigned int index) const; + /** + * Given flags, + * determines the values which + * must be computed only for the + * reference cell. Make sure, + * that #mapping_type is set by + * the derived class, such that + * this function can operate + * correctly. + */ + virtual UpdateFlags update_once (const UpdateFlags flags) const; + /** + * Given flags, + * determines the values which + * must be computed in each cell + * cell. Make sure, that + * #mapping_type is set by the + * derived class, such that this + * function can operate + * correctly. + */ + virtual UpdateFlags update_each (const UpdateFlags flags) const; protected: /** @@ -150,11 +177,36 @@ class FE_PolyTensor : public FiniteElement * possible to avoid the mapping. */ enum MappingType { - /// Shape functions do not depend on actual mesh cell + /** + * No mapping has been + * selected, throw an error + * if needed. + */ + no_mapping, + /** + * Shape functions do not + * depend on actual mesh + * cell + */ independent, - /// Shape functions are transformed covariant. + /** + * Shape functions do not + * depend on actual mesh + * cell. The mapping class + * must be + * MappingCartesian. + */ + independent_on_cartesian, + /** + * Shape functions are + * transformed covariant. + */ covariant, - /// Shape functions are transformed contravariant. + /** + * Shape functions are + * transformed + * contravariant. + */ contravariant }; @@ -240,17 +292,13 @@ class FE_PolyTensor : public FiniteElement 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 diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h index d9b52ead79..3a8cf09a52 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -14,9 +14,11 @@ #define __deal2__fe_raviart_thomas_h #include +#include #include #include #include +#include template class MappingQ; @@ -258,16 +260,12 @@ class FE_RaviartThomas : public FiniteElement /** * Check whether a shape function - * is non-zero on a face. + * may be non-zero on a face. * * Right now, this is only * implemented for RT0 in * 1D. Otherwise, returns always * @p true. - * - * Implementation of the - * interface in - * FiniteElement */ virtual bool has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const; @@ -291,31 +289,15 @@ class FE_RaviartThomas : public FiniteElement DeclException0 (ExcNotUsefulInThisDimension); 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 - * FiniteElement. - */ virtual void fill_fe_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, @@ -324,11 +306,6 @@ class FE_RaviartThomas : public FiniteElement typename Mapping::InternalDataBase &fe_internal, FEValuesData& data) const; - /** - * Implementation of the same - * function in - * FiniteElement. - */ virtual void fill_fe_face_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, @@ -338,11 +315,6 @@ class FE_RaviartThomas : public FiniteElement typename Mapping::InternalDataBase &fe_internal, FEValuesData& data) const; - /** - * Implementation of the same - * function in - * FiniteElement. - */ virtual void fill_fe_subface_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, @@ -570,6 +542,112 @@ class FE_RaviartThomas : public FiniteElement template friend class FE_RaviartThomas; }; + + +/** + * The Raviart-Thomas elements with node functionals defined as point + * values in Gauss points. + * + *

Description of node values

+ * + * For this Raviart-Thomas element, the node values are not cell and + * face moments with respect to certain polynomials, but the values in + * quadrature points. + * + * For an RT-element of degree k, we choose + * k+1d-1 Gauss points on each face. This way, the + * normal component which is in Qk is uniquely + * determined. Furthermore, since this Gauss-formula is exact on + * Q2k+1, these node values correspond to the exact + * integration of the moments of the RT-space. + * + * In the interior of the cells, the moments are with respect to an + * anisotropic Qk space, where the test functions + * are one degree lower in the direction corresponding to the vector + * component under consideration. This can be emulated by using an + * anisotropic Gauss formula for integration. + * + * @warning The degree stored in the member variable + * FiniteElementData::degree is higher by one than the + * constructor argument! + * + * @author Guido Kanschat, 2005 + */ +template +class FE_RaviartThomasNodal + : + public FE_PolyTensor, dim> +{ + /** + * Constructor for the Raviart-Thomas + * element of degree @p p. + */ + FE_RaviartThomasNodal (const unsigned int p); + + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * FE_RaviartThomasNodal(degree), with + * @p dim and @p degree + * replaced by appropriate + * values. + */ + virtual std::string get_name () const; + + virtual FiniteElement* clone () const; + + /** + * Check whether a shape function + * may be non-zero on a face. + * + * Right now, always returns + * @p true. + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) 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 (const unsigned int degree); + + /** + * Compute the vector used for + * the + * @p restriction_is_additive + * field passed to the base + * class's constructor. + */ + static std::vector + get_ria_vector (const unsigned int degree); + /** + * Initialize the + * FiniteElementBase::unit_support_points + * and FiniteElementBase::unit_face_support_points + * fields. Called from the + * constructor. + */ + void initialize_unit_support_points (const unsigned int degree); + + /** + * Initialize the + * #inverse_node_matrix + * field. Called from the + * constructor. + */ + void initialize_node_matrix (); +}; + + /*@}*/ /* -------------- declaration of explicit specializations ------------- */ diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc index 4c9cebdd16..b6fb9f3a8f 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -26,7 +26,6 @@ FE_PolyTensor::FE_PolyTensor ( FiniteElement (fe_data, restriction_is_additive_flags, nonzero_components), - degree(degree), poly_space(POLY(degree)) { cached_point(0) = -1; @@ -252,12 +251,13 @@ FE_PolyTensor::get_data (const UpdateFlags update_flags, 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 +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 @@ -269,7 +269,8 @@ FE_PolyTensor::fill_fe_values (const Mapping &m const unsigned int n_quad = quadrature.n_quadrature_points; const UpdateFlags flags(fe_data.current_update_flags()); - Assert(mapping_type == independent, ExcNotImplemented()); + Assert(mapping_type == independent || mapping_type == independent_on_cartesian, + ExcNotImplemented()); Assert(!(flags & update_values) || fe_data.shape_values.size() == n_dofs, ExcDimensionMismatch(fe_data.shape_values.size(), n_dofs)); @@ -307,13 +308,14 @@ FE_PolyTensor::fill_fe_values (const Mapping &m 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 +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 @@ -334,7 +336,8 @@ FE_PolyTensor::fill_fe_face_values (const Mapping const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - Assert(mapping_type == independent, ExcNotImplemented()); + Assert(mapping_type == independent || mapping_type == independent_on_cartesian, + ExcNotImplemented()); //TODO: Size assertions for (unsigned int i=0; i::fill_fe_face_values (const Mapping 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 +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 @@ -393,7 +397,8 @@ FE_PolyTensor::fill_fe_subface_values (const Mapping const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - Assert(mapping_type == independent, ExcNotImplemented()); + Assert(mapping_type == independent || mapping_type == independent_on_cartesian, + ExcNotImplemented()); //TODO: Size assertions for (unsigned int i=0; i::element_multiplicity (const unsigned int index) const } +template +UpdateFlags +FE_PolyTensor::update_once (const UpdateFlags flags) const +{ + Assert (mapping_type != no_mapping, ExcNotInitialized()); + const bool values_once = (mapping_type == independent); + + UpdateFlags out = update_default; + + if (values_once && (flags & update_values)) + out |= update_values; + + return out; +} + + +template +UpdateFlags +FE_PolyTensor::update_each (const UpdateFlags flags) const +{ + Assert (mapping_type != no_mapping, ExcNotInitialized()); + const bool values_once = (mapping_type == independent); + + UpdateFlags out = update_default; + + if (!values_once && (flags & update_values)) + out |= update_values | update_covariant_transformation; + if (flags & update_gradients) + out |= update_gradients | update_covariant_transformation; + if (flags & update_second_derivatives) + out |= update_second_derivatives | update_covariant_transformation; + + return out; +} + template class FE_PolyTensor,deal_II_dimension>; diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc new file mode 100644 index 0000000000..1570bf7df0 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc @@ -0,0 +1,294 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003, 2004, 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 +#include +#include +#include +#include +#include +#include + +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + + +template +FE_RaviartThomasNodal::FE_RaviartThomasNodal (const unsigned int deg) + : + FE_PolyTensor, dim> ( + deg, + FiniteElementData(get_dpo_vector(deg), + dim, deg+1), + get_ria_vector (deg), + std::vector >( + FiniteElementData(get_dpo_vector(deg), + dim,deg+1).dofs_per_cell, + std::vector(dim,true))) +{ + Assert (dim >= 2, ExcImpossibleInDim(dim)); + + this->mapping_type = independent_on_cartesian; + + for (unsigned int i=0; i::children_per_cell; ++i) + this->prolongation[i].reinit (this->dofs_per_cell, + this->dofs_per_cell); + FETools::compute_embedding_matrices (*this, &this->prolongation[0]); + // finally fill in support points + // on cell and face + initialize_unit_support_points (deg); + initialize_node_matrix(); +} + + + +template +std::string +FE_RaviartThomasNodal::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 + +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FE_RaviartThomasNodal<" << dim << ">(" << degree-1 << ")"; + +#ifndef HAVE_STD_STRINGSTREAM + namebuf << std::ends; +#endif + return namebuf.str(); +} + + + +template +bool +FE_RaviartThomasNodal::has_support_on_face (unsigned int, unsigned int) const +{ + return true; +} + + + +template +FiniteElement * +FE_RaviartThomasNodal::clone() const +{ + return new FE_RaviartThomasNodal(degree-1); +} + + +#if deal_II_dimension == 1 + +template <> +std::vector +FE_RaviartThomasNodal<1>::get_dpo_vector (const unsigned int deg) +{ + std::vector dpo(2); + dpo[0] = 1; + dpo[1] = deg; + return dpo; +} + +#endif + + +template +std::vector +FE_RaviartThomasNodal::get_dpo_vector (const unsigned int deg) +{ + // the element is face-based and we have + // (deg+1)^(dim-1) DoFs per face + unsigned int dofs_per_face = 1; + for (unsigned int d=0; d dpo(dim+1); + dpo[dim-1] = dofs_per_face; + dpo[dim] = interior_dofs; + + return dpo; +} + + + +#if deal_II_dimension == 1 + +template <> +std::vector +FE_RaviartThomasNodal<1>::get_ria_vector (const unsigned int) +{ + Assert (false, ExcImpossibleInDim(1)); + return std::vector(); +} + +#endif + + +template +std::vector +FE_RaviartThomasNodal::get_ria_vector (const unsigned int deg) +{ + unsigned int dofs_per_cell, dofs_per_face; + switch (dim) + { + case 2: + dofs_per_face = deg+1; + dofs_per_cell = 2*(deg+1)*(deg+2); + break; + case 3: + dofs_per_face = (deg+1)*(deg+1); + dofs_per_cell = 3*(deg+1)*(deg+1)*(deg+2); + break; + default: + Assert (false, ExcNotImplemented()); + } + Assert (FiniteElementData(get_dpo_vector(deg),dim).dofs_per_cell == + dofs_per_cell, + ExcInternalError()); + Assert (FiniteElementData(get_dpo_vector(deg),dim).dofs_per_face == + dofs_per_face, + ExcInternalError()); + + // all face dofs need to be + // non-additive, since they have + // continuity requirements. + // however, the interior dofs are + // made additive + std::vector ret_val(dofs_per_cell,false); + for (unsigned int i=GeometryInfo::faces_per_cell*dofs_per_face; + i < dofs_per_cell; ++i) + ret_val[i] = true; + + return ret_val; +} + + +template +void +FE_RaviartThomasNodal::initialize_unit_support_points (const unsigned int deg) +{ + this->unit_support_points.resize (this->dofs_per_cell); + this->unit_face_support_points.resize (this->dofs_per_face); + + unsigned int current = 0; + + // On the faces, we choose as many + // Gauss points as necessary to + // determine the normal component + // uniquely. This is the deg of + // the Raviart-Thomas element plus + // one. + if (dim>1) + { + QGauss face_points (deg+1); + Assert (face_points.n_quadrature_points == this->dofs_per_face, + ExcInternalError()); + for (unsigned int k=0;kdofs_per_face;++k) + this->unit_face_support_points[k] = face_points.point(k); + Quadrature faces = QProjector::project_to_all_faces(face_points); + for (unsigned int k=0;kunit_support_points[k] = faces.point(k); + + current = faces.n_quadrature_points; + } + // In the interior, we need + // anisotropic Gauss quadratures, + // different for each direction. + QGauss<1> high(deg+1); + QGauss<1> low(deg); + + for (unsigned int d=0;d* quadrature; + if (dim == 1) quadrature = new QAnisotropic(high); + if (dim == 2) quadrature = new QAnisotropic(((d==0) ? low : high), + ((d==1) ? low : high)); + if (dim == 3) quadrature = new QAnisotropic(((d==0) ? low : high), + ((d==1) ? low : high), + ((d==2) ? low : high)); + for (unsigned int k=0;kn_quadrature_points;++k) + this->unit_support_points[current++] = quadrature->point(k); + delete quadrature; + } + Assert (current == this->dofs_per_cell, ExcInternalError()); +} + + +template +void +FE_RaviartThomasNodal::initialize_node_matrix () +{ + const unsigned int n_dofs = this->dofs_per_cell; + + // We use an auxiliary matrix in + // this function. Therefore, + // inverse_node_matrix is still + // empty and shape_value_component + // returns the 'raw' shape values. + FullMatrix N(n_dofs, n_dofs); + + // The curent node functional index + unsigned int current = 0; + + // For each face and all quadrature + // points on it, the node value is + // the normal component of the + // shape function, possibly + // pointing in negative direction. + for (unsigned int face = 0; face::faces_per_cell; ++face) + for (unsigned int k=0;kdofs_per_face;++k) + { + for (unsigned int i=0;ishape_value_component( + i, unit_support_point(current), + GeometryInfo< dim >::unit_normal_direction[face]) + * GeometryInfo< dim >::unit_normal_orientation[face]; + ++current; + } + // Interior degrees of freedom in each direction + const unsigned int n_cell = (n_dofs - current) / dim; + + for (unsigned int d=0;dshape_value_component(i, unit_support_point(current), d); + ++current; + } + Assert (current == n_dofs, ExcInternalError()); + + inverse_node_matrix.invert(N); +} + + + +template FE_RaviartThomasNodal; -- 2.39.5