From 21e4f76f838a3e20efba6c5d87562c9c32e923c7 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 25 Oct 2005 08:19:16 +0000 Subject: [PATCH] implement Raviart-Thomas basis functions for standard RT interpolation git-svn-id: https://svn.dealii.org/trunk@11657 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/fe/fe_raviart_thomas.h | 285 +-- .../deal.II/source/fe/fe_raviart_thomas.cc | 1569 +---------------- .../source/fe/fe_raviart_thomas_nodal.cc | 19 +- 3 files changed, 79 insertions(+), 1794 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 4b08664d9f..70f0206611 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -136,7 +136,9 @@ template class MappingQ; * @author Wolfgang Bangerth, 2003 */ template -class FE_RaviartThomas : public FiniteElement +class FE_RaviartThomas + : + public FE_PolyTensor, dim> { public: /** @@ -156,92 +158,6 @@ class FE_RaviartThomas : public FiniteElement */ virtual std::string get_name () const; - /** - * Return the value of the - * @p componentth vector - * component of the @p ith shape - * function at the point - * @p p. See the - * FiniteElement base - * class for more information - * about the semantics of this - * function. - */ - virtual double shape_value_component (const unsigned int i, - const Point &p, - const unsigned int component) const; - - /** - * Return the gradient of the - * @p componentth vector - * component of the @p ith shape - * function at the point - * @p p. See the - * FiniteElement base - * class for more information - * about the semantics of this - * function. - */ - virtual Tensor<1,dim> shape_grad_component (const unsigned int i, - const Point &p, - const unsigned int component) const; - - /** - * Return the second derivative - * of the @p componentth vector - * component of the @p ith shape - * function at the point - * @p p. See the - * FiniteElement base - * class for more information - * about the semantics of this - * function. - */ - virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, - const Point &p, - const unsigned int component) const; - - /** - * Return the order - * of this finite element, - * i.e. the value passed to the - * constructor. - * - * Note that for this element, - * the order is actually one - * lower than the maximal - * polynomial degree, unlike for - * most of the other - * elements. For example, the RT0 - * element as piecewise linear - * shape functions, even though - * the normal component of them - * is piecewise constant on each - * face (the latter being the - * property that defines the - * order of the element). - */ - unsigned int get_degree () const; - - /** - * Return the matrix - * interpolating from the given - * finite element to the present - * one. The size of the matrix is - * then @p dofs_per_cell times - * source.dofs_per_cell. - * - * These matrices are only - * available if the source - * element is also a Raviart - * Thomas element. Otherwise, an - * exception of type - * FiniteElement::ExcInterpolationNotImplemented - * is thrown. - */ - virtual void - get_interpolation_matrix (const FiniteElement &source, - FullMatrix &matrix) const; /** * Number of base elements in a @@ -282,66 +198,18 @@ 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) - * of this object. - * - * This function is made virtual, - * since finite element objects - * are usually accessed through - * pointers to their base class, - * rather than the class itself. - */ virtual unsigned int memory_consumption () const; - - protected: - virtual FiniteElement * clone() const; - - 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; - private: /** * The order of the @@ -354,45 +222,6 @@ class FE_RaviartThomas : public FiniteElement */ const unsigned int rt_order; - /** - * Spaces describing the - * anisotropic polynomial spaces - * for each vector component, - * i.e. there are @p dim - * elements of this field. The - * values for this member are - * created in - * create_polynomials(). - */ - const std::vector > polynomials; - - /** - * For each shape function, store - * to which vector component (on - * the unit cell, they are mixed - * on the real cell by the - * transformation) they belong, - * and which index they have - * within the anisotropic tensor - * product polynomial space - * describing this vector - * component. - * - * These values are computed by - * the compute_renumber() - * function. - */ - const std::vector > renumber; - - - /** - * Generate the polynomial spaces - * for the polynomials() - * member. - */ - static std::vector > - create_polynomials (const unsigned int degree); - /** * Only for internal use. Its * full name is @@ -416,34 +245,6 @@ class FE_RaviartThomas : public FiniteElement static std::vector get_ria_vector (const unsigned int degree); - /** - * Compute the values of the - * @p renumber field. - */ - static std::vector > - compute_renumber (const unsigned int); - - /** - * Initialize the hanging node - * constraints matrices. Called - * from the constructor. - */ - void initialize_constraints (); - - /** - * Initialize the embedding - * matrices. Called from the - * constructor. - */ - void initialize_embedding (); - - /** - * Initialize the restriction - * matrices. Called from the - * constructor. - */ - void initialize_restriction (); - /** * Initialize the @p * generalized_support_points @@ -456,14 +257,6 @@ class FE_RaviartThomas : public FiniteElement */ void initialize_support_points (const unsigned int rt_degree); - /** - * Initialize the - * @p unit_face_support_points field - * of the FiniteElement - * class. Called from the - * constructor. - */ - void initialize_face_support_points (); /** * Given a set of flags indicating * what quantities are requested @@ -658,15 +451,6 @@ class FE_RaviartThomasNodal 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; virtual void interpolate(std::vector& local_dofs, const std::vector& values) const; virtual void interpolate(std::vector& local_dofs, @@ -715,70 +499,9 @@ class FE_RaviartThomasNodal #ifndef DOXYGEN -template <> -void FE_RaviartThomas<1>::initialize_face_support_points (); - template <> std::vector FE_RaviartThomas<1>::get_dpo_vector (const unsigned int); -template <> -std::vector > -FE_RaviartThomas<1>::create_polynomials (const unsigned int); - -template <> -std::vector > -FE_RaviartThomas<2>::create_polynomials (const unsigned int); - -template <> -std::vector > -FE_RaviartThomas<3>::create_polynomials (const unsigned int); - -template <> -std::vector > -FE_RaviartThomas<1>::compute_renumber (const unsigned int); - -template <> -std::vector > -FE_RaviartThomas<2>::compute_renumber (const unsigned int); - -template <> -std::vector > -FE_RaviartThomas<3>::compute_renumber (const unsigned int); - -template <> -void -FE_RaviartThomas<1>::initialize_constraints (); - -template <> -void -FE_RaviartThomas<2>::initialize_constraints (); - -template <> -void -FE_RaviartThomas<3>::initialize_constraints (); - -template <> -void -FE_RaviartThomas<1>::initialize_embedding (); - -template <> -void -FE_RaviartThomas<1>::initialize_restriction (); - -template <> -void -FE_RaviartThomas<2>::initialize_restriction (); - -template <> -void -FE_RaviartThomas<3>::initialize_restriction (); - -template <> -void -FE_RaviartThomas<1>:: -get_interpolation_matrix (const FiniteElement<1> &, - FullMatrix &) const; - #endif // DOXYGEN #endif 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 f69681b45a..e8c0513042 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -30,147 +30,66 @@ # include #endif - -// namespace for some functions that are used in this file. they are -// specific to numbering conventions used for the FE_RT element, and -// are thus not very interesting to the outside world -namespace -{ - // auxiliary type to allow for some - // kind of explicit template - // specialization of the following - // functions - template struct int2type {}; - - - // generate the j-th out of a total - // of N points on the unit square - // in 2d. N needs not be a square - // number, but must be the product - // of two integers - // - // there is one complication: we - // want to generate interpolation - // points on the unit square for - // the shape functions for this - // element, but for that we need to - // make sure that these - // interpolation points make the - // resulting matrix rows linearly - // independent. this is a problem - // since we have anisotropic - // polynomials, so for example for - // the lowest order elements, we - // have as polynomials in for the - // x-component of the shape - // functions only "x" and "1-x", - // i.e. no y-dependence. if we - // select as interpolation points - // the points (.5,0) and (.5,1), - // we're hosed! - // - // thus, the third parameter gives - // the coordinate direction in - // which the polynomial degree is - // highest. we use this to select - // interpolation points primarily - // in this direction then - inline - Point<2> generate_unit_point (const unsigned int j, - const unsigned int N, - const unsigned int d, - const int2type<2> &) - { - Assert (d<2, ExcInternalError()); - - // factorize N int N1*N2. note - // that we always have N1<=N2, - // since the square root is - // rounded down - const unsigned int N1 = static_cast(std::sqrt(1.*N)); - const unsigned int N2 = N/N1; - Assert (N1*N2 == N, ExcInternalError()); - - const unsigned int Nx = (d==0 ? N2 : N1), - Ny = (d==1 ? N2 : N1); - - return Point<2> (Nx == 1 ? .5 : 1.*(j%Nx)/(Nx-1), - Ny == 1 ? .5 : 1.*(j/Nx)/(Ny-1)); - } - - - // generate the j-th out of a total - // of N points on the unit cube - // in 3d. N needs not be a cube - // number, but must be the product - // of three integers - // - // the same applies as above for - // the meaning of the parameter "d" - inline - Point<3> generate_unit_point (const unsigned int /*j*/, - const unsigned int N, - const unsigned int d, - const int2type<3> &) - { - Assert (d<3, ExcInternalError()); - - const unsigned int N1 = static_cast(std::pow(1.*N, 1./3.)); - const unsigned int N2 = static_cast(std::sqrt(1.*N/N1)); - const unsigned int N3 = N/(N1*N2); - Assert (N1*N2*N3 == N, ExcInternalError()); - - Assert (false, ExcNotImplemented()); - - return Point<3> (); - } - -} - +#include +using namespace std; template -FE_RaviartThomas::FE_RaviartThomas (const unsigned int rt_order) +FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) : - FiniteElement (FiniteElementData(get_dpo_vector(rt_order), - dim, rt_order+1, FiniteElementData::Hdiv), - get_ria_vector (rt_order), - std::vector >(FiniteElementData(get_dpo_vector(rt_order),dim,rt_order+1).dofs_per_cell, - std::vector(dim,true))), - rt_order(rt_order), - polynomials (create_polynomials(rt_order)), - renumber (compute_renumber(rt_order)) + FE_PolyTensor, dim> ( + deg, + FiniteElementData(get_dpo_vector(deg), + dim, deg+1, FiniteElementData::Hdiv), + get_ria_vector (deg), + std::vector >( + FiniteElementData(get_dpo_vector(deg), + dim,deg+1).dofs_per_cell, + std::vector(dim,true))), + rt_order(deg) { Assert (dim >= 2, ExcImpossibleInDim(dim)); - - // check formula (III.3.22) in the - // book by Brezzi & Fortin about - // the number of degrees of freedom - // per cell - Assert (((dim==2) && - (this->dofs_per_cell == 2*(rt_order+1)*(rt_order+2))) - || - ((dim==3) && - (this->dofs_per_cell == 3*(rt_order+1)*(rt_order+1)*(rt_order+2))), - ExcInternalError()); - Assert (renumber.size() == this->dofs_per_cell, - ExcInternalError()); + const unsigned int n_dofs = this->dofs_per_cell; + + // First, initialize the + // generalized support points and + // quadrature weights, since they + // are required for interpolation. + initialize_support_points(deg); + // Now compute the inverse node + //matrix, generating the correct + //basis functions from the raw + //ones. + 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(M); + // From now on, the shape functions + // will be the correct ones, not + // the raw shape functions anymore. + // initialize the various matrices - initialize_constraints (); - for (unsigned int i=0; i::children_per_cell; ++i) - this->prolongation[i].reinit (this->dofs_per_cell, - this->dofs_per_cell); + this->prolongation[i].reinit (n_dofs, + n_dofs); FETools::compute_embedding_matrices (*this, &this->prolongation[0]); - initialize_restriction (); - - // finally fill in support points - // on cell and face - initialize_support_points (rt_order); - initialize_face_support_points (); - + std::vector > + face_embeddings(1<<(dim-1), FullMatrix(this->dofs_per_face, + this->dofs_per_face)); + FETools::compute_face_embedding_matrices(*this, &face_embeddings[0], 0, 0); + this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face, + this->dofs_per_face); + unsigned int target_row=0; + for (unsigned int d=0;dinterface_constraints(target_row,j) = face_embeddings[d](i,j); + ++target_row; + } +//TODO:[WB] What is this? // then make // system_to_component_table // invalid, since this has no @@ -217,599 +136,11 @@ FE_RaviartThomas::clone() const } -template -double -FE_RaviartThomas::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)); - - // check whether this shape - // function has a contribution in - // this component at all, and if so - // delegate to the respective - // polynomial - if (component == renumber[i].first) - return polynomials[component].compute_value(renumber[i].second, p); - else - return 0; -} - - - -template -Tensor<1,dim> -FE_RaviartThomas::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)); - - // check whether this shape - // function has a contribution in - // this component at all, and if so - // delegate to the respective - // polynomial - if (component == renumber[i].first) - return polynomials[component].compute_grad(renumber[i].second, p); - else - return Tensor<1,dim>(); -} - - - -template -Tensor<2,dim> -FE_RaviartThomas::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)); - - // check whether this shape - // function has a contribution in - // this component at all, and if so - // delegate to the respective - // polynomial - if (component == renumber[i].first) - return polynomials[component].compute_grad_grad(renumber[i].second, p); - else - return Tensor<2,dim>(); -} - - - -#if deal_II_dimension == 1 - -template <> -void -FE_RaviartThomas<1>:: -get_interpolation_matrix (const FiniteElement<1> &, - FullMatrix &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} - -#endif - - -template -void -FE_RaviartThomas:: -get_interpolation_matrix (const FiniteElement &x_source_fe, - FullMatrix &interpolation_matrix) const -{ - // this is only implemented, if the - // source FE is also a - // Raviart-Thomas element, - // otherwise throw an exception, as - // the documentation says - AssertThrow ((x_source_fe.get_name().find ("FE_RaviartThomas<") == 0) - || - (dynamic_cast*>(&x_source_fe) != 0), - typename FiniteElement:: - ExcInterpolationNotImplemented()); - - // ok, source is a RT element, so - // we will be able to do the work - const FE_RaviartThomas &source_fe - = dynamic_cast&>(x_source_fe); - - Assert (interpolation_matrix.m() == this->dofs_per_cell, - ExcDimensionMismatch (interpolation_matrix.m(), - this->dofs_per_cell)); - Assert (interpolation_matrix.n() == source_fe.dofs_per_cell, - ExcDimensionMismatch (interpolation_matrix.m(), - source_fe.dofs_per_cell)); - - - // compute the interpolation - // matrices in much the same way as - // we do for the embedding matrices - // from mother to child. - const unsigned int dofs_per_coordinate = this->dofs_per_cell/dim; - Assert (dofs_per_coordinate*dim == this->dofs_per_cell, - ExcInternalError()); - for (unsigned int d=0; d cell_interpolation (dofs_per_coordinate, - dofs_per_coordinate); - FullMatrix source_interpolation (dofs_per_coordinate, - source_dofs_per_coordinate); - FullMatrix tmp (dofs_per_coordinate, - source_dofs_per_coordinate); - for (unsigned int d=0; d p = generate_unit_point (j, dofs_per_coordinate, - d, int2type()); - for (unsigned int i=0; idofs_per_cell; ++i) - if (renumber[i].first == d) - for (unsigned int j=0; j 1e-15) - interpolation_matrix(i,j) = tmp(renumber[i].second, - source_fe.renumber[j].second); - } - - // if this were a Lagrange - // interpolation element, we could - // make sure that the row sum of - // each of the matrices is 1 at - // this point. note that this won't - // work here, since we are working - // with hierarchical elements for - // which the shape functions don't - // sum up to 1 - // - // however, we can make sure that - // only components couple that have - // the same vector component - for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; j -void -FE_RaviartThomas<1>::initialize_constraints () -{ - Assert (false, ExcImpossibleInDim(1)); -} - -#endif - -#if deal_II_dimension == 2 - -template <> -void -FE_RaviartThomas<2>::initialize_constraints () -{ - const unsigned int dim = 2; - - this->interface_constraints. - TableBase<2,double>::reinit (this->interface_constraints_size()); - - // this case is too easy, so - // special case it - if (rt_order == 0) - { - this->interface_constraints(0,0) = this->interface_constraints(1,0) = .5; - return; - } - - // for higher orders of the - // Raviart-Thomas element: - - // restricted to each face, the - // normal component of the shape - // functions is an element of P_{k} - // (in 2d), or Q_{k} (in 3d), where - // k is the degree of the element - // - // from this, we interpolate - // between mother and cell - // face. this is slightly - // complicated by the fact that we - // don't use Lagrange interpolation - // polynomials, but rather - // hierarchical polynomials, so we - // can't just use point - // interpolation. what we do - // instead is to evaluate at a - // number of points and then invert - // the interpolation matrix - - // mathematically speaking, this - // works in the following way: on - // each subface, we want that - // finite element solututions from - // both sides coincide. i.e. if a - // and b are expansion coefficients - // for the shape functions from - // both sides, we seek a relation - // between x and y such that - // sum_i a_i phi^c_i(x) - // == sum_j b_j phi_j(x) - // for all points x on the - // interface. here, phi^c_i are the - // shape functions on the small - // cell on one side of the face, - // and phi_j those on the big cell - // on the other side. To get this - // relation, it suffices to look at - // a sufficient number of points - // for which this has to hold. if - // there are n functions, then we - // need n evaluation points, and we - // choose them equidistantly. - // - // what one then gets is a matrix - // system - // a A == b B - // where - // A_ij = phi^c_i(x_j) - // B_ij = phi_i(x_j) - // and the relation we are looking for - // is - // a = (A^T)^-1 B^T b - // - // below, we build up these - // matrices, but rather than - // transposing them after the - // fact, we do so while building - // them. A will be - // subface_interpolation, B will be - // face_interpolation. note that we - // build up these matrices for all - // faces at once, rather than - // considering them separately. the - // reason is that we finally will - // want to have them in this order - // anyway, as this is the format we - // need inside deal.II - const std::vector > - face_polynomials (Polynomials::Hierarchical:: - generate_complete_basis (rt_order)); - Assert (face_polynomials.size() == this->dofs_per_face, ExcInternalError()); - - FullMatrix face_interpolation (2*this->dofs_per_face, this->dofs_per_face); - FullMatrix subface_interpolation (2*this->dofs_per_face, 2*this->dofs_per_face); - - // generate the matrix for the - // evaluation points on the big - // face, and the corresponding - // points in the coordinate system - // of the small face. order the - // shape functions in the same way - // we want to have them in the - // final matrix. extend shape - // functions on the small faces by - // zero to the other face on which - // they are not defined (we do this - // by simply not considering these - // entries in the matrix) - // - // note the agreeable fact that for - // this element, all the shape - // functions we presently care for - // are face-based (i.e. not vertex - // shape functions); thus, for this - // element, we can skip the - // annoying index shifting for the - // constraints matrix due to its - // weird format - for (unsigned int subface=0; subface::subfaces_per_face; ++subface) - for (unsigned int i=0; idofs_per_face; ++i) - { - const double p_face (1.*i/rt_order/2 + (subface == 0 ? 0. : .5)); - const double p_subface (1.*i/rt_order); - - for (unsigned int j=0; jdofs_per_face; ++j) - { - face_interpolation(subface*this->dofs_per_face+i, - j) - = face_polynomials[j].value(p_face); - subface_interpolation(subface*this->dofs_per_face+i, - subface*this->dofs_per_face+j) - = face_polynomials[j].value(p_subface); - } - } - - subface_interpolation.gauss_jordan (); - subface_interpolation.mmult (this->interface_constraints, - face_interpolation); - - // there is one additional thing to - // be considered: since the shape - // functions on the real cell - // contain the Jacobian (actually, - // the determinant of the inverse), - // there is an additional factor of - // 2 when going from the big to the - // small cell: - this->interface_constraints *= 1./2; - - // finally: constraints become - // really messy if the matrix in - // question has some entries that - // are almost zero, but not - // quite. this will happen in the - // above procedure due to - // round-off. let us simply delete - // these entries - for (unsigned int i=0; iinterface_constraints.m(); ++i) - for (unsigned int j=0; jinterface_constraints.n(); ++j) - if (std::fabs(this->interface_constraints(i,j)) < 1e-14) - this->interface_constraints(i,j) = 0.; -} - -#endif - -#if deal_II_dimension == 3 - -template <> -void -FE_RaviartThomas<3>::initialize_constraints () -{ - Assert (false, ExcNotImplemented()); -} - -#endif - - -#if deal_II_dimension == 1 - -template <> -void -FE_RaviartThomas<1>::initialize_restriction () -{} - -#endif - - -#if deal_II_dimension == 2 - -template <> -void -FE_RaviartThomas<2>::initialize_restriction () -{ - const unsigned int dim = 2; - switch (rt_order) - { - case 0: - { - // this is a strange element, - // since it is both additive - // and then it is also - // not. ideally, we would - // like to have the value of - // the shape function on the - // coarse line to be the mean - // value of that on the two - // child ones. thus, one - // should make it - // additive. however, - // additivity only works if - // an element does not have - // any continuity - // requirements, since - // otherwise degrees of - // freedom are shared between - // adjacent elements, and - // when we make the element - // additive, that would mean - // that we end up adding up - // contributions not only - // from the child cells of - // this cell, but also from - // the child cells of the - // neighbor, and since we - // cannot know whether there - // even exists a neighbor we - // cannot simply make the - // element additive. - // - // so, until someone comes - // along with a better - // alternative, we do the - // following: make the - // element non-additive, and - // simply pick the value of - // one of the child lines for - // the value of the mother - // line (note that we have to - // multiply by two, since the - // shape functions scale with - // the inverse Jacobian). we - // thus throw away the - // information of one of the - // child lines, but there - // seems to be no other way - // than that... - // - // note: to make things - // consistent, and - // restriction independent of - // the order in which we - // travel across the cells of - // the coarse grid, we have - // to make sure that we take - // the same small line when - // visiting its two - // neighbors, to get the - // value for the mother - // line. we take the first - // line always, in the - // canonical direction of - // lines - for (unsigned int c=0; c::children_per_cell; ++c) - this->restriction[c].reinit (this->dofs_per_cell, - this->dofs_per_cell); - - this->restriction[0](0,0) = 2.; - this->restriction[1](1,1) = 2.; - this->restriction[3](2,2) = 2.; - this->restriction[0](3,3) = 2.; - - break; - }; - - - case 1: - { - for (unsigned int c=0; c::children_per_cell; ++c) - this->restriction[c].reinit (this->dofs_per_cell, - this->dofs_per_cell); - - // first set the corner - // nodes. note that they are - // non-additive - this->restriction[0](0,0) = 2.; - this->restriction[0](6,6) = 2.; - - this->restriction[1](1,1) = 2.; - this->restriction[1](2,2) = 2.; - - this->restriction[2](3,3) = 2.; - this->restriction[2](5,5) = 2.; - - this->restriction[3](4,4) = 2.; - this->restriction[3](7,7) = 2.; - - // then also set the bubble - // nodes. they _are_ - // additive. to understand - // what's going on, recall - // that the bubble shape - // functions have value -1 - // (!) at the center point, - // by construction of the - // polynomials, and that the - // corner nodes have values - // 1/2 there since they are - // just the linears, and not - // some interpolating - // polynomial - // - // (actually, the - // additive/non-additive - // business shouldn't make - // that much of a difference: - // node 4 on cell 0 and node - // 0 on cell 3 must have the - // same value, since normal - // components are - // continuous. so we could - // pick either and make these - // shape functions - // non-additive as well. we - // choose to take the mean - // value, which should be the - // same as either value, and - // make the shape function - // additive) - this->restriction[0](10,0) = 1.; - this->restriction[0](10,4) = -1.; - this->restriction[3](10,0) = -1.; - this->restriction[3](10,4) = 1.; - - this->restriction[1](11,1) = 1.; - this->restriction[1](11,5) = -1.; - this->restriction[2](11,1) = -1.; - this->restriction[2](11,5) = 1.; - - this->restriction[0](8,6) = 1.; - this->restriction[0](8,2) = -1.; - this->restriction[1](8,6) = -1.; - this->restriction[1](8,2) = 1.; - - this->restriction[3](9,7) = 1.; - this->restriction[3](9,3) = -1.; - this->restriction[2](9,7) = -1.; - this->restriction[2](9,3) = 1.; - - break; - }; - - // in case we don't have the - // matrices (yet), leave them - // empty. this does not - // prevent the use of this FE, - // but will prevent the use of - // these matrices - }; -} - -#endif - -#if deal_II_dimension == 3 - -template <> -void -FE_RaviartThomas<3>::initialize_restriction () -{ - Assert (false, ExcNotImplemented()); -} - -#endif - - #if deal_II_dimension == 1 template @@ -873,8 +204,8 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) boundary_weights.reinit(n_face_points, legendre.n()); - Assert (face_points.n_quadrature_points == this->dofs_per_face, - ExcInternalError()); +// Assert (face_points.n_quadrature_points == this->dofs_per_face, +// ExcInternalError()); for (unsigned int k=0;k::initialize_support_points (const unsigned int deg) * legendre.compute_value(i, face_points.point(k)); } } - + Quadrature faces = QProjector::project_to_all_faces(face_points); for (;current::faces_per_cell*n_face_points; ++current) @@ -935,34 +266,6 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) #endif -#if deal_II_dimension == 1 - -template <> -void FE_RaviartThomas<1>::initialize_face_support_points () -{ - // no faces in 1d, so nothing to do -} - -#endif - - -template -void FE_RaviartThomas::initialize_face_support_points () -{ - this->unit_face_support_points.resize (this->dofs_per_face); - - // like with cell - // unit_support_points: - // associate all of the in - // the face mid-point, since - // there is no other useful - // way - for (unsigned int i=0; idofs_per_face; ++i) - this->unit_face_support_points[i] = (dim == 2 ? - Point(.5) : - Point(.5,.5)); -} - #if deal_II_dimension == 1 @@ -1055,263 +358,6 @@ FE_RaviartThomas::get_ria_vector (const unsigned int rt_order) } -#if deal_II_dimension == 1 - -template <> -std::vector > -FE_RaviartThomas<1>::create_polynomials (const unsigned int) -{ - Assert (false, ExcImpossibleInDim(1)); - return std::vector > (); -} - -#endif - - -#if deal_II_dimension == 2 - -template <> -std::vector > -FE_RaviartThomas<2>::create_polynomials (const unsigned int rt_order) -{ - const unsigned int dim = 2; - - // use the fact that the RT(k) - // spaces are spanned by the - // functions - // P_{k+1,k} \times P_{k,k+1}, - // see the book by Brezzi and - // Fortin - const std::vector > pols[2] - = { Polynomials::Hierarchical::generate_complete_basis (rt_order+1), - Polynomials::Hierarchical::generate_complete_basis (rt_order)}; - - // create spaces (k+1,k) and (k,k+1) - std::vector > > - pols_vector_1(dim), pols_vector_2(dim); - pols_vector_1[0] = pols[0]; - pols_vector_1[1] = pols[1]; - - pols_vector_2[0] = pols[1]; - pols_vector_2[1] = pols[0]; - - const AnisotropicPolynomials anisotropic[dim] - = { AnisotropicPolynomials (pols_vector_1), - AnisotropicPolynomials (pols_vector_2) }; - - // work around a stupid bug in - // gcc2.95 where the compiler - // complains about reaching the end - // of a non-void function when we - // simply return the following - // object unnamed, rather than - // first creating a named object - // and then returning it... - const std::vector > - ret_val (&anisotropic[0], &anisotropic[dim]); - return ret_val; -} - -#endif - - -#if deal_II_dimension == 3 - -template <> -std::vector > -FE_RaviartThomas<3>::create_polynomials (const unsigned int rt_order) -{ - const unsigned int dim = 3; - - // use the fact that the RT(k) - // spaces are spanned by the - // functions - // P_{k+1,k,k} \times P_{k,k+1,k} - // \times P_{k,k,k+1}, - // see the book by Brezzi and - // Fortin - const std::vector > pols[2] - = { Polynomials::Hierarchical::generate_complete_basis (rt_order+1), - Polynomials::Hierarchical::generate_complete_basis (rt_order)}; - - // create spaces (k+1,k,k), - // (k,k+1,k) and (k,k,k+1) - std::vector > > - pols_vector_1(dim), pols_vector_2(dim), pols_vector_3(dim); - pols_vector_1[0] = pols[0]; - pols_vector_1[1] = pols[1]; - pols_vector_1[2] = pols[1]; - - pols_vector_2[0] = pols[1]; - pols_vector_2[1] = pols[0]; - pols_vector_2[2] = pols[1]; - - pols_vector_3[0] = pols[1]; - pols_vector_3[1] = pols[1]; - pols_vector_3[2] = pols[0]; - - const AnisotropicPolynomials anisotropic[dim] - = { AnisotropicPolynomials (pols_vector_1), - AnisotropicPolynomials (pols_vector_2), - AnisotropicPolynomials (pols_vector_3) }; - - // work around a stupid bug in - // gcc2.95 where the compiler - // complains about reaching the end - // of a non-void function when we - // simply return the following - // object unnamed, rather than - // first creating a named object - // and then returning it... - const std::vector > - ret_val (&anisotropic[0], &anisotropic[dim]); - return ret_val; -} - -#endif - - - -#if deal_II_dimension == 1 - -template <> -std::vector > -FE_RaviartThomas<1>::compute_renumber (const unsigned int) -{ - Assert (false, ExcImpossibleInDim(1)); - return std::vector > (); -} - -#endif - - -#if deal_II_dimension == 2 - -template <> -std::vector > -FE_RaviartThomas<2>::compute_renumber (const unsigned int rt_order) -{ - const unsigned int dim = 2; - - std::vector > ret_val; - - // to explain the following: the - // first (rt_order+1) shape functions - // are on face 0, and point in - // y-direction, so are for the - // second vector component. then - // there are (rt_order+1) shape - // functions on face 1, which is - // for the x vector component, and - // so on. since the order of face - // rt_orders of freedom is arbitrary, - // we simply use the same order as - // that provided by the 1d - // polynomial class on which this - // element is based. after - // 4*(rt_order+1), the remaining - // shape functions are all bubbles, - // so we can number them in any way - // we want. we do so by first - // numbering the x-vectors, then - // the y-vectors - // - // now, we have to find a mapping - // from the above ordering to: - // first which vector component - // they belong to (easy), and - // second the index within this - // component as provided by the - // AnisotropicPolynomials class - // - // this is mostly a counting - // argument, tedious and error - // prone, and so boring to explain - // that we rather not try to do so - // here (it's simple, but boring, - // as said), aside from a few - // comments below - - // face 0 - for (unsigned int i=0; i test[dim] = { std::vector(ret_val.size()/dim, false), - std::vector(ret_val.size()/dim, false) }; - for (unsigned int i=0; i -std::vector > -FE_RaviartThomas<3>::compute_renumber (const unsigned int /*rt_order*/) -{ - Assert (false, ExcNotImplemented()); - return std::vector > (); -} - -#endif - - - template UpdateFlags @@ -1341,489 +387,10 @@ FE_RaviartThomas::update_each (const UpdateFlags flags) const return out; } - - //--------------------------------------------------------------------------- // Data field initialization //--------------------------------------------------------------------------- -template -typename Mapping::InternalDataBase * -FE_RaviartThomas::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; - - // initialize fields only if really - // necessary. otherwise, don't - // allocate memory - if (flags & update_values) - data->shape_values.resize (this->dofs_per_cell, - std::vector > (n_q_points)); - - if (flags & update_gradients) - data->shape_gradients.resize (this->dofs_per_cell, - std::vector > (n_q_points)); - - // if second derivatives through - // finite differencing is required, - // then initialize some objects for - // that - if (flags & update_second_derivatives) - data->initialize_2nd (this, mapping, quadrature); - - // next already fill those fields - // of which we have information by - // now. note that the shape values - // and gradients are only those on - // the unit cell, and need to be - // transformed when visiting an - // actual cell - for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int q=0; qshape_values[i][q][c] - = shape_value_component(i,quadrature.point(q),c); - - if (flags & update_gradients) - for (unsigned int c=0; cshape_gradients[i][q][c] - = shape_grad_component(i,quadrature.point(q),c); - } - - return data; -} - - - - -//--------------------------------------------------------------------------- -// Fill data of FEValues -//--------------------------------------------------------------------------- - -template -void -FE_RaviartThomas::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); - - // get the flags indicating the - // fields that have to be filled - const UpdateFlags flags(fe_data.current_update_flags()); - - const unsigned int n_q_points = quadrature.n_quadrature_points; - - // fill shape function - // values. these are vector-valued, - // so we have to transform - // them. since the output format - // (in data.shape_values) is a - // sequence of doubles (one for - // each non-zero shape function - // value, and for each quadrature - // point, rather than a sequence of - // small vectors, we have to use a - // number of conversions - if (flags & update_values) - { - std::vector > shape_values (n_q_points); - - Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim, - ExcInternalError()); - Assert (data.shape_values.n_cols() == n_q_points, - ExcInternalError()); - - for (unsigned int k=0; kdofs_per_cell; ++k) - { - // first transform shape - // values... - Assert (fe_data.shape_values[k].size() == n_q_points, - ExcInternalError()); - mapping.transform_covariant(fe_data.shape_values[k], 0, - shape_values, - mapping_data); - - // then copy over to target: - for (unsigned int q=0; q > shape_grads1 (n_q_points); - std::vector > shape_grads2 (n_q_points); - - Assert (data.shape_gradients.size() == this->dofs_per_cell * dim, - ExcInternalError()); - Assert (data.shape_gradients[0].size() == n_q_points, - ExcInternalError()); - - // loop over all shape - // functions, and treat the - // gradients of each shape - // function at all quadrature - // points - for (unsigned int k=0; kdofs_per_cell; ++k) - { - // treat the gradients of - // this particular shape - // function at all - // q-points. if Dv is the - // gradient of the shape - // function on the unit - // cell, then - // (J^-T)Dv(J^-1) is the - // value we want to have on - // the real cell. so, we - // will have to apply a - // covariant transformation - // to Dv twice. since the - // interface only allows - // multiplication with - // (J^-1) from the right, - // we have to trick a - // little in between - Assert (fe_data.shape_gradients[k].size() == n_q_points, - ExcInternalError()); - // do first transformation - mapping.transform_covariant(fe_data.shape_gradients[k], 0, - shape_grads1, - mapping_data); - // transpose matrix - for (unsigned int q=0; qcompute_2nd (mapping, cell, - QProjector::DataSetDescriptor::cell(), - mapping_data, fe_data, data); -} - - - -template -void -FE_RaviartThomas::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); - - // offset determines which data set - // to take (all data sets for all - // faces are stored contiguously) - const typename QProjector::DataSetDescriptor offset - = (QProjector::DataSetDescriptor:: - face (face, cell->face_orientation(face), - quadrature.n_quadrature_points)); - - // get the flags indicating the - // fields that have to be filled - const UpdateFlags flags(fe_data.current_update_flags()); - - const unsigned int n_q_points = quadrature.n_quadrature_points; - - // fill shape function - // values. these are vector-valued, - // so we have to transform - // them. since the output format - // (in data.shape_values) is a - // sequence of doubles (one for - // each non-zero shape function - // value, and for each quadrature - // point, rather than a sequence of - // small vectors, we have to use a - // number of conversions - if (flags & update_values) - { - Assert (fe_data.shape_values.size() == this->dofs_per_cell, - ExcInternalError()); - Assert (fe_data.shape_values[0].size() == - GeometryInfo::faces_per_cell * n_q_points * - (dim == 3 ? 2 : 1), - ExcInternalError()); - - std::vector > shape_values (n_q_points); - - Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim, - ExcInternalError()); - Assert (data.shape_values.n_cols() == n_q_points, - ExcInternalError()); - - for (unsigned int k=0; kdofs_per_cell; ++k) - { - // first transform shape - // values... - mapping.transform_covariant(fe_data.shape_values[k], offset, - shape_values, - mapping_data); - - // then copy over to target: - for (unsigned int q=0; qdofs_per_cell, - ExcInternalError()); - Assert (fe_data.shape_gradients[0].size() == - GeometryInfo::faces_per_cell * n_q_points * - (dim == 3 ? 2 : 1), - ExcInternalError()); - - std::vector > shape_grads1 (n_q_points); - std::vector > shape_grads2 (n_q_points); - - Assert (data.shape_gradients.size() == this->dofs_per_cell * dim, - ExcInternalError()); - Assert (data.shape_gradients[0].size() == n_q_points, - ExcInternalError()); - - // loop over all shape - // functions, and treat the - // gradients of each shape - // function at all quadrature - // points - for (unsigned int k=0; kdofs_per_cell; ++k) - { - // treat the gradients of - // this particular shape - // function at all - // q-points. if Dv is the - // gradient of the shape - // function on the unit - // cell, then - // (J^-T)Dv(J^-1) is the - // value we want to have on - // the real cell. so, we - // will have to apply a - // covariant transformation - // to Dv twice. since the - // interface only allows - // multiplication with - // (J^-1) from the right, - // we have to trick a - // little in between - // - // do first transformation - mapping.transform_covariant(fe_data.shape_gradients[k], offset, - shape_grads1, - mapping_data); - // transpose matrix - for (unsigned int q=0; qcompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); -} - - - -template -void -FE_RaviartThomas::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); - - // offset determines which data set - // to take (all data sets for all - // faces are stored contiguously) - const typename QProjector::DataSetDescriptor offset - = (QProjector::DataSetDescriptor:: - sub_face (face, subface, cell->face_orientation(face), - quadrature.n_quadrature_points)); - - // get the flags indicating the - // fields that have to be filled - const UpdateFlags flags(fe_data.current_update_flags()); - - const unsigned int n_q_points = quadrature.n_quadrature_points; - - // fill shape function - // values. these are vector-valued, - // so we have to transform - // them. since the output format - // (in data.shape_values) is a - // sequence of doubles (one for - // each non-zero shape function - // value, and for each quadrature - // point, rather than a sequence of - // small vectors, we have to use a - // number of conversions - if (flags & update_values) - { - Assert (fe_data.shape_values[0].size() == - GeometryInfo::faces_per_cell * - GeometryInfo::subfaces_per_face * - n_q_points, - ExcInternalError()); - - std::vector > shape_values (n_q_points); - - Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim, - ExcInternalError()); - Assert (data.shape_values.n_cols() == n_q_points, - ExcInternalError()); - - for (unsigned int k=0; kdofs_per_cell; ++k) - { - // first transform shape - // values... - mapping.transform_covariant(fe_data.shape_values[k], offset, - shape_values, - mapping_data); - - // then copy over to target: - for (unsigned int q=0; q::faces_per_cell * - GeometryInfo::subfaces_per_face * - n_q_points, - ExcInternalError()); - - std::vector > shape_grads1 (n_q_points); - std::vector > shape_grads2 (n_q_points); - - Assert (data.shape_gradients.size() == this->dofs_per_cell * dim, - ExcInternalError()); - Assert (data.shape_gradients[0].size() == n_q_points, - ExcInternalError()); - - // loop over all shape - // functions, and treat the - // gradients of each shape - // function at all quadrature - // points - for (unsigned int k=0; kdofs_per_cell; ++k) - { - // treat the gradients of - // this particular shape - // function at all - // q-points. if Dv is the - // gradient of the shape - // function on the unit - // cell, then - // (J^-T)Dv(J^-1) is the - // value we want to have on - // the real cell. so, we - // will have to apply a - // covariant transformation - // to Dv twice. since the - // interface only allows - // multiplication with - // (J^-1) from the right, - // we have to trick a - // little in between - // - // do first transformation - mapping.transform_covariant(fe_data.shape_gradients[k], offset, - shape_grads1, - mapping_data); - // transpose matrix - for (unsigned int q=0; qcompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); -} @@ -1935,9 +502,17 @@ FE_RaviartThomas::interpolate( for (unsigned int k=0;k::unit_normal_direction[face]); + local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i) + * values[face*n_face_points+k](GeometryInfo::unit_normal_direction[face]+offset); } + + const unsigned start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; + const unsigned start_cell_points = GeometryInfo::faces_per_cell*n_face_points; + + for (unsigned int k=0;k::interpolate( for (unsigned int k=0;kdofs_per_face] += boundary_weights(k,i) * values[GeometryInfo::unit_normal_direction[face]][face*n_face_points+k]; } + + const unsigned start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; + const unsigned start_cell_points = GeometryInfo::faces_per_cell*n_face_points; + + for (unsigned int k=0;k::memory_consumption () const } - -template -unsigned int -FE_RaviartThomas::get_degree () const -{ - return rt_order; -} - - - template class FE_RaviartThomas; 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 3c978fa9d2..a786989222 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 @@ -29,8 +29,6 @@ # include #endif -#include - template FE_RaviartThomasNodal::FE_RaviartThomasNodal (const unsigned int deg) : @@ -55,7 +53,6 @@ FE_RaviartThomasNodal::FE_RaviartThomasNodal (const unsigned int deg) // Set up the generalized support // points initialize_support_points (deg); - //Now compute the inverse node //matrix, generating the correct //basis functions from the raw @@ -123,15 +120,15 @@ FE_RaviartThomasNodal::get_name () const } - template -bool -FE_RaviartThomasNodal::has_support_on_face (unsigned int, unsigned int) const +FiniteElement * +FE_RaviartThomasNodal::clone() const { - return true; + return new FE_RaviartThomasNodal(this->degree-1); } + template void FE_RaviartThomasNodal::interpolate( @@ -238,14 +235,6 @@ FE_RaviartThomasNodal::interpolate( -template -FiniteElement * -FE_RaviartThomasNodal::clone() const -{ - return new FE_RaviartThomasNodal(this->degree-1); -} - - #if deal_II_dimension == 1 template <> -- 2.39.5