From: Daniel Arndt Date: Wed, 12 Aug 2015 09:50:50 +0000 (-0500) Subject: Implement FE_Q_Bubbles that is based on TensorProductPolynomialsBubbles X-Git-Tag: v8.4.0-rc2~629^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=eafc440daf71643c862971f302624aaab1b88cbc;p=dealii.git Implement FE_Q_Bubbles that is based on TensorProductPolynomialsBubbles --- diff --git a/doc/doxygen/images/fe_q_bubbles_conditioning.png b/doc/doxygen/images/fe_q_bubbles_conditioning.png new file mode 100644 index 0000000000..a0618721d8 Binary files /dev/null and b/doc/doxygen/images/fe_q_bubbles_conditioning.png differ diff --git a/doc/news/changes.h b/doc/news/changes.h index 2db0318d82..2253787358 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -66,6 +66,11 @@ inconvenience this causes.
    +
  1. New: FE_Q_Bubbles describes a FiniteElement based on FE_Q + enriched by bubble functions. +
    + (Daniel Arndt, 2015/08/12) +
  2. New: MultithreadInfo::set_thread_limit() can now be called more than once and the environment variable DEAL_II_NUM_THREADS will be respected diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h index d93c690e58..ce2a7bf4d6 100644 --- a/include/deal.II/fe/fe_poly.templates.h +++ b/include/deal.II/fe/fe_poly.templates.h @@ -18,6 +18,7 @@ #include #include #include +#include #include #include diff --git a/include/deal.II/fe/fe_q_base.h b/include/deal.II/fe/fe_q_base.h index b89de91721..a5a9f83377 100644 --- a/include/deal.II/fe/fe_q_base.h +++ b/include/deal.II/fe/fe_q_base.h @@ -27,8 +27,8 @@ DEAL_II_NAMESPACE_OPEN /*@{*/ /** - * This class collects the basic methods used in FE_Q and FE_Q_DG0. There is - * no public constructor for this class as it is not functional as a stand- + * This class collects the basic methods used in FE_Q, FE_Q_DG0 and FE_Q_Bubbles. + * There is no public constructor for this class as it is not functional as a stand- * alone. The completion of definitions is left to the derived classes. * * @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann, @@ -324,6 +324,13 @@ private: * Mutex for protecting initialization of restriction and embedding matrix. */ mutable Threads::Mutex mutex; + + /* + * The highest polynomial degree of the underlying tensor product space + * without any enrichment. For FE_Q*(p) this is p. Note that enrichments + * may lead to a difference to degree. + */ + const unsigned int q_degree; }; diff --git a/include/deal.II/fe/fe_q_bubbles.h b/include/deal.II/fe/fe_q_bubbles.h new file mode 100644 index 0000000000..1e6b9cd015 --- /dev/null +++ b/include/deal.II/fe/fe_q_bubbles.h @@ -0,0 +1,205 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2012 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#ifndef __deal2__fe_q_bubbles_h +#define __deal2__fe_q_bubbles_h + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +/*!@addtogroup fe */ +/*@{*/ + +/** + * Implementation of a scalar Lagrange finite element @p Q_p^+ that yields the + * finite element space of continuous, piecewise polynomials of degree @p p in + * each coordinate direction plus some bubble enrichment space spanned by + * $(2x_j-1)^{p-1}\prod_{i=0}^{dim-1}(x_i(1-x_i))$. Therefore the highest polynomial + * degree is $p+1$. + * This class is realized using tensor product polynomials based on equidistant + * or given support points. + * + * The standard constructor of this class takes the degree @p p of this finite + * element. Alternatively, it can take a quadrature formula @p points defining + * the support points of the Lagrange interpolation in one coordinate direction. + * + * For more information about the spacedim template parameter check + * the documentation of FiniteElement or the one of Triangulation. + * + * Due to the fact that the enrichments are small almost everywhere + * for large p, the condition number for the mass and stiffness matrix fastly + * increaseses with increasing p. + * Below you see a comparison with FE_Q(QGaussLobatto(p+1)) for dim=1. + * + *

    + * @image html fe_q_bubbles_conditioning.png + *

    + * + * Therefore, this element should be used with care for $p>3$. + * + *

    Implementation

    + * + * The constructor creates a TensorProductPolynomials object that includes the + * tensor product of @p LagrangeEquidistant polynomials of degree @p p plus the + * bubble enrichments. This @p TensorProductPolynomialsBubbles object + * provides all values and derivatives of the shape functions. In case a + * quadrature rule is given, the constructor creates a + * TensorProductPolynomialsBubbles object that includes the tensor product of @p + * Lagrange polynomials with the support points from @p points and the bubble enrichments + * as defined above. + * + * Furthermore the constructor fills the @p interface_constrains, the @p + * prolongation (embedding) and the @p restriction matrices. + * + *

    Numbering of the degrees of freedom (DoFs)

    + * + * The original ordering of the shape functions represented by the + * TensorProductPolynomialsBubbles is a tensor product + * numbering. However, the shape functions on a cell are renumbered + * beginning with the shape functions whose support points are at the + * vertices, then on the line, on the quads, and finally (for 3d) on + * the hexes. Finally, there are support points for the bubble enrichments + * in the middle of the cell. + * + */ +template +class FE_Q_Bubbles : public FE_Q_Base,dim,spacedim> +{ +public: + /** + * Constructor for tensor product polynomials of degree @p p plus bubble enrichments + * + */ + FE_Q_Bubbles (const unsigned int p); + + /** + * Constructor for tensor product polynomials with support points @p points + * plus bubble enrichments based on a one-dimensional quadrature + * formula. The degree of the finite element is points.size(). + * Note that the first point has to be 0 and the last one 1. + */ + FE_Q_Bubbles (const Quadrature<1> &points); + + /** + * Return a string that uniquely identifies a finite element. This class + * returns FE_Q_Bubbles(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + virtual std::string get_name () const; + + /** + * Interpolate a set of scalar values, computed in the generalized support + * points. + */ + virtual void interpolate(std::vector &local_dofs, + const std::vector &values) const; + + /** + * Interpolate a set of vector values, computed in the generalized support + * points. + * + * Since a finite element often only interpolates part of a vector, + * offset is used to determine the first component of the vector to + * be interpolated. Maybe consider changing your data structures to use the + * next function. + */ + virtual void interpolate(std::vector &local_dofs, + const std::vector > &values, + unsigned int offset = 0) const; + + /** + * Interpolate a set of vector values, computed in the generalized support + * points. + */ + virtual void interpolate( + std::vector &local_dofs, + const VectorSlice > > &values) 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 @p + * FE_Q_Bubbles element. Otherwise, an exception of type + * FiniteElement::ExcInterpolationNotImplemented is thrown. + */ + virtual void + get_interpolation_matrix (const FiniteElement &source, + FullMatrix &matrix) const; + + virtual const FullMatrix & + get_prolongation_matrix (const unsigned int child, + const RefinementCase &refinement_case) const; + + virtual const FullMatrix & + get_restriction_matrix (const unsigned int child, + const RefinementCase &refinement_case) const; + + /** + * Check for non-zero values on a face. + * + * This function returns @p true, if the shape function @p shape_index has + * non-zero values on the face @p face_index. + * + * Implementation of the interface in FiniteElement + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const; + +protected: + /** + * @p clone function instead of a copy constructor. + * + * This function is needed by the constructors of @p FESystem. + */ + virtual FiniteElement *clone() const; + +private: + + /** + * Returns the restriction_is_additive flags. + * Only the last components for the bubble enrichments are true. + */ + static std::vector get_riaf_vector(const unsigned int degree); + + /** + * 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); + + /** + * Number of additional bubble functions + */ + const unsigned int n_bubbles; +}; + + + +/*@}*/ + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index b4130f5407..73fee37d1e 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -35,6 +35,7 @@ SET(_src fe_poly_tensor.cc fe_q_base.cc fe_q.cc + fe_q_bubbles.cc fe_q_dg0.cc fe_q_hierarchical.cc fe_q_iso_q1.cc @@ -72,6 +73,7 @@ SET(_inst fe_poly.inst.in fe_poly_tensor.inst.in fe_q_base.inst.in + fe_q_bubbles.inst.in fe_q_dg0.inst.in fe_q_hierarchical.inst.in fe_q.inst.in diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc index 147b4eb984..ebd25e7906 100644 --- a/source/fe/fe_poly.cc +++ b/source/fe/fe_poly.cc @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include diff --git a/source/fe/fe_poly.inst.in b/source/fe/fe_poly.inst.in index 6331952dd9..dcd33ed098 100644 --- a/source/fe/fe_poly.inst.in +++ b/source/fe/fe_poly.inst.in @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 1998 - 2014 by the deal.II authors +// Copyright (C) 1998 - 2015 by the deal.II authors // // This file is part of the deal.II library. // @@ -20,6 +20,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS #if deal_II_dimension <= deal_II_space_dimension template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; + template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly >, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 866290fac4..f6eb9fdbdc 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -124,6 +125,10 @@ struct FE_Q_Base::Implementation { const unsigned int dim = 2; + unsigned int q_deg = fe.degree; + if (types_are_equal >::value) + q_deg = fe.degree-1; + // restricted to each face, the traces 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. @@ -172,10 +177,10 @@ struct FE_Q_Base::Implementation // Add midpoint constraint_points.push_back (Point (0.5)); - if (fe.degree>1) + if (q_deg>1) { - const unsigned int n=fe.degree-1; - const double step=1./fe.degree; + const unsigned int n=q_deg-1; + const double step=1./q_deg; // subface 0 for (unsigned int i=1; i<=n; ++i) constraint_points.push_back ( @@ -197,13 +202,13 @@ struct FE_Q_Base::Implementation const std::vector &index_map_inverse = fe.poly_space.get_numbering_inverse(); const std::vector face_index_map = - FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(fe.degree); + FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point()) - 1.) < 1e-14, ExcInternalError()); for (unsigned int i=0; i p; p[0] = constraint_points[i](0); @@ -227,6 +232,10 @@ struct FE_Q_Base::Implementation { const unsigned int dim = 3; + unsigned int q_deg = fe.degree; + if (types_are_equal >::value) + q_deg = fe.degree-1; + // For a detailed documentation of the interpolation see the // FE_Q_Base<2>::initialize_constraints function. @@ -253,10 +262,10 @@ struct FE_Q_Base::Implementation constraint_points.push_back (Point (0.5, 0)); constraint_points.push_back (Point (0.5, 1)); - if (fe.degree>1) + if (q_deg>1) { - const unsigned int n=fe.degree-1; - const double step=1./fe.degree; + const unsigned int n=q_deg-1; + const double step=1./q_deg; std::vector > line_support_points(n); for (unsigned int i=0; i::Implementation // Now construct relation between destination (child) and source (mother) // dofs. - const unsigned int pnts=(fe.degree+1)*(fe.degree+1); + const unsigned int pnts=(q_deg+1)*(q_deg+1); // use that the element evaluates to 1 at index 0 and along the line at // zero const std::vector &index_map_inverse = fe.poly_space.get_numbering_inverse(); const std::vector face_index_map = - FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(fe.degree); + FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point()) - 1.) < 1e-14, ExcInternalError()); @@ -327,7 +336,7 @@ struct FE_Q_Base::Implementation for (unsigned int i=0; i constraint_point; @@ -373,14 +382,14 @@ struct FE_Q_Base::Implementation for (unsigned int j=0; j::FE_Q_Base (const POLY &poly_space, const std::vector &restriction_is_additive_flags) : FE_Poly(poly_space, fe_data, restriction_is_additive_flags, - std::vector(1, std::vector(1,true))) + std::vector(1, std::vector(1,true))), + q_degree (types_are_equal >::value + ?this->degree-1 + :this->degree) {} @@ -422,17 +434,18 @@ FE_Q_Base::initialize (const std::vector > &points) // distinguish q/q_dg0 case: need to be flexible enough to allow more // degrees of freedom than there are FE_Q degrees of freedom for derived // class FE_Q_DG0 that otherwise shares 95% of the code. - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + const unsigned int q_dofs_per_cell = Utilities::fixed_power(q_degree+1); Assert(q_dofs_per_cell == this->dofs_per_cell || - q_dofs_per_cell+1 == this->dofs_per_cell, ExcInternalError()); + q_dofs_per_cell+1 == this->dofs_per_cell || + q_dofs_per_cell+dim == this->dofs_per_cell , ExcInternalError()); { std::vector renumber(q_dofs_per_cell); - const FiniteElementData fe(get_dpo_vector(this->degree),1, - this->degree); + const FiniteElementData fe(get_dpo_vector(q_degree),1, + q_degree); FETools::hierarchic_to_lexicographic_numbering (fe, renumber); - if (this->dofs_per_cell > q_dofs_per_cell) - renumber.push_back(q_dofs_per_cell); + for (unsigned int i= q_dofs_per_cell; idofs_per_cell; ++i) + renumber.push_back(i); this->poly_space.set_numbering(renumber); } @@ -471,7 +484,7 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, x_source_fe.dofs_per_cell)); // only evaluate Q dofs - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + const unsigned int q_dofs_per_cell = Utilities::fixed_power(q_degree+1); const unsigned int source_q_dofs_per_cell = Utilities::fixed_power(source_fe->degree+1); // evaluation is simply done by evaluating the other FE's basis functions on @@ -504,7 +517,7 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, } // cut off very small values - const double eps = 2e-13*this->degree*dim; + const double eps = 2e-13*q_degree*dim; for (unsigned int i=0; idofs_per_cell; ++i) for (unsigned int j=0; jdofs_per_cell; ++j) if (std::fabs(interpolation_matrix(i,j)) < eps) @@ -596,7 +609,7 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe // Rule of thumb for FP accuracy, that can be expected for a given // polynomial degree. This value is used to cut off values close to // zero. - double eps = 2e-13*this->degree*(dim-1); + double eps = 2e-13*q_degree*(dim-1); // compute the interpolation matrix by simply taking the value at the // support points. @@ -888,11 +901,11 @@ void FE_Q_Base::initialize_unit_face_support_points (const st return; const unsigned int codim = dim-1; - this->unit_face_support_points.resize(Utilities::fixed_power(this->degree+1)); + this->unit_face_support_points.resize(Utilities::fixed_power(q_degree+1)); // find renumbering of faces and assign from values of quadrature std::vector face_index_map = - FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(this->degree); + FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(q_degree); Quadrature<1> support_1d(points); Quadrature support_quadrature(support_1d); this->unit_face_support_points.resize(support_quadrature.size()); @@ -914,7 +927,7 @@ FE_Q_Base::initialize_quad_dof_index_permutation () Assert (this->adjust_quad_dof_index_for_face_orientation_table.n_elements()==8*this->dofs_per_quad, ExcInternalError()); - const unsigned int n=this->degree-1; + const unsigned int n=q_degree-1; Assert(n*n==this->dofs_per_quad, ExcInternalError()); // alias for the table to fill @@ -1137,7 +1150,7 @@ FE_Q_Base return this->prolongation[refinement_case-1][child]; // distinguish q/q_dg0 case: only treat Q dofs first - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + const unsigned int q_dofs_per_cell = Utilities::fixed_power(q_degree+1); // compute the interpolation matrices in much the same way as we do for // the constraints. it's actually simpler here, since we don't have this @@ -1145,7 +1158,7 @@ FE_Q_Base // interpolation matrix is formed by a permutation of the indices of the // cell matrix. The value eps is used a threshold to decide when certain // evaluations of the Lagrange polynomials are zero or one. - const double eps = 1e-15*this->degree*dim; + const double eps = 1e-15*q_degree*dim; #ifdef DEBUG // in DEBUG mode, check that the evaluation of support points in the @@ -1167,7 +1180,7 @@ FE_Q_Base // the tensor product structure of this element and only evaluate 1D // information from the polynomial. This makes the cost of this function // almost negligible also for high order elements - const unsigned int dofs1d = this->degree+1; + const unsigned int dofs1d = q_degree+1; std::vector > subcell_evaluations (dim, Table<2,double>(dofs1d, dofs1d)); const std::vector &index_map_inverse = @@ -1315,7 +1328,7 @@ FE_Q_Base FullMatrix restriction(this->dofs_per_cell, this->dofs_per_cell); // distinguish q/q_dg0 case - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + const unsigned int q_dofs_per_cell = Utilities::fixed_power(q_degree+1); // for these Lagrange interpolation polynomials, construction of the // restriction matrices is relatively simple. the reason is that the @@ -1339,11 +1352,11 @@ FE_Q_Base // (compute on one child) by the same value (compute on a later child), // so we don't have to care about this - const double eps = 1e-15*this->degree*dim; + const double eps = 1e-15*q_degree*dim; const std::vector &index_map_inverse = this->poly_space.get_numbering_inverse(); - const unsigned int dofs1d = this->degree+1; + const unsigned int dofs1d = q_degree+1; std::vector > evaluations1d (dofs1d); restriction.reinit(this->dofs_per_cell, this->dofs_per_cell); @@ -1531,9 +1544,12 @@ std::pair, std::vector > FE_Q_Base::get_constant_modes () const { Table<2,bool> constant_modes(1, this->dofs_per_cell); - // FE_Q_DG0 should not use this function... - AssertDimension(this->dofs_per_cell, Utilities::fixed_power(this->degree+1)); - constant_modes.fill(true); + // We here just care for the constant mode due to the polynomial space + // without any enrichments + // As there may be more constant modes derived classes may to implement this + // themselves. An example for this is FE_Q_DG0. + for (unsigned int i=0; i(q_degree+1); ++i) + constant_modes(0, i) = true; return std::pair, std::vector > (constant_modes, std::vector(1, 0)); } diff --git a/source/fe/fe_q_base.inst.in b/source/fe/fe_q_base.inst.in index f0b13ce5c1..1e36e6bee3 100644 --- a/source/fe/fe_q_base.inst.in +++ b/source/fe/fe_q_base.inst.in @@ -20,6 +20,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS #if deal_II_dimension <= deal_II_space_dimension template class FE_Q_Base, deal_II_dimension, deal_II_space_dimension>; template class FE_Q_Base, deal_II_dimension, deal_II_space_dimension>; + template class FE_Q_Base, deal_II_dimension, deal_II_space_dimension>; template class FE_Q_Base >, deal_II_dimension, deal_II_space_dimension>; #endif } diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc new file mode 100644 index 0000000000..7616ac57a3 --- /dev/null +++ b/source/fe/fe_q_bubbles.cc @@ -0,0 +1,527 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2012 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace FE_Q_Bubbles_Helper +{ + namespace + { + template + inline + void + compute_embedding_matrices(const FE_Q_Bubbles &fe, + std::vector > > &matrices, + const bool isotropic_only) + { + const unsigned int dpc = fe.dofs_per_cell; + const unsigned int degree = fe.degree; + + // Initialize quadrature formula on fine cells + MappingQ1 mapping; + + Quadrature *q_fine; + Quadrature<1> q_dummy(std::vector >(1), std::vector (1,1.)); + switch (dim) + { + case 1: + if (spacedim==1) + q_fine = new QGauss (degree+1); + else if (spacedim==2) + q_fine = new QAnisotropic(QGauss<1>(degree+1), q_dummy); + else + q_fine = new QAnisotropic(QGauss<1>(degree+1), q_dummy, q_dummy); + break; + case 2: + if (spacedim==2) + q_fine = new QGauss (degree+1); + else + q_fine = new QAnisotropic(QGauss<1>(degree+1), QGauss<1>(degree+1), q_dummy); + break; + case 3: + q_fine = new QGauss (degree+1); + break; + default: + AssertThrow(false, ExcInternalError()); + } + + const unsigned int nq = q_fine->size(); + + // loop over all possible refinement cases + unsigned int ref_case = (isotropic_only) + ? RefinementCase::isotropic_refinement + : RefinementCase::cut_x; + for (; ref_case <= RefinementCase::isotropic_refinement; ++ref_case) + { + const unsigned int nc + = GeometryInfo::n_children(RefinementCase(ref_case)); + + for (unsigned int i=0; i tr; + GridGenerator::hyper_cube (tr, 0, 1); + tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); + tr.execute_coarsening_and_refinement(); + + DoFHandler dh(tr); + dh.distribute_dofs(fe); + + FEValues fine (mapping, fe, *q_fine, + update_quadrature_points + | update_JxW_values | update_values); + + const unsigned int n_dofs = dh.n_dofs(); + + FullMatrix fine_mass(n_dofs); + FullMatrix coarse_rhs_matrix(n_dofs, dpc); + + std::vector > child_ldi + (nc, std::vector(fe.dofs_per_cell)); + + //now create the mass matrix and all the right_hand sides + unsigned int child_no = 0; + typename dealii::DoFHandler::active_cell_iterator cell + = dh.begin_active(); + for (; cell!=dh.end(); ++cell, ++child_no) + { + fine.reinit(cell); + cell->get_dof_indices(child_ldi[child_no]); + + for (unsigned int q=0; q quad_tmp; + for (unsigned int k=0; k solution (n_dofs, dpc); + fine_mass.gauss_jordan(); + fine_mass.mmult(solution, coarse_rhs_matrix); + + //and distribute to the fine cell matrices + for (unsigned int child_no=0; child_no 1.e-12) + matrices[ref_case-1][child_no](i,j)=solution(gdi, j); + } + } + } + } +} + + +template +FE_Q_Bubbles::FE_Q_Bubbles (const unsigned int q_degree) + : + FE_Q_Base, dim, spacedim> ( + TensorProductPolynomialsBubbles(Polynomials::LagrangeEquidistant::generate_complete_basis(q_degree)), + FiniteElementData(get_dpo_vector(q_degree), + 1, q_degree+1, + FiniteElementData::H1), + get_riaf_vector(q_degree)), + n_bubbles((q_degree<=1)?1:dim) +{ + Assert (q_degree > 0, + ExcMessage ("This element can only be used for polynomial degrees " + "greater than zero")); + + std::vector > support_points_1d(q_degree+1); + for (unsigned int i=0; i<=q_degree; ++i) + support_points_1d[i][0] = static_cast(i)/q_degree; + + this->initialize(support_points_1d); + + // adjust unit support point for discontinuous node + Point point; + for (unsigned int d=0; dunit_support_points.push_back(point); + AssertDimension(this->dofs_per_cell, this->unit_support_points.size()); + + this->reinit_restriction_and_prolongation_matrices(); + if (dim == spacedim) + { + FE_Q_Bubbles_Helper::compute_embedding_matrices + (*this, this->prolongation, false); + // Fill restriction matrices with L2-projection + FETools::compute_projection_matrices (*this, this->restriction); + } + +} + + + +template +FE_Q_Bubbles::FE_Q_Bubbles (const Quadrature<1> &points) + : + FE_Q_Base, dim, spacedim> ( + TensorProductPolynomialsBubbles(Polynomials::generate_complete_Lagrange_basis(points.get_points())), + FiniteElementData(get_dpo_vector(points.size()-1), + 1, points.size(), + FiniteElementData::H1), + get_riaf_vector(points.size()-1)), + n_bubbles((points.size()-1<=1)?1:dim) +{ + const unsigned int q_degree = points.size()-1; + (void) q_degree; + Assert (q_degree > 0, + ExcMessage ("This element can only be used for polynomial degrees " + "at least one")); + + this->initialize(points.get_points()); + + // adjust unit support point for discontinuous node + Point point; + for (unsigned int d=0; dunit_support_points.push_back(point); + AssertDimension(this->dofs_per_cell, this->unit_support_points.size()); + + this->reinit_restriction_and_prolongation_matrices(); + if (dim == spacedim) + { + FE_Q_Bubbles_Helper::compute_embedding_matrices + (*this, this->prolongation, false); + // Fill restriction matrices with L2-projection + FETools::compute_projection_matrices (*this, this->restriction); + } +} + + + +template +std::string +FE_Q_Bubbles::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 + + std::ostringstream namebuf; + bool type = true; + const unsigned int n_points = this->degree; + std::vector points(n_points); + const unsigned int dofs_per_cell = this->dofs_per_cell; + const std::vector > &unit_support_points = this->unit_support_points; + unsigned int index = 0; + + // Decode the support points in one coordinate direction. + for (unsigned int j=0; j1) ? (unit_support_points[j](1)==0 && + ((dim>2) ? unit_support_points[j](2)==0: true)) : true) + { + if (index == 0) + points[index] = unit_support_points[j](0); + else if (index == 1) + points[n_points-1] = unit_support_points[j](0); + else + points[index-1] = unit_support_points[j](0); + + index++; + } + } + // Do not consider the discontinuous node for dimension 1 + Assert (index == n_points || (dim==1 && index == n_points+n_bubbles), + ExcMessage ("Could not decode support points in one coordinate direction.")); + + // Check whether the support points are equidistant. + for (unsigned int j=0; jdegree-1)) > 1e-15) + { + type = false; + break; + } + + if (type == true) + namebuf << "FE_Q_Bubbles<" << dim << ">(" << this->degree-1 << ")"; + else + { + + // Check whether the support points come from QGaussLobatto. + const QGaussLobatto<1> points_gl(n_points); + type = true; + for (unsigned int j=0; j(QGaussLobatto(" << this->degree << "))"; + else + namebuf << "FE_Q_Bubbles<" << dim << ">(QUnknownNodes(" << this->degree-1 << "))"; + } + return namebuf.str(); +} + + + +template +FiniteElement * +FE_Q_Bubbles::clone() const +{ + return new FE_Q_Bubbles(*this); +} + + + +template +void +FE_Q_Bubbles::interpolate(std::vector &local_dofs, + const std::vector &values) const +{ + Assert (values.size() == this->unit_support_points.size(), + ExcDimensionMismatch(values.size(), + this->unit_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + Assert (this->n_components() == 1, + ExcDimensionMismatch(this->n_components(), 1)); + + std::copy(values.begin(), values.end(), local_dofs.begin()); + + // We don't use the bubble functions for local interpolation + for (unsigned int i = 0; i +void +FE_Q_Bubbles::interpolate(std::vector &local_dofs, + const std::vector > &values, + unsigned int offset) const +{ + Assert (values.size() == this->unit_support_points.size(), + ExcDimensionMismatch(values.size(), + this->unit_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())); + + for (unsigned int i=0; idofs_per_cell-1; ++i) + { + const std::pair index + = this->system_to_component_index(i); + local_dofs[i] = values[i](offset+index.first); + } + + // We don't use the bubble functions for local interpolation + for (unsigned int i = 0; i +void +FE_Q_Bubbles::interpolate( + std::vector &local_dofs, + const VectorSlice > > &values) const +{ + Assert (values[0].size() == this->unit_support_points.size(), + ExcDimensionMismatch(values.size(), + this->unit_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + Assert (values.size() == this->n_components(), + ExcDimensionMismatch(values.size(), this->n_components())); + + for (unsigned int i=0; idofs_per_cell-1; ++i) + { + const std::pair index + = this->system_to_component_index(i); + local_dofs[i] = values[index.first][i]; + } + + // We don't use the bubble functions for local interpolation + for (unsigned int i = 0; i +void +FE_Q_Bubbles:: +get_interpolation_matrix (const FiniteElement &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + // We don't know how to do this properly, yet. + // However, for SolutionTransfer to work we need to provide an implementation + // for the case that the x_source_fe is identical to this FE + typedef FE_Q_Bubbles FEQBUBBLES; + typedef FiniteElement FEL; + + AssertThrow ((x_source_fe.get_name().find ("FE_Q_Bubbles<") == 0) + || + (dynamic_cast(&x_source_fe) != 0) + , + typename FEL:: + ExcInterpolationNotImplemented()); + Assert (interpolation_matrix.m() == this->dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + this->dofs_per_cell)); + Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + x_source_fe.dofs_per_cell)); + + //Provide a short cut in case we are just inquiring the identity + if (dynamic_cast(&x_source_fe)->degree == this->degree) + for (unsigned int i=0; i +std::vector +FE_Q_Bubbles::get_riaf_vector(const unsigned int q_deg) +{ + unsigned int n_cont_dofs = Utilities::fixed_power(q_deg+1); + const unsigned int n_bubbles = (q_deg<=1?1:dim); + return std::vector (n_cont_dofs+n_bubbles,true); +} + + + +template +std::vector +FE_Q_Bubbles::get_dpo_vector(const unsigned int q_deg) +{ + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i +bool +FE_Q_Bubbles::has_support_on_face +(const unsigned int shape_index, + const unsigned int face_index) const +{ + // discontinuous functions have no support on faces + if (shape_index >= this->n_dofs_per_cell()-n_bubbles) + return false; + else + return FE_Q_Base,dim,spacedim>::has_support_on_face(shape_index, face_index); +} + + + +template +const FullMatrix & +FE_Q_Bubbles::get_prolongation_matrix +(const unsigned int child, + const RefinementCase &refinement_case) const +{ + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Prolongation matrices are only available for refined cells!")); + Assert (child::n_children(refinement_case), + ExcIndexRange(child,0,GeometryInfo::n_children(refinement_case))); + + Assert (this->prolongation[refinement_case-1][child].n() != 0, + ExcMessage("This prolongation matrix has not been computed yet!")); + // finally return the matrix + return this->prolongation[refinement_case-1][child]; +} + + + +template +const FullMatrix & +FE_Q_Bubbles::get_restriction_matrix +(const unsigned int child, + const RefinementCase &refinement_case) const +{ + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Restriction matrices are only available for refined cells!")); + Assert (child::n_children(refinement_case), + ExcIndexRange(child,0,GeometryInfo::n_children(refinement_case))); + + Assert(this->restriction[refinement_case-1][child].n() != 0, + ExcMessage("This restriction matrix has not been computed yet!")); + + //finally return the matrix + return this->restriction[refinement_case-1][child]; +} + +// explicit instantiations +#include "fe_q_bubbles.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_q_bubbles.inst.in b/source/fe/fe_q_bubbles.inst.in new file mode 100644 index 0000000000..7cf26de8ed --- /dev/null +++ b/source/fe/fe_q_bubbles.inst.in @@ -0,0 +1,24 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template class FE_Q_Bubbles; +#endif + }