From 16c86939d43892cd5ea951e07d891e626791fcf7 Mon Sep 17 00:00:00 2001 From: Ivy Weber Date: Tue, 10 Jan 2023 13:08:37 +0100 Subject: [PATCH] base: implement Hermite polynomials --- doc/doxygen/references.bib | 9 ++ include/deal.II/base/polynomials_hermite.h | 152 +++++++++++++++++++++ source/base/CMakeLists.txt | 1 + source/base/polynomials_hermite.cc | 127 +++++++++++++++++ 4 files changed, 289 insertions(+) create mode 100644 include/deal.II/base/polynomials_hermite.h create mode 100644 source/base/polynomials_hermite.cc diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index aaf6dde279..78f904cd4d 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -2076,6 +2076,15 @@ url = {https://doi.org/10.11588/heidok.00005743} } +@article{CiarletRiavart1972interpolation, + author = {P. Ciarlet and P. Raviart}, + title = {General Langrange and Hermite interpolation in $\mathbb{R}^{n}$ with applications to finite element methods}, + journal = {Archive for rational mechanics and analysis}, + volume = {46}, + pages = {177-199}, + year = {1972} +} + @article{Chronopoulos1989, author = {A. T. Chronopoulos and C. W. Gear}, title = {S-step Iterative Methods for Symmetric Linear Systems}, diff --git a/include/deal.II/base/polynomials_hermite.h b/include/deal.II/base/polynomials_hermite.h new file mode 100644 index 0000000000..79a33410e8 --- /dev/null +++ b/include/deal.II/base/polynomials_hermite.h @@ -0,0 +1,152 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 - 2023 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +#ifndef dealii_polynomials_hermite_h +#define dealii_polynomials_hermite_h + +#include + +#include +#include +#include +#include + +#include + + + +DEAL_II_NAMESPACE_OPEN + +/** + * @addtogroup Polynomials + * @{ + */ + +namespace Polynomials +{ + /** + * This class implements Hermite interpolation polynomials (see + * @cite CiarletRiavart1972interpolation) enforcing the maximum + * possible level of regularity $r$ in the FEM basis given a + * polynomial degree of $2r+1$. The polynomials all represent + * either a non-zero shape value or derivative at $x=0$ and $x=1$ + * on the reference interval $x \in [0,1]$. + * + * Indices $j = 0, 1, \dots, r$ refer to polynomials corresponding + * to a non-zero derivative (or shape value for $j=0$) of + * order $j$ at $x=0$, and indices $j = r+1, r+2, \dots, 2r+1$ + * refer to polynomials with a non-zero derivative of order + * $j-(r+1)$ (or value for $j=r+1$) at $x=1$. In particular, the + * $0^{th}$ function has a value of $1$ at $x=0$, and the + * $(r+1)^{th}$ function has a value of $1$ at $x=1$.The basis is + * rescaled such that a function corresponding to a non-zero $j^{th}$ + * derivative has derivative value $j! 4^{j}$ at the corresponding + * node. This is done to prevent the $L^{2}$-norm of the basis functions + * from reducing exponentially with the chosen regularity. + */ + class PolynomialsHermite : public Polynomial + { + public: + /** + * Constructor for an individual Hermite polynomial. We write $f_{j}$ + * for a polynomial that has a non-zero $j^{th}$ derivative at $x=0$ + * and $g_{j}$ for a polynomial with a non-zero $j^{th}$ derivative + * at $x=1$, meaning $f_{j}$ will have @p index $=j$ and $g_{j}$ will + * have @p index $= j + \mathtt{regularity} + 1$. The resulting + * polynomials will be degree $2\times \mathtt{regularity} +1$ + * and obey the following conditions: + * @f{align*}{ + * &\begin{matrix} + * \left. \frac{d^{i}}{dx^{i}} f_{j}(x) \right\vert_{x=0} + * = i! 4^{i} \delta_{i, j}, \hfill + * &\qquad \hfill 0 \leq i \leq \mathtt{regularity}, \\ + * \left. \frac{d^{i}}{dx^{i}} f_{j}(x) \right\vert_{x=1} + * = 0, \hfill &\qquad \hfill 0 \leq i \leq \mathtt{regularity}, + * \end{matrix} \qquad 0 \leq j \leq \mathtt{regularity}, \\ + * &\begin{matrix} + * \left. \frac{d^{i}}{dx^{i}} g_{j}(x) \right\vert_{x=0} + * = 0, \hfill &\qquad \hfill 0 \leq i \leq \mathtt{regularity}, \\ + * \left. \frac{d^{i}}{dx^{i}} g_{j}(x) \right\vert_{x=1} + * = i! 4^{i} \delta_{i, j}, \hfill + * &\qquad \hfill 0 \leq i \leq \mathtt{regularity}, + * \end{matrix} \qquad 0 \leq j \leq \mathtt{regularity}, + * @f} + * where $\delta_{i,j}$ is equal to $1$ whenever $i=j$, + * and equal to $0$ otherwise. These polynomials have explicit + * formulas given by + * @f{align*}{ + * f_{j}(x) &= 4^{j} x^{j} (1-x)^{\mathtt{regularity}+1} + * \sum_{k=0}^{\mathtt{regularity} - j} \;^{\mathtt{regularity} + k} C_{k} + * x^{k}, \\ g_{j}(x) &= 4^{j} x^{\mathtt{regularity}+1} (x-1)^{j} + * \sum_{k=0}^{\mathtt{regularity} - j} \;^{\mathtt{regularity} + k} C_{k} + * (1-x)^{k}, + * @f} + * where $^{n} C_{r} = \frac{n!}{r!(n-r)!}$ is the $r^{th}$ binomial + * coefficient of degree $n, \; 0 \leq r \leq n$. + * + * @param regularity The highest derivative for which the basis + * is used to enforce regularity. + * @param index The local index of the generated polynomial in the + * Hermite basis. + */ + PolynomialsHermite(const unsigned int regularity, const unsigned int index); + + /** + * This function generates a vector of Polynomial objects + * representing a complete basis of degree $2\times\mathtt{regularity} +1$ + * on the reference interval $[0,1]$. + * + * @param regularity The generated basis can be used to strongly + * enforce continuity in all derivatives up to and including this + * order. + */ + static std::vector> + generate_complete_basis(const unsigned int regularity); + + protected: + /** + * Degree of the polynomial basis being used. + */ + unsigned int degree; + + /** + * The order of the highest derivative in which the Hermite + * basis can be used to impose continuity across element + * boundaries. It's related to the degree $p$ by + * $p = 2 \times\mathtt{regularity} +1$. + */ + unsigned int regularity; + + /** + * This variable stores the derivative that the shape function + * corresponds to at the element boundary given by side. + */ + unsigned int side_index; + + /** + * This stores whether the shape function corresponds to a non-zero + * value or derivative at $x=0$ on the reference interval + * ($\mathtt{side} =0$) or at $x=1$ ($\mathtt{side} =1$). + */ + unsigned int side; + }; +} // namespace Polynomials +/** @} */ + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index 1a134b2051..941800eb0b 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -71,6 +71,7 @@ set(_unity_include_src polynomials_bernstein.cc polynomials_bdm.cc polynomials_bernardi_raugel.cc + polynomials_hermite.cc polynomials_nedelec.cc polynomials_integrated_legendre_sz.cc polynomial_space.cc diff --git a/source/base/polynomials_hermite.cc b/source/base/polynomials_hermite.cc new file mode 100644 index 0000000000..c4c2dbb145 --- /dev/null +++ b/source/base/polynomials_hermite.cc @@ -0,0 +1,127 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 - 2023 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Polynomials +{ + namespace + { + std::vector + hermite_poly_coeffs(const unsigned int regularity, const unsigned int index) + { + AssertIndexRange(index, 2 * regularity + 2); + + const unsigned int curr_index = index % (regularity + 1); + const unsigned int side = (index > regularity) ? 1 : 0; + + // Signed ints are used here to protect against underflow errors + const int loop_control_1 = static_cast(regularity + 1 - curr_index); + const int loop_control_2 = (side == 1) ? + static_cast(curr_index + 1) : + static_cast(regularity + 2); + + std::vector poly_coeffs(2 * regularity + 2, 0.0); + + if (side == 1) // right side: g polynomials + { + int binomial_1 = (curr_index % 2) ? -1 : 1; + + for (int i = 0; i < loop_control_2; ++i) + { + int inv_binomial = 1; + + for (int j = 0; j < loop_control_1; ++j) + { + int binomial_2 = 1; + + for (int k = 0; k < j + 1; ++k) + { + poly_coeffs[regularity + i + k + 1] += + binomial_1 * inv_binomial * binomial_2; + binomial_2 *= k - j; + binomial_2 /= k + 1; + } + inv_binomial *= regularity + j + 1; + inv_binomial /= j + 1; + } + // ints used here to protect against underflow errors + binomial_1 *= -static_cast(curr_index - i); + binomial_1 /= i + 1; + } + } + else // left side: f polynomials + { + int binomial = 1; + + for (int i = 0; i < loop_control_2; ++i) + { + int inv_binomial = 1; + + for (int j = 0; j < loop_control_1; ++j) + { + poly_coeffs[curr_index + i + j] += binomial * inv_binomial; + inv_binomial *= regularity + j + 1; + inv_binomial /= j + 1; + } + // Protection needed here against underflow errors + binomial *= -static_cast(regularity + 1 - i); + binomial /= i + 1; + } + } + + // rescale coefficients by a factor of 4^curr_index to account for reduced + // L2-norms + double precond_factor = Utilities::pow(4, curr_index); + for (auto &it : poly_coeffs) + it *= precond_factor; + + return poly_coeffs; + } + } // namespace + + + + PolynomialsHermite::PolynomialsHermite(const unsigned int regularity, + const unsigned int index) + : Polynomial(hermite_poly_coeffs(regularity, index)) + , degree(2 * regularity + 1) + , regularity(regularity) + , side_index(index % (regularity + 1)) + , side((index >= regularity + 1) ? 1 : 0) + { + AssertIndexRange(index, 2 * (regularity + 1)); + } + + + + std::vector> + PolynomialsHermite::generate_complete_basis(const unsigned int regularity) + { + std::vector> polys; + const unsigned int sz = 2 * regularity + 2; + polys.reserve(sz); + + for (unsigned int i = 0; i < sz; ++i) + polys.emplace_back(PolynomialsHermite(regularity, i)); + + return polys; + } +} // namespace Polynomials + +DEAL_II_NAMESPACE_CLOSE -- 2.39.5