From b7dbfee287ea7642762d6ecbf77d2626df94ff5e Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 5 Jan 2021 10:13:23 -0700 Subject: [PATCH] Clarify a number of comments around FE_Q_Bubbles. --- .../base/tensor_product_polynomials_bubbles.h | 14 ++++-- include/deal.II/fe/fe_q_bubbles.h | 46 ++++++++++++------- .../tensor_product_polynomials_bubbles.cc | 8 ++-- source/fe/fe_q_bubbles.cc | 7 +-- 4 files changed, 50 insertions(+), 25 deletions(-) diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h index 7b45a6b9f8..a2042ac3d1 100644 --- a/include/deal.II/base/tensor_product_polynomials_bubbles.h +++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h @@ -36,10 +36,18 @@ DEAL_II_NAMESPACE_OPEN */ /** - * Tensor product of given polynomials and bubble functions of form - * $(2*x_j-1)^{degree-1}\prod_{i=0}^{dim-1}(x_i(1-x_i))$. This class inherits + * A class that represents a space of tensor product polynomials, augmented + * by $dim$ (non-normalized) bubble functions of form + * $\varphi_j(\mathbf x) + * = 2^{\text{degree}-1}\left(x_j-frac 12\right)^{\text{degree}-1} + * \left[\prod_{i=0}^{dim-1}(x_i(1-x_i))\right]$ + * for $j=0,\ldots,dim-1$. If `degree` is one, then the first factor + * disappears and one receives the usual bubble function centered + * at the mid-point of the cell. + * + * This class inherits * most of its functionality from TensorProductPolynomials. The bubble - * enrichments are added for the last indices. index. + * enrichments are added for the last index. */ template class TensorProductPolynomialsBubbles : public ScalarPolynomialsBase diff --git a/include/deal.II/fe/fe_q_bubbles.h b/include/deal.II/fe/fe_q_bubbles.h index 6ae5e84445..fd54bfd838 100644 --- a/include/deal.II/fe/fe_q_bubbles.h +++ b/include/deal.II/fe/fe_q_bubbles.h @@ -30,24 +30,30 @@ DEAL_II_NAMESPACE_OPEN /*@{*/ /** - * Implementation of a scalar Lagrange finite element @p Q_p^+ that yields the + * Implementation of a scalar Lagrange finite element $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. + * each coordinate direction plus some (non-normalized) bubble enrichment space + * spanned by the additional shape function + * $\varphi_j(\mathbf x) + * = 2^{p-1}\left(x_j-frac 12\right)^{p-1} + * \left[\prod_{i=0}^{dim-1}(x_i(1-x_i))\right]$. + * for $j=0,\ldots,dim-1$. If $p$ is one, then the first factor + * disappears and one receives the usual bubble function centered + * at the mid-point of the cell. + * Because these last shape functions have polynomial degree is $p+1$, the + * overall polynomial degree of the shape functions in the space described + * by this class is $p+1$. * - * 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. + * This class is realized using tensor product + * polynomials based on equidistant or given support points, in the same way as + * one can provide support points to the FE_Q class's constructors. * * For more information about the spacedim template parameter check - * the documentation of FiniteElement or the one of Triangulation. + * the documentation of the FiniteElement class, 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 + * $p$, the condition number for the mass and stiffness matrix quickly + * increaseses with increasing $p$. Below you see a comparison with * FE_Q(QGaussLobatto(p+1)) for dim=1. * *

@@ -56,6 +62,7 @@ DEAL_II_NAMESPACE_OPEN * * Therefore, this element should be used with care for $p>3$. * + * *

Implementation

* * The constructor creates a TensorProductPolynomials object that includes the @@ -70,6 +77,7 @@ DEAL_II_NAMESPACE_OPEN * 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 @@ -91,10 +99,16 @@ public: 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. + * 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 then + * points.size(), the plus one compared to the + * corresponding case for the FE_Q class coming from the additional + * bubble function. See the documentation of the FE_Q constructors + * for more information. + * + * Note that the first point has to be 0 + * and the last one 1. */ FE_Q_Bubbles(const Quadrature<1> &points); diff --git a/source/base/tensor_product_polynomials_bubbles.cc b/source/base/tensor_product_polynomials_bubbles.cc index 6ad931730f..a31abcc417 100644 --- a/source/base/tensor_product_polynomials_bubbles.cc +++ b/source/base/tensor_product_polynomials_bubbles.cc @@ -80,13 +80,15 @@ TensorProductPolynomialsBubbles::compute_value(const unsigned int i, const unsigned int comp = i - tensor_polys.n(); - // compute \prod_{i=1}^d 4*(1-x_i^2)(p) + // Compute \prod_{i=1}^d 4*x_i*(1-x_i) double value = 1.; for (unsigned int j = 0; j < dim; ++j) value *= 4 * p(j) * (1 - p(j)); - // and multiply with (2x_i-1)^{r-1} + + // Then multiply with (2x_i-1)^{r-1}. Since q_degree is generally a + // small integer, using a loop is likely faster than using std::pow. for (unsigned int i = 0; i < q_degree - 1; ++i) - value *= 2 * p(comp) - 1; + value *= (2 * p(comp) - 1); return value; } diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc index 3e06400cdf..3f2f9ab232 100644 --- a/source/fe/fe_q_bubbles.cc +++ b/source/fe/fe_q_bubbles.cc @@ -422,7 +422,7 @@ template 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_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); } @@ -437,8 +437,9 @@ FE_Q_Bubbles::get_dpo_vector(const unsigned int q_deg) for (unsigned int i = 1; i < dpo.size(); ++i) dpo[i] = dpo[i - 1] * (q_deg - 1); - dpo[dim] += - (q_deg <= 1 ? 1 : dim); // all the bubble functions are discontinuous + // Then add the bubble functions; they are all associated with the + // cell interior + dpo[dim] += (q_deg <= 1 ? 1 : dim); return dpo; } -- 2.39.5