From: Guido Kanschat Date: Mon, 17 Oct 2005 15:02:30 +0000 (+0000) Subject: working on interpolation X-Git-Tag: v8.0.0~12990 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8b9c34da5d442369f9d3cac9ae8b76078576abf0;p=dealii.git working on interpolation git-svn-id: https://svn.dealii.org/trunk@11613 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 6b50cae56b..31224d5ed6 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -15,11 +15,14 @@ #include #include +#include #include #include #include #include +#include + template class MappingQ; @@ -325,15 +328,15 @@ 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 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; + virtual void interpolate( + std::vector& local_dofs, + const VectorSlice > >& values) const; private: /** * The order of the @@ -548,6 +551,23 @@ class FE_RaviartThomas : public FiniteElement */ std::vector > > shape_gradients; }; + + /** + * 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. + */ + QGauss face_quadrature; + + /** + * Legendre polynomials are used + * for computing the moments. + */ + std::vector > legendre_polynomials; /** * 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 28aeebe84c..d0a3dc4bbf 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -12,6 +12,7 @@ //--------------------------------------------------------------------------- #include +#include #include #include #include @@ -138,7 +139,9 @@ 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)) + renumber (compute_renumber(rt_order)), + face_quadrature((dim>1) ? rt_order+1 : 0), + legendre_polynomials(rt_order) { Assert (dim >= 2, ExcImpossibleInDim(dim)); @@ -809,20 +812,54 @@ FE_RaviartThomas<3>::initialize_restriction () #endif -//TODO: Unit support points should go away according to new definition. template void -FE_RaviartThomas::initialize_support_points (const unsigned int) +FE_RaviartThomas::initialize_support_points (const unsigned int deg) { -// unsigned int n_face_points = (dim>1) ? 1 : 0; -// // compute (deg+1)^(dim-1) -// for (unsigned int d=1;d cell_quadrature(deg+1); + const unsigned int n_interior_points + = (deg>1) ? cell_quadrature.n_quadrature_points : 0; + + unsigned int n_face_points = (dim>1) ? 1 : 0; + // compute (deg+1)^(dim-1) + for (unsigned int d=1;dgeneralized_support_points.resize (GeometryInfo::faces_per_cell*n_face_points + +cell_quadrature.n_quadrature_points); + this->generalized_face_support_points.resize (n_face_points); + + // Number of the point being entered + unsigned int current = 0; + + 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->generalized_face_support_points[k] = face_points.point(k); + Quadrature 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; + } -// this->unit_generalized_support_points.resize (this->dofs_per_cell); -// this->generalized_face_support_points.resize (this->dofs_per_face); + 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()); + +//TODO: Unit support points should go away according to new definition. this->unit_support_points.resize (this->dofs_per_cell); switch (dim) @@ -1869,6 +1906,43 @@ FE_RaviartThomas::has_support_on_face (const unsigned int shape_index, +template +void +FE_RaviartThomas::interpolate( + std::vector&, + const std::vector&) const +{ + Assert(false, ExcNotImplemented()); +} + + + +template +void +FE_RaviartThomas::interpolate( + std::vector& local_dofs, + const std::vector >& values, + unsigned int offset) const +{ + Assert (values.size() == this->generalized_support_points.size(), + ExcDimensionMismatch(values.size(), this->generalized_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + Assert (values[0].size() >= offset+this->n_components(), + ExcDimensionMismatch(values[0].size(),offset+this->n_components())); + + std::fill(local_dofs.begin(), local_dofs.end(), 0.); +} + + +template +void +FE_RaviartThomas::interpolate( + std::vector& /*local_dofs*/, + const VectorSlice > >& /*values*/) const +{} + + template unsigned int FE_RaviartThomas::memory_consumption () const