From 4160007e96c9078ca16e65480806b3c1c6ed7326 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 9 Mar 2018 00:25:36 +0100 Subject: [PATCH] Extend Hermite-like interpolation to case degree=2. --- include/deal.II/base/polynomial.h | 16 +++--- source/base/polynomial.cc | 88 +++++++++++++++++++++---------- source/fe/fe_dgq.cc | 12 ++--- 3 files changed, 73 insertions(+), 43 deletions(-) diff --git a/include/deal.II/base/polynomial.h b/include/deal.II/base/polynomial.h index c78532899f..3fdee6c474 100644 --- a/include/deal.II/base/polynomial.h +++ b/include/deal.II/base/polynomial.h @@ -591,9 +591,7 @@ namespace Polynomials /** * Polynomials for a variant of Hermite polynomials with better condition - * number in the interpolation than the basis from - * HermiteInterpolation. This class is only implemented for degree at least - * three, $n\geq 3$. + * number in the interpolation than the basis from HermiteInterpolation. * * In analogy to the actual Hermite polynomials this basis evaluates the * first polynomial $p_0$ to 1 at $x=0$ and has both a zero value and zero @@ -616,7 +614,14 @@ namespace Polynomials * Lagrange polynomials with double roots at $x=0$ and $x=1$. For example at * $n=4$, all of $p_0, p_1, p_3, p_4$ get an additional root at $x=0.5$ * through the factor $(x-0.5)$. - + * + * The basis only contains Hermite information at degree>=3, + * but it is also implemented for degrees between 0 and two. For the linear + * case, the usual hat functions are implemented, whereas the polynomials + * for degree=2 are $p_0(x)=(1-x)^2$, $p_1(x)=4x(x-1)$, and + * $p_2(x)=x^2$, in accordance with the construction principle for degree 3 + * that allows a non-zero of $p_0$ and $p_2$. + * * These two relaxations improve the condition number of the mass matrix * (i.e., interpolation) significantly, as can be seen from the following * table: @@ -696,8 +701,7 @@ namespace Polynomials /** * Return the polynomials with index 0 up to degree+1 in - * a space of degree up to degree. Here, degree has to - * be at least 3. + * a space of degree up to degree. */ static std::vector > generate_complete_basis (const unsigned int degree); diff --git a/source/base/polynomial.cc b/source/base/polynomial.cc index c495c34d91..3fb971418b 100644 --- a/source/base/polynomial.cc +++ b/source/base/polynomial.cc @@ -1286,9 +1286,6 @@ namespace Polynomials : Polynomial(0) { - Assert(degree>=3, - ExcNotImplemented("Hermite interpolation makes no sense for " - "degrees less than three")); AssertIndexRange(index, degree+1); this->coefficients.clear(); @@ -1296,16 +1293,52 @@ namespace Polynomials this->lagrange_support_points.resize(degree); - // 4 Polynomials with degree 3 - // entries (1,0) and (3,2) of the mass matrix will be equal to 0 - // - // | x 0 x x | - // | 0 x x x | - // M = | x x x 0 | - // | x x 0 x | - // - if (degree==3) + if (degree == 0) + this->lagrange_weight = 1.; + else if (degree == 1) { + if (index == 0) + { + this->lagrange_support_points[0] = 1.; + this->lagrange_weight = -1.; + } + else + { + this->lagrange_support_points[0] = 0.; + this->lagrange_weight = 1.; + } + } + else if (degree == 2) + { + if (index == 0) + { + this->lagrange_support_points[0] = 1.; + this->lagrange_support_points[1] = 1.; + this->lagrange_weight = 1.; + } + else if (index == 1) + { + this->lagrange_support_points[0] = 0; + this->lagrange_support_points[1] = 1; + this->lagrange_weight = 4.; + } + else + { + this->lagrange_support_points[0] = 0.; + this->lagrange_support_points[1] = 0.; + this->lagrange_weight = 1.; + } + } + else if (degree==3) + { + // 4 Polynomials with degree 3 + // entries (1,0) and (3,2) of the mass matrix will be equal to 0 + // + // | x 0 x x | + // | 0 x x x | + // M = | x x x 0 | + // | x x 0 x | + // if (index==0) { this->lagrange_support_points[0] = 2./7.; @@ -1335,23 +1368,22 @@ namespace Polynomials this->lagrange_weight = 3.5; } } - - // Higher order Polynomials degree>=4: the entries (1,0) and - // (degree,degree-1) of the mass matrix will be equal to 0 - // - // | x 0 x x x x x | - // | 0 x x x . . . x x x | - // | x x x x x x x | - // | x x x x x x x | - // | . . . | - // M = | . . . | - // | . . . | - // | x x x x x x x | - // | x x x x . . . x x 0 | - // | x x x x x 0 x | - // - if (degree >= 4) + else { + // Higher order Polynomials degree>=4: the entries (1,0) and + // (degree,degree-1) of the mass matrix will be equal to 0 + // + // | x 0 x x x x x | + // | 0 x x x . . . x x x | + // | x x x x x x x | + // | x x x x x x x | + // | . . . | + // M = | . . . | + // | . . . | + // | x x x x x x x | + // | x x x x . . . x x 0 | + // | x x x x x 0 x | + // We find the inner points as the zeros of the Jacobi polynomials // with alpha = beta = 2 which is the polynomial with the kernel // (1-x)^2 (1+x)^2, the two polynomials achieving zero value and zero diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc index 9abded35b8..321944f493 100644 --- a/source/fe/fe_dgq.cc +++ b/source/fe/fe_dgq.cc @@ -62,10 +62,8 @@ FE_DGQ::FE_DGQ (const unsigned int degree) { // Compute support points, which are the tensor product of the Lagrange // interpolation points in the constructor. - Quadrature support_quadrature(internal::FE_DGQ::get_QGaussLobatto_points(degree)); - Assert (support_quadrature.get_points().size() > 0, - (typename FiniteElement::ExcFEHasNoSupportPoints ())); - this->unit_support_points = support_quadrature.get_points(); + this->unit_support_points = + Quadrature(internal::FE_DGQ::get_QGaussLobatto_points(degree)).get_points(); // do not initialize embedding and restriction here. these matrices are // initialized on demand in get_restriction_matrix and @@ -877,11 +875,7 @@ FE_DGQLegendre::clone() const template FE_DGQHermite::FE_DGQHermite (const unsigned int degree) - : FE_DGQ(degree < 3 ? - Polynomials::generate_complete_Lagrange_basis - (internal::FE_DGQ::get_QGaussLobatto_points(degree)) - : - Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree)) + : FE_DGQ(Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree)) {} -- 2.39.5