]> https://gitweb.dealii.org/ - dealii.git/commitdiff
base: implement Hermite polynomials 14657/head
authorIvy Weber <ivy.weber@it.uu.se>
Tue, 10 Jan 2023 12:08:37 +0000 (13:08 +0100)
committerMatthias Maier <tamiko@43-1.org>
Fri, 9 Jun 2023 17:55:24 +0000 (12:55 -0500)
doc/doxygen/references.bib
include/deal.II/base/polynomials_hermite.h [new file with mode: 0644]
source/base/CMakeLists.txt
source/base/polynomials_hermite.cc [new file with mode: 0644]

index aaf6dde27945fa7d5a8d2489a5bf369a77cd9bac..78f904cd4d5630af1b9ffc71079a19bae4b9be5f 100644 (file)
   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 (file)
index 0000000..79a3341
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/subscriptor.h>
+
+#include <vector>
+
+
+
+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<double>
+  {
+  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<Polynomial<double>>
+    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 <code>side</code>.
+     */
+    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
index 1a134b2051dcbbbed87476a6c0198372cfedfe60..941800eb0b660c195932abcae5252ae556d64927 100644 (file)
@@ -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 (file)
index 0000000..c4c2dbb
--- /dev/null
@@ -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 <deal.II/base/polynomials_hermite.h>
+#include <deal.II/base/utilities.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Polynomials
+{
+  namespace
+  {
+    std::vector<double>
+    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<int>(regularity + 1 - curr_index);
+      const int loop_control_2 = (side == 1) ?
+                                   static_cast<int>(curr_index + 1) :
+                                   static_cast<int>(regularity + 2);
+
+      std::vector<double> 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<int>(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<int>(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<double>(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<Polynomial<double>>
+  PolynomialsHermite::generate_complete_basis(const unsigned int regularity)
+  {
+    std::vector<Polynomial<double>> 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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.