From 2d88e8704183a59f1ed00561e31e8ec5bbea2a23 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Sat, 24 Sep 2005 11:26:56 +0000 Subject: [PATCH] new interpolation functions; documentation git-svn-id: https://svn.dealii.org/trunk@11526 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 41 +++++ .../deal.II/include/fe/fe_raviart_thomas.h | 17 ++ .../source/fe/fe_raviart_thomas_nodal.cc | 152 +++++++++++++----- 3 files changed, 172 insertions(+), 38 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 9c0115d5ed..19040d23ed 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -90,6 +90,10 @@ template class FECollection; * *

Notes on the implementation of derived classes

* + * The following sections list the information to be provided by + * derived classes, depending on the space dimension. They are + * followed by a list of functions helping to generate these values. + * *

Finite elements in one dimension

* * Finite elements in one dimension need only set the #restriction @@ -213,6 +217,43 @@ template class FECollection; * constraints that are entered more than once (as is necessary for the case * above), it insists that the weights are exactly the same. * + *

Helper functions

+ * + * Construction of a finite element and computation of the matrices + * described above may be a tedious task, in particular if it has to + * be performed for several space dimensions. Therefore, some + * functions in FETools have been provided to help with these tasks. + * + * First, aready the basis of the shape function space may be + * difficult to implement for arbitrary order and dimension. On the + * other hand, if the @ref GlossNodes "node values" are given, then + * the duality relation between node functionals and basis functions + * defines the basis. As a result, the shape function space may be + * defined with arbitrary "raw" basis functions, such that the actual + * finite element basis is computed from linear combinations of + * them. The coefficients of these combinations are determined by the + * duality of node values. + * + * Using this matrix allows the construction of the basis of shape + * functions in two steps. + *
    + * + *
  1. Define the space of shape functions using an arbitrary basis + * wj and compute the matrix M of node + * functionals Ni applied to these basis functions, + * such that its entries are mij = + * Ni(wj). + * + *
  2. Compute the basis vj of the finite element + * shape function space by applying M-1 to the basis + * wj. + *
+ * + * The function computing the matrix M for you is + * FETools::compute_node_matrix(). It relies on the existence of + * #generalized_support_points and implementation of interpolate() + * with VectorSlice argument. + * * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1998, 2000, 2001, 2005 */ template 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 d01d299ac0..bfb805ff94 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -334,6 +334,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 VectorSlice > >& values) const; private: /** * The order of the @@ -630,6 +639,14 @@ class FE_RaviartThomasNodal */ 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; private: /** * Only for internal use. Its 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 index e22c31beb4..08fd82fd0d 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc @@ -52,7 +52,7 @@ FE_RaviartThomasNodal::FE_RaviartThomasNodal (const unsigned int deg) // basis functions initialize_unit_support_points (deg); initialize_node_matrix(); - + for (unsigned int i=0; i::children_per_cell; ++i) this->prolongation[i].reinit (this->dofs_per_cell, this->dofs_per_cell); @@ -96,6 +96,111 @@ FE_RaviartThomasNodal::has_support_on_face (unsigned int, unsigned int) con } +template +void +FE_RaviartThomasNodal::interpolate( + std::vector&, + const std::vector&) const +{ + Assert(false, ExcNotImplemented()); +} + + +template +void +FE_RaviartThomasNodal::interpolate( + std::vector& local_dofs, + const std::vector >& values, + unsigned int offset) const +{ + Assert (values.size() == generalized_support_points.size(), + ExcDimensionMismatch(values.size(), 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())); + + // First do interpolation on + // faces. There, the component + // evaluated depends on the face + // direction and orientation. + unsigned int fbase = 0; + unsigned int f=0; + for (;f::faces_per_cell; + ++f, fbase+=this->dofs_per_face) + { + for (unsigned int i=0;idofs_per_face;++i) + { + local_dofs[fbase+i] = values[fbase+i](offset+GeometryInfo::unit_normal_direction[f]); + } + } + + // The remaining points form dim + // chunks, one for each component. + const unsigned int istep = (this->dofs_per_cell - fbase) / dim; + Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError()); + + f = 0; + while (fbase < this->dofs_per_cell) + { + for (unsigned int i=0;idofs_per_cell, ExcInternalError()); +} + + + +template +void +FE_RaviartThomasNodal::interpolate( + std::vector& local_dofs, + const VectorSlice > >& values) const +{ + Assert (values.size() == this->n_components(), + ExcDimensionMismatch(values.size(), this->n_components())); + Assert (values[0].size() == generalized_support_points.size(), + ExcDimensionMismatch(values.size(), generalized_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + // First do interpolation on + // faces. There, the component + // evaluated depends on the face + // direction and orientation. + unsigned int fbase = 0; + unsigned int f=0; + for (;f::faces_per_cell; + ++f, fbase+=this->dofs_per_face) + { + for (unsigned int i=0;idofs_per_face;++i) + { + local_dofs[fbase+i] = values[GeometryInfo::unit_normal_direction[f]][fbase+i]; + } + } + // The remaining points form dim + // chunks, one for each component. + const unsigned int istep = (this->dofs_per_cell - fbase) / dim; + Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError()); + + f = 0; + while (fbase < this->dofs_per_cell) + { + for (unsigned int i=0;idofs_per_cell, ExcInternalError()); +} + + + template FiniteElement * @@ -199,8 +304,8 @@ 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); + this->generalized_support_points.resize (this->dofs_per_cell); + this->generalized_face_support_points.resize (this->dofs_per_face); unsigned int current = 0; @@ -216,11 +321,11 @@ FE_RaviartThomasNodal::initialize_unit_support_points (const unsigned int d 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); + 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->unit_support_points[k] = faces.point(k); + this->generalized_support_points[k] = faces.point(k); current = this->dofs_per_face*GeometryInfo::faces_per_cell; } @@ -242,7 +347,7 @@ FE_RaviartThomasNodal::initialize_unit_support_points (const unsigned int d ((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); + this->generalized_support_points[current++] = quadrature->point(k); delete quadrature; } Assert (current == this->dofs_per_cell, ExcInternalError()); @@ -254,44 +359,15 @@ 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, this->unit_support_points[current], - GeometryInfo< dim >::unit_normal_direction[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, this->unit_support_points[current], d); - ++current; - } - Assert (current == n_dofs, ExcInternalError()); + FullMatrix M(n_dofs, n_dofs); + FETools::compute_node_matrix(M, *this); this->inverse_node_matrix.reinit(n_dofs, n_dofs); - this->inverse_node_matrix.invert(N); + this->inverse_node_matrix.invert(M); } -- 2.39.5