From: kanschat Date: Wed, 12 Sep 2012 12:35:57 +0000 (+0000) Subject: introduce Hermite interpolation polynomials X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=25b8b18a6c460726ce5e8f7392c9dc61361aabbb;p=dealii-svn.git introduce Hermite interpolation polynomials git-svn-id: https://svn.dealii.org/trunk@26310 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/base/polynomial.h b/deal.II/include/deal.II/base/polynomial.h index 611b7a7884..3d92dc98c4 100644 --- a/deal.II/include/deal.II/base/polynomial.h +++ b/deal.II/include/deal.II/base/polynomial.h @@ -672,9 +672,64 @@ namespace Polynomials */ static std::vector > > recursive_coefficients; }; -} +/** + * Polynomials for Hermite interpolation condition. + * + * This is the set of polynomials of degree at least three, such that + * the following interpolation conditions are met: the polynomials and + * their first derivatives vanish at the values x=0 and + * x=1, with the exceptions p0(0)=1, + * p1(1)=1, p'2(0)=1, + * p'3(1)=1. + * + * For degree three, we obtain the standard four Hermitian + * interpolation polynomials, see for instance Wikipedia. + * For higher degrees, these are augmented + * first, by the polynomial of degree four with vanishing values and + * derivatives at x=0 and x=1, then by the product of + * this fourth order polynomial with Legendre polynomials of + * increasing order. The implementation is + * @f{align*}{ + * p_0(x) &= 2x^3-3x^2+1 \\ + * p_1(x) &= -2x^2+3x^2 \\ + * p_2(x) &= x^3-2x^2+x \\ + * p_3(x) &= x^3-x^2 \\ + * p_4(x) &= 16x^2(x-1)^2 \\ + * \ldots & \ldots \\ + * p_k(x) &= x^2(x-1)^2 L_{k-4}(x) + * @f} + * + * @author Guido Kanschat + * @date 2012 + */ + class HermiteInterpolation : public Polynomial + { + public: + /** + * Constructor for polynomial + * with index p. See + * the class documentation on + * the definition of the + * sequence of polynomials. + */ + HermiteInterpolation (const unsigned int p); + + /** + * Return the polynomials with index + * 0 up to + * p+1 in a space of + * degree up to + * p. Here, p + * has to be at least 3. + */ + static std::vector > + generate_complete_basis (const unsigned int p); + }; +} + /** @} */ diff --git a/deal.II/source/base/polynomial.cc b/deal.II/source/base/polynomial.cc index 652c8023c3..652d9d972c 100644 --- a/deal.II/source/base/polynomial.cc +++ b/deal.II/source/base/polynomial.cc @@ -1323,9 +1323,61 @@ std::vector > Lobatto::generate_complete_basis (const unsigne return v; } } -} + +// ------------------ HermiteInterpolation --------------- // + + HermiteInterpolation::HermiteInterpolation (const unsigned int p) + : + Polynomial((p<4) ? 3 : p+1) + { + if (p==0) + { + this->coefficients[0] = 1.; + this->coefficients[2] = -3.; + this->coefficients[3] = 2.; + } + else if (p==1) + { + this->coefficients[2] = 3.; + this->coefficients[3] = -2.; + } + else if (p==2) + { + this->coefficients[1] = 1.; + this->coefficients[2] = -2.; + this->coefficients[3] = 1.; + } + else if (p==3) + { + this->coefficients[2] = -1.; + this->coefficients[3] = 1.; + } + else + { + this->coefficients[4] = 16.; + this->coefficients[3] = -32.; + this->coefficients[2] = 16.; + + if (p>4) + { + Legendre legendre(p-4); + (*this) *= legendre; + } + } + } + std::vector > + HermiteInterpolation::generate_complete_basis (const unsigned int n) + { + std::vector > basis (n + 1); + + for (unsigned int i = 0; i <= n; ++i) + basis[i] = HermiteInterpolation (i); + + return basis; + } +} // ------------------ explicit instantiations --------------- //