From b9d23ae88d8be6bbc9d43a465328d87cacd6ec4c Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Wed, 19 Oct 2005 04:42:07 +0000 Subject: [PATCH] implement interpolation with some surprise git-svn-id: https://svn.dealii.org/trunk@11619 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/fe/fe_raviart_thomas.h | 70 +++--- .../deal.II/source/fe/fe_raviart_thomas.cc | 212 ++++++++++-------- 2 files changed, 162 insertions(+), 120 deletions(-) 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 31224d5ed6..4b08664d9f 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -14,6 +14,7 @@ #define __deal2__fe_raviart_thomas_h #include +#include #include #include #include @@ -34,6 +35,10 @@ template class MappingQ; * space Hdiv. These elements generate vector fields with * normel components continuous between mesh cells. * + * @todo Right now, the description of node values and interpolation + * does not match the actual state of this class. This will be + * improved after some discussion. + * * We follow the usual definition of the degree of RT elements, which * denotes the polynomial degree of the largest complete polynomial * subspace contained in the RT space. Then, approciamtion order of @@ -278,6 +283,15 @@ class FE_RaviartThomas : public FiniteElement virtual bool has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const; + virtual void interpolate(std::vector& local_dofs, + const std::vector& values) const; + virtual void interpolate(std::vector& local_dofs, + const std::vector >& values, + unsigned int offset = 0) const; + + virtual void interpolate( + std::vector& local_dofs, + const VectorSlice > >& values) const; /** * Determine an estimate for the * memory consumption (in bytes) @@ -328,15 +342,6 @@ class FE_RaviartThomas : public FiniteElement typename Mapping::InternalDataBase &fe_internal, FEValuesData& data) const; - virtual void interpolate(std::vector& local_dofs, - const std::vector& values) const; - virtual void interpolate(std::vector& local_dofs, - const std::vector >& values, - unsigned int offset = 0) const; - - virtual void interpolate( - std::vector& local_dofs, - const VectorSlice > >& values) const; private: /** * The order of the @@ -440,11 +445,14 @@ class FE_RaviartThomas : public FiniteElement void initialize_restriction (); /** - * Initialize the - * @p unit_support_points field - * of the FiniteElement - * class. Called from the - * constructor. + * Initialize the @p + * generalized_support_points + * field of the FiniteElement + * class and fill the tables with + * interpolation weights + * (#boundary_weights and + * #interior_weights). Called + * from the constructor. */ void initialize_support_points (const unsigned int rt_degree); @@ -553,21 +561,31 @@ class FE_RaviartThomas : public FiniteElement }; /** - * The quadrature formula used to - * generate support points on - * faces and computing the - * moments on faces. Its number - * of points is one order higher - * than the degree of the RT - * space. + * These are the factors + * multiplied to a function in + * the + * #generalized_face_support_points + * when computing the + * integration. They are + * organized such that there is + * one row for each generalized + * face support point and one + * column for each degree of + * freedom on the face. */ - QGauss face_quadrature; - + Table<2, double> boundary_weights; /** - * Legendre polynomials are used - * for computing the moments. + * Precomputed factors for + * interpolation of interior + * degrees of freedom. The + * rationale for this Table is + * the same as for + * #boundary_weights. Only, this + * table has a third coordinate + * for the space direction of the + * component evaluated. */ - std::vector > legendre_polynomials; + Table<3, double> interior_weights; /** * Allow access from other diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc index d0a3dc4bbf..f69681b45a 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -139,9 +139,7 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int rt_order) std::vector(dim,true))), rt_order(rt_order), polynomials (create_polynomials(rt_order)), - renumber (compute_renumber(rt_order)), - face_quadrature((dim>1) ? rt_order+1 : 0), - legendre_polynomials(rt_order) + renumber (compute_renumber(rt_order)) { Assert (dim >= 2, ExcImpossibleInDim(dim)); @@ -812,14 +810,47 @@ FE_RaviartThomas<3>::initialize_restriction () #endif +#if deal_II_dimension == 1 template void FE_RaviartThomas::initialize_support_points (const unsigned int deg) { + return; + + Assert (false, ExcNotImplemented()); + QGauss cell_quadrature(deg+1); const unsigned int n_interior_points - = (deg>1) ? cell_quadrature.n_quadrature_points : 0; + = (deg>0) ? cell_quadrature.n_quadrature_points : 0; + + this->generalized_support_points.resize (2 + n_interior_points); + + // Number of the point being entered + unsigned int current = 0; + + + if (deg==0) return; + + interior_weights.reinit(TableIndices<3>(2+n_interior_points, 0, dim)); + + for (unsigned int k=0;kgeneralized_support_points[current++] = cell_quadrature.point(k); + + Assert (current == this->generalized_support_points.size(), + ExcInternalError()); +} + +#else + +// Version for 2d and higher. See above for 1d version +template +void +FE_RaviartThomas::initialize_support_points (const unsigned int deg) +{ + QGauss cell_quadrature(deg+1); + const unsigned int n_interior_points + = (deg>0) ? cell_quadrature.n_quadrature_points : 0; unsigned int n_face_points = (dim>1) ? 1 : 0; // compute (deg+1)^(dim-1) @@ -828,7 +859,7 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) this->generalized_support_points.resize (GeometryInfo::faces_per_cell*n_face_points - +cell_quadrature.n_quadrature_points); + + n_interior_points); this->generalized_face_support_points.resize (n_face_points); // Number of the point being entered @@ -837,106 +868,72 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) if (dim>1) { QGauss face_points (deg+1); + TensorProductPolynomials legendre + = Polynomials::Legendre::generate_complete_basis(deg); + + boundary_weights.reinit(n_face_points, legendre.n()); + Assert (face_points.n_quadrature_points == this->dofs_per_face, ExcInternalError()); - for (unsigned int k=0;kdofs_per_face;++k) - this->generalized_face_support_points[k] = face_points.point(k); + for (unsigned int k=0;kgeneralized_face_support_points[k] = face_points.point(k); + // Compute its quadrature + // contribution for each + // moment. + for (unsigned int i=0;i faces = QProjector::project_to_all_faces(face_points); - for (unsigned int k=0; - kdofs_per_face*GeometryInfo::faces_per_cell; - ++k) - this->generalized_support_points[k] = faces.point(k); - - current = this->dofs_per_face*GeometryInfo::faces_per_cell; + for (;current::faces_per_cell*n_face_points; + ++current) + { + // Enter the support point + // into the vector + this->generalized_support_points[current] = faces.point(current); + } } if (deg==0) return; - for (unsigned int k=0;kgeneralized_support_points[current++] = cell_quadrature.point(k); - - Assert (current == this->generalized_support_points.size(), - ExcInternalError()); + // Create Legendre basis for the + // space D_xi Q_k + std::vector* > polynomials(dim); + for (unsigned int dd=0;dd > > poly(dim); + for (unsigned int d=0;d(poly); + } + + interior_weights.reinit(TableIndices<3>(n_interior_points, polynomials[0]->n(), dim)); - this->unit_support_points.resize (this->dofs_per_cell); - switch (dim) + for (unsigned int k=0;kdofs_per_face, ExcInternalError()); - - // associate support points - // with mid-face points if a - // shape function has a - // non-zero normal component - // there, otherwise with the - // cell center. the reason - // for this non-unique - // support point is that we - // use hierarchical shape - // functions, rather than - // Lagrange functions, for - // which we get into the same - // trouble as in the - // FE_Q_Hierarchical element; - // see the respective - // function there - - // start with the face shape - // functions - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[0*this->dofs_per_face+i] = Point(.5, .0); - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[1*this->dofs_per_face+i] = Point(1., .5); - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[2*this->dofs_per_face+i] = Point(.5, 1.); - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[3*this->dofs_per_face+i] = Point(.0, .5); - - // associate the rest with - // the cell center - for (unsigned int i=4*this->dofs_per_face; idofs_per_cell; ++i) - this->unit_support_points[i] = Point(.5, .5); - - break; - } - - case 3: - { - // same as in 2d - Assert ((rt_order+1)*(rt_order+1) == this->dofs_per_face, ExcInternalError()); - - // start with the face shape - // functions - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[0*this->dofs_per_face+i] = Point(.5, .0, .5); - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[1*this->dofs_per_face+i] = Point(.5, 1., .5); - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[2*this->dofs_per_face+i] = Point(.5, .5, 0.); - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[3*this->dofs_per_face+i] = Point(1., .5, .5); - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[4*this->dofs_per_face+i] = Point(.5, .5, 1.); - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_support_points[5*this->dofs_per_face+i] = Point(.0, .5, .5); - - // associate the rest with - // the cell center - for (unsigned int i=6*this->dofs_per_face; idofs_per_cell; ++i) - this->unit_support_points[i] = Point(.5, .5, .5); - - break; - } + this->generalized_support_points[current++] = cell_quadrature.point(k); + for (unsigned int i=0;in();++i) + for (unsigned int d=0;dcompute_value(i,cell_quadrature.point(k)); + } - default: - Assert (false, ExcNotImplemented()); - }; + for (unsigned int d=0;dgeneralized_support_points.size(), + ExcInternalError()); } +#endif #if deal_II_dimension == 1 @@ -1932,15 +1929,42 @@ FE_RaviartThomas::interpolate( ExcDimensionMismatch(values[0].size(),offset+this->n_components())); std::fill(local_dofs.begin(), local_dofs.end(), 0.); + + const unsigned int n_face_points = boundary_weights.size(0); + for (unsigned int face=0;face::faces_per_cell;++face) + for (unsigned int k=0;k::unit_normal_direction[face]); + } } template void FE_RaviartThomas::interpolate( - std::vector& /*local_dofs*/, - const VectorSlice > >& /*values*/) const -{} + std::vector& local_dofs, + const VectorSlice > >& values) const +{ + Assert (values.size() == this->n_components(), + ExcDimensionMismatch(values.size(), this->n_components())); + Assert (values[0].size() == this->generalized_support_points.size(), + ExcDimensionMismatch(values[0].size(), this->generalized_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + + std::fill(local_dofs.begin(), local_dofs.end(), 0.); + + const unsigned int n_face_points = boundary_weights.size(0); + for (unsigned int face=0;face::faces_per_cell;++face) + for (unsigned int k=0;k::unit_normal_direction[face]][face*n_face_points+k]; + } +} template -- 2.39.5