]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement Simplex::BarycentricPolynomials<dim>.
authorDavid Wells <drwells@email.unc.edu>
Thu, 14 Jan 2021 21:52:56 +0000 (16:52 -0500)
committerDavid Wells <drwells@email.unc.edu>
Thu, 21 Jan 2021 03:02:33 +0000 (22:02 -0500)
This gives us a much more convenient way to set up polynomials on simplices.

include/deal.II/simplex/barycentric_polynomials.h [new file with mode: 0644]
source/simplex/CMakeLists.txt
source/simplex/barycentric_polynomials.cc [new file with mode: 0644]
source/simplex/fe_lib.cc
tests/simplex/barycentric_01.cc [new file with mode: 0644]
tests/simplex/barycentric_01.output [new file with mode: 0644]

diff --git a/include/deal.II/simplex/barycentric_polynomials.h b/include/deal.II/simplex/barycentric_polynomials.h
new file mode 100644 (file)
index 0000000..6fe8626
--- /dev/null
@@ -0,0 +1,703 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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_simplex_barycentric_polynomials_h
+#define dealii_simplex_barycentric_polynomials_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/scalar_polynomials_base.h>
+#include <deal.II/base/table.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+  /**
+   * Polynomial implemented in barycentric coordinates.
+   *
+   * Barycentric coordinates are a coordinate system defined on simplices that
+   * are particularly easy to work with since they express coordinates in the
+   * simplex as convex combinations of the vertices. For example, any point in a
+   * triangle can be written as
+   *
+   * @f[
+   *   (x, y) = c_0 (x_0, y_0) + c_1 (x_1, y_1) + c_2 (x_2, y_2).
+   * @f]
+   *
+   * where each value $c_i$ is the relative weight of each vertex (so the
+   * centroid is, in 2D, where each $c_i = 1/3$). Since we only consider convex
+   * combinations we can rewrite this equation as
+   *
+   * @f[
+   *   (x, y) = (1 - c_1 - c_2) (x_0, y_0) + c_1 (x_1, y_1) + c_2 (x_2, y_2).
+   * @f]
+   *
+   * This results in three polynomials that are equivalent to $P^1$ in 2D. More
+   * exactly, this class implements a polynomial space defined with the basis,
+   * in 2D, of
+   * @f{align*}{
+   * t_0(x, y) &= 1 - x - y \\
+   * t_1(x, y) &= x \\
+   * t_2(x, y) &= y
+   * @f}
+   * and, in 3D,
+   * @f{align*}{
+   * t_0(x, y) &= 1 - x - y - z \\
+   * t_1(x, y) &= x             \\
+   * t_2(x, y) &= y             \\
+   * t_2(x, y) &= z
+   * @f}
+   *
+   * which is, in practice, a very convenient basis for defining simplex
+   * polynomials: for example, the fourth basis function of a TRI6 element is
+   *
+   * @f[
+   * 4 * t_1(x, y) * t_2(x, y).
+   * @f]
+   *
+   * Barycentric polynomials in <code>dim</code>-dimensional space have
+   * <code>dim + 1</code> variables in since <code>t_0</code> can be written in
+   * terms of the other monomials.
+   *
+   * Monomials can be conveniently constructed with
+   * BarycentricPolynomial::monomial().
+   *
+   * @ingroup Polynomials
+   */
+  template <int dim, typename Number = double>
+  class BarycentricPolynomial
+  {
+  public:
+    /**
+     * Constructor for the zero polynomial.
+     */
+    BarycentricPolynomial();
+
+    /**
+     * Constructor for a monomial.
+     */
+    BarycentricPolynomial(const TableIndices<dim + 1> &powers,
+                          const Number                 coefficient);
+
+    /**
+     * Return the specified monomial.
+     */
+    static BarycentricPolynomial<dim, Number>
+    monomial(const unsigned int d);
+
+    /**
+     * Print the polynomial to the output stream with lowest-order terms first.
+     * For example, the first P6 basis function is printed as
+     * <code>-1 * t0^1 + 2 * t0^2</code>, where <code>t0</code> is the first
+     * barycentric variable, <code>t1</code> is the second, etc.
+     */
+    void
+    print(std::ostream &out) const;
+
+    /**
+     * Degree of each barycentric polynomial.
+     */
+    TableIndices<dim + 1>
+    degrees() const;
+
+    /**
+     * Unary minus.
+     */
+    BarycentricPolynomial<dim, Number>
+    operator-() const;
+
+    /**
+     * Add a scalar.
+     */
+    template <typename Number2>
+    BarycentricPolynomial<dim, Number>
+    operator+(const Number2 &a) const;
+
+    /**
+     * Subtract a scalar.
+     */
+    template <typename Number2>
+    BarycentricPolynomial<dim, Number>
+    operator-(const Number2 &a) const;
+
+    /**
+     * Multiply by a scalar.
+     */
+    template <typename Number2>
+    BarycentricPolynomial<dim, Number> operator*(const Number2 &a) const;
+
+    /**
+     * Divide by a scalar.
+     */
+    template <typename Number2>
+    BarycentricPolynomial<dim, Number>
+    operator/(const Number2 &a) const;
+
+    /**
+     * Add another barycentric polynomial.
+     */
+    BarycentricPolynomial<dim, Number>
+    operator+(const BarycentricPolynomial<dim, Number> &augend) const;
+
+    /**
+     * Subtract another barycentric polynomial.
+     */
+    BarycentricPolynomial<dim, Number>
+    operator-(const BarycentricPolynomial<dim, Number> &augend) const;
+
+    /**
+     * Multiply by another barycentric polynomial.
+     */
+    BarycentricPolynomial<dim, Number>
+    operator*(const BarycentricPolynomial<dim, Number> &multiplicand) const;
+
+    /**
+     * Differentiate in barycentric coordinates.
+     */
+    BarycentricPolynomial<dim, Number>
+    barycentric_derivative(const unsigned int coordinate) const;
+
+    /**
+     * Differentiate in Cartesian coordinates.
+     */
+    BarycentricPolynomial<dim, Number>
+    derivative(const unsigned int coordinate) const;
+
+    /**
+     * Evaluate the polynomial.
+     */
+    Number
+    value(const Point<dim> &point) const;
+
+    /**
+     * Return an estimate, in bytes, of the memory usage of the object.
+     */
+    std::size_t
+    memory_consumption() const;
+
+  protected:
+    /**
+     * Coefficients of the polynomial. The exponents are the integer indexes.
+     */
+    Table<dim + 1, Number> coefficients;
+
+    /**
+     * Utility function for barycentric polynomials - its convenient to loop
+     * over all the indices at once in a dimension-independent way, but we also
+     * need to access the actual indices of the underlying Table object. This
+     * utility function converts an integral index into the equivalent
+     * TableIndices array (which are also the implicitly stored polynomial
+     * exponents).
+     */
+    static TableIndices<dim + 1>
+    index_to_indices(const std::size_t &          index,
+                     const TableIndices<dim + 1> &extent);
+  };
+
+  /**
+   * Scalar polynomial space based on barycentric polynomials.
+   */
+  template <int dim>
+  class BarycentricPolynomials : public ScalarPolynomialsBase<dim>
+  {
+  public:
+    /**
+     * Make the dimension available to the outside.
+     */
+    static const unsigned int dimension = dim;
+
+    /**
+     * Get the standard Lagrange basis for a specified degree.
+     */
+    static BarycentricPolynomials<dim>
+    get_fe_p_basis(const unsigned int degree);
+
+    /*
+     * Constructor taking the polynomial @p degree as input.
+     */
+    BarycentricPolynomials(
+      const std::vector<BarycentricPolynomial<dim>> &polynomials);
+
+    /**
+     * @copydoc ScalarPolynomialsBase::evaluate()
+     */
+    void
+    evaluate(const Point<dim> &           unit_point,
+             std::vector<double> &        values,
+             std::vector<Tensor<1, dim>> &grads,
+             std::vector<Tensor<2, dim>> &grad_grads,
+             std::vector<Tensor<3, dim>> &third_derivatives,
+             std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_value()
+     */
+    double
+    compute_value(const unsigned int i, const Point<dim> &p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_1st_derivative()
+     */
+    Tensor<1, dim>
+    compute_1st_derivative(const unsigned int i,
+                           const Point<dim> & p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_2nd_derivative()
+     */
+    Tensor<2, dim>
+    compute_2nd_derivative(const unsigned int i,
+                           const Point<dim> & p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_3rd_derivative()
+     */
+    Tensor<3, dim>
+    compute_3rd_derivative(const unsigned int i,
+                           const Point<dim> & p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_4th_derivative()
+     */
+    Tensor<4, dim>
+    compute_4th_derivative(const unsigned int i,
+                           const Point<dim> & p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_grad()
+     */
+    Tensor<1, dim>
+    compute_grad(const unsigned int i, const Point<dim> &p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::compute_grad_grad()
+     */
+    Tensor<2, dim>
+    compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::memory_consumption()
+     */
+    virtual std::size_t
+    memory_consumption() const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::name()
+     */
+    std::string
+    name() const override;
+
+    /**
+     * @copydoc ScalarPolynomialsBase::clone()
+     */
+    virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
+    clone() const override;
+
+  protected:
+    std::vector<BarycentricPolynomial<dim>> polys;
+
+    Table<2, BarycentricPolynomial<dim>> poly_grads;
+
+    Table<3, BarycentricPolynomial<dim>> poly_hessians;
+
+    Table<4, BarycentricPolynomial<dim>> poly_third_derivatives;
+
+    Table<5, BarycentricPolynomial<dim>> poly_fourth_derivatives;
+  };
+
+  // non-member template functions for algebra
+
+  /**
+   * Multiply a Simplex::BarycentricPolynomial by a constant.
+   */
+  template <int dim, typename Number1, typename Number2>
+  BarycentricPolynomial<dim, Number1>
+  operator*(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
+  {
+    return bp * Number1(a);
+  }
+
+  /**
+   * Add a constant to a Simplex::BarycentricPolynomial.
+   */
+  template <int dim, typename Number1, typename Number2>
+  BarycentricPolynomial<dim, Number1>
+  operator+(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
+  {
+    return bp + Number1(a);
+  }
+
+  /**
+   * Subtract a Simplex::BarycentricPolynomial from a constant.
+   */
+  template <int dim, typename Number1, typename Number2>
+  BarycentricPolynomial<dim, Number1>
+  operator-(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
+  {
+    return bp - Number1(a);
+  }
+
+  /**
+   * Write a Simplex::BarycentricPolynomial to the provided output stream.
+   */
+  template <int dim, typename Number>
+  std::ostream &
+  operator<<(std::ostream &out, const BarycentricPolynomial<dim, Number> &bp)
+  {
+    bp.print(out);
+    return out;
+  }
+} // namespace Simplex
+
+// Template function definitions
+
+namespace Simplex
+{
+  // BarycentricPolynomial:
+  template <int dim, typename Number>
+  BarycentricPolynomial<dim, Number>::BarycentricPolynomial()
+  {
+    TableIndices<dim + 1> extents;
+    for (unsigned int d = 0; d < dim + 1; ++d)
+      extents[d] = 1;
+    coefficients.reinit(extents);
+
+    coefficients(TableIndices<dim + 1>{}) = Number();
+  }
+
+
+
+  template <int dim, typename Number>
+  BarycentricPolynomial<dim, Number>::BarycentricPolynomial(
+    const TableIndices<dim + 1> &powers,
+    const Number                 coefficient)
+  {
+    TableIndices<dim + 1> extents;
+    for (unsigned int d = 0; d < dim + 1; ++d)
+      extents[d] = powers[d] + 1;
+    coefficients.reinit(extents);
+
+    coefficients(powers) = coefficient;
+  }
+
+
+
+  template <int dim, typename Number>
+  BarycentricPolynomial<dim, Number>
+  BarycentricPolynomial<dim, Number>::monomial(const unsigned int d)
+  {
+    AssertIndexRange(d, dim + 1);
+    TableIndices<dim + 1> indices;
+    indices[d] = 1;
+    return BarycentricPolynomial<dim, Number>(indices, Number(1));
+  }
+
+
+
+  template <int dim, typename Number>
+  void
+  BarycentricPolynomial<dim, Number>::print(std::ostream &out) const
+  {
+    const auto &coeffs     = this->coefficients;
+    auto        first      = index_to_indices(0, coeffs.size());
+    bool        print_plus = false;
+    if (coeffs(first) != Number())
+      {
+        out << coeffs(first);
+        print_plus = true;
+      }
+    for (std::size_t i = 1; i < coeffs.n_elements(); ++i)
+      {
+        const auto indices = index_to_indices(i, coeffs.size());
+        if (coeffs(indices) == Number())
+          continue;
+        if (print_plus)
+          out << " + ";
+        out << coeffs(indices);
+        for (unsigned int d = 0; d < dim + 1; ++d)
+          {
+            if (indices[d] != 0)
+              out << " * t" << d << '^' << indices[d];
+          }
+        print_plus = true;
+      }
+
+    if (!print_plus)
+      out << Number();
+  }
+
+
+
+  template <int dim, typename Number>
+  TableIndices<dim + 1>
+  BarycentricPolynomial<dim, Number>::degrees() const
+  {
+    auto deg = coefficients.size();
+    for (unsigned int d = 0; d < dim + 1; ++d)
+      deg[d] -= 1;
+    return deg;
+  }
+
+
+
+  template <int dim, typename Number>
+  BarycentricPolynomial<dim, Number>
+  BarycentricPolynomial<dim, Number>::operator-() const
+  {
+    return *this * Number(-1);
+  }
+
+
+
+  template <int dim, typename Number>
+  template <typename Number2>
+  BarycentricPolynomial<dim, Number>
+  BarycentricPolynomial<dim, Number>::operator+(const Number2 &a) const
+  {
+    BarycentricPolynomial<dim, Number> result(*this);
+    result.coefficients(index_to_indices(0, result.coefficients.size())) += a;
+
+    return result;
+  }
+
+
+
+  template <int dim, typename Number>
+  template <typename Number2>
+  BarycentricPolynomial<dim, Number>
+  BarycentricPolynomial<dim, Number>::operator-(const Number2 &a) const
+  {
+    return *this + (-a);
+  }
+
+
+
+  template <int dim, typename Number>
+  template <typename Number2>
+  BarycentricPolynomial<dim, Number> BarycentricPolynomial<dim, Number>::
+                                     operator*(const Number2 &a) const
+  {
+    if (a == Number2())
+      {
+        return BarycentricPolynomial<dim, Number>();
+      }
+
+    BarycentricPolynomial<dim, Number> result(*this);
+    for (std::size_t i = 0; i < result.coefficients.n_elements(); ++i)
+      {
+        const auto index = index_to_indices(i, result.coefficients.size());
+        result.coefficients(index) *= a;
+      }
+
+    return result;
+  }
+
+
+
+  template <int dim, typename Number>
+  template <typename Number2>
+  BarycentricPolynomial<dim, Number>
+  BarycentricPolynomial<dim, Number>::operator/(const Number2 &a) const
+  {
+    Assert(a != Number2(), ExcDivideByZero());
+    return *this * (Number(1) / Number(a));
+  }
+
+
+
+  template <int dim, typename Number>
+  BarycentricPolynomial<dim, Number>
+  BarycentricPolynomial<dim, Number>::
+  operator+(const BarycentricPolynomial<dim, Number> &augend) const
+  {
+    TableIndices<dim + 1> deg;
+    for (unsigned int d = 0; d < dim + 1; ++d)
+      {
+        deg[d] = std::max(degrees()[d], augend.degrees()[d]);
+      }
+
+    BarycentricPolynomial<dim, Number> result(deg, Number());
+
+    auto add_coefficients = [&](const Table<dim + 1, Number> &in) {
+      for (std::size_t i = 0; i < in.n_elements(); ++i)
+        {
+          const auto index = index_to_indices(i, in.size());
+          result.coefficients(index) += in(index);
+        }
+    };
+
+    add_coefficients(this->coefficients);
+    add_coefficients(augend.coefficients);
+    return result;
+  }
+
+
+
+  template <int dim, typename Number>
+  BarycentricPolynomial<dim, Number>
+  BarycentricPolynomial<dim, Number>::
+  operator-(const BarycentricPolynomial<dim, Number> &augend) const
+  {
+    return *this + (-augend);
+  }
+
+
+
+  template <int dim, typename Number>
+  BarycentricPolynomial<dim, Number> BarycentricPolynomial<dim, Number>::
+                                     operator*(const BarycentricPolynomial<dim, Number> &multiplicand) const
+  {
+    TableIndices<dim + 1> deg;
+    for (unsigned int d = 0; d < dim + 1; ++d)
+      {
+        deg[d] = multiplicand.degrees()[d] + degrees()[d];
+      }
+
+    BarycentricPolynomial<dim, Number> result(deg, Number());
+
+    const auto &coef_1   = this->coefficients;
+    const auto &coef_2   = multiplicand.coefficients;
+    auto &      coef_out = result.coefficients;
+
+    for (std::size_t i1 = 0; i1 < coef_1.n_elements(); ++i1)
+      {
+        const auto index_1 = index_to_indices(i1, coef_1.size());
+        for (std::size_t i2 = 0; i2 < coef_2.n_elements(); ++i2)
+          {
+            const auto index_2 = index_to_indices(i2, coef_2.size());
+
+            TableIndices<dim + 1> index_out;
+            for (unsigned int d = 0; d < dim + 1; ++d)
+              index_out[d] = index_1[d] + index_2[d];
+            coef_out(index_out) += coef_1(index_1) * coef_2(index_2);
+          }
+      }
+
+    return result;
+  }
+
+
+
+  template <int dim, typename Number>
+  BarycentricPolynomial<dim, Number>
+  BarycentricPolynomial<dim, Number>::barycentric_derivative(
+    const unsigned int coordinate) const
+  {
+    AssertIndexRange(coordinate, dim + 1);
+
+    if (degrees()[coordinate] == 0)
+      return BarycentricPolynomial<dim, Number>();
+
+    auto deg = degrees();
+    deg[coordinate] -= 1;
+    BarycentricPolynomial<dim, Number> result(
+      deg, std::numeric_limits<Number>::max());
+    const auto &coeffs_in  = coefficients;
+    auto &      coeffs_out = result.coefficients;
+    for (std::size_t i = 0; i < coeffs_out.n_elements(); ++i)
+      {
+        const auto out_index   = index_to_indices(i, coeffs_out.size());
+        auto       input_index = out_index;
+        input_index[coordinate] += 1;
+
+        coeffs_out(out_index) =
+          coeffs_in(input_index) * input_index[coordinate];
+      }
+
+    return result;
+  }
+
+
+
+  template <int dim, typename Number>
+  BarycentricPolynomial<dim, Number>
+  BarycentricPolynomial<dim, Number>::derivative(
+    const unsigned int coordinate) const
+  {
+    AssertIndexRange(coordinate, dim);
+    return -barycentric_derivative(0) + barycentric_derivative(coordinate + 1);
+  }
+
+
+
+  template <int dim, typename Number>
+  Number
+  BarycentricPolynomial<dim, Number>::value(const Point<dim> &point) const
+  {
+    // TODO: this is probably not numerically stable for higher order.
+    // We really need some version of Horner's method.
+    Number result = {};
+
+    // Begin by converting point (which is in Cartesian coordinates) to
+    // barycentric coordinates:
+    std::array<Number, dim + 1> b_point;
+    b_point[0] = 1.0;
+    for (unsigned int d = 0; d < dim; ++d)
+      {
+        b_point[0] -= point[d];
+        b_point[d + 1] = point[d];
+      }
+
+    // Now evaluate the polynomial at the computed barycentric point:
+    for (std::size_t i = 0; i < coefficients.n_elements(); ++i)
+      {
+        const auto indices = index_to_indices(i, coefficients.size());
+        const auto coef    = coefficients(indices);
+        if (coef == Number())
+          continue;
+
+        auto temp = Number(1);
+        for (unsigned int d = 0; d < dim + 1; ++d)
+          temp *= std::pow(b_point[d], indices[d]);
+        result += coef * temp;
+      }
+
+    return result;
+  }
+
+  template <int dim, typename Number>
+  std::size_t
+  BarycentricPolynomial<dim, Number>::memory_consumption() const
+  {
+    return coefficients.memory_consumption();
+  }
+
+  template <int dim, typename Number>
+  TableIndices<dim + 1>
+  BarycentricPolynomial<dim, Number>::index_to_indices(
+    const std::size_t &          index,
+    const TableIndices<dim + 1> &extent)
+  {
+    TableIndices<dim + 1> result;
+    auto                  temp = index;
+
+    for (unsigned int n = 0; n < dim + 1; ++n)
+      {
+        std::size_t slice_size = 1;
+        for (unsigned int n2 = n + 1; n2 < dim + 1; ++n2)
+          slice_size *= extent[n2];
+        result[n] = temp / slice_size;
+        temp %= slice_size;
+      }
+    return result;
+  }
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index f547691c5aae70c170468633a559cc4c9ab47cc2..398736607dda611e94920a390208405533f3ca9d 100644 (file)
@@ -18,6 +18,7 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 SET(_unity_include_src
   fe_lib.cc
   grid_generator.cc
+  barycentric_polynomials.cc
   polynomials.cc
   quadrature_lib.cc
   )
diff --git a/source/simplex/barycentric_polynomials.cc b/source/simplex/barycentric_polynomials.cc
new file mode 100644 (file)
index 0000000..8e4e541
--- /dev/null
@@ -0,0 +1,341 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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/simplex/barycentric_polynomials.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+  namespace internal
+  {
+    /**
+     * Get the highest degree of the barycentric polynomial (in Cartesian
+     * coordinates).
+     */
+    template <int dim>
+    unsigned int
+    get_degree(const std::vector<BarycentricPolynomial<dim>> &polys)
+    {
+      // Since the first variable in a simplex polynomial is, e.g., in 2D,
+      //
+      // t0 = 1 - x - y
+      //
+      // (that is, it depends on the Cartesian variables), we have to compute
+      // its degree separately. An example: t0*t1*t2 has degree 1 in the affine
+      // polynomial basis but is degree 2 in the Cartesian polynomial basis.
+      std::size_t max_degree = 0;
+      for (const auto &poly : polys)
+        {
+          const TableIndices<dim + 1> degrees = poly.degrees();
+
+          const auto  degree_0 = degrees[0];
+          std::size_t degree_d = 0;
+          for (unsigned int d = 1; d < dim + 1; ++d)
+            degree_d = std::max(degree_d, degrees[d]);
+
+          max_degree = std::max(max_degree, degree_d + degree_0);
+        }
+
+      return max_degree;
+    }
+  } // namespace internal
+
+
+  template <int dim>
+  BarycentricPolynomials<dim>
+  BarycentricPolynomials<dim>::get_fe_p_basis(const unsigned int degree)
+  {
+    std::vector<BarycentricPolynomial<dim>> polys;
+
+    auto M = [](const unsigned int d) {
+      return BarycentricPolynomial<dim, double>::monomial(d);
+    };
+    switch (degree)
+      {
+        case 0:
+          polys.push_back(0 * M(0) + 1);
+          break;
+        case 1:
+          {
+            for (unsigned int d = 0; d < dim + 1; ++d)
+              polys.push_back(M(d));
+            break;
+          }
+        case 2:
+          {
+            for (unsigned int d = 0; d < dim + 1; ++d)
+              polys.push_back(M(d) * (2 * M(d) - 1));
+            polys.push_back(4 * M(1) * M(0));
+            if (dim >= 2)
+              {
+                polys.push_back(4 * M(1) * M(2));
+                polys.push_back(4 * M(2) * M(0));
+              }
+            if (dim == 3)
+              {
+                polys.push_back(4 * M(3) * M(0));
+                polys.push_back(4 * M(1) * M(3));
+                polys.push_back(4 * M(2) * M(3));
+              }
+            break;
+          }
+        default:
+          Assert(false, ExcNotImplemented());
+      }
+
+    return BarycentricPolynomials<dim>(polys);
+  }
+
+
+
+  template <int dim>
+  BarycentricPolynomials<dim>::BarycentricPolynomials(
+    const std::vector<BarycentricPolynomial<dim>> &polynomials)
+    : ScalarPolynomialsBase<dim>(internal::get_degree(polynomials),
+                                 polynomials.size())
+  {
+    polys = polynomials;
+
+    poly_grads.reinit({polynomials.size(), dim});
+    poly_hessians.reinit({polynomials.size(), dim, dim});
+    poly_third_derivatives.reinit({polynomials.size(), dim, dim, dim});
+    poly_fourth_derivatives.reinit({polynomials.size(), dim, dim, dim, dim});
+
+    for (std::size_t i = 0; i < polynomials.size(); ++i)
+      {
+        // gradients
+        for (unsigned int d = 0; d < dim; ++d)
+          poly_grads[i][d] = polynomials[i].derivative(d);
+
+        // hessians
+        for (unsigned int d0 = 0; d0 < dim; ++d0)
+          for (unsigned int d1 = 0; d1 < dim; ++d1)
+            poly_hessians[i][d0][d1] = poly_grads[i][d0].derivative(d1);
+
+        // third derivatives
+        for (unsigned int d0 = 0; d0 < dim; ++d0)
+          for (unsigned int d1 = 0; d1 < dim; ++d1)
+            for (unsigned int d2 = 0; d2 < dim; ++d2)
+              poly_third_derivatives[i][d0][d1][d2] =
+                poly_hessians[i][d0][d1].derivative(d2);
+
+        // fourth derivatives
+        for (unsigned int d0 = 0; d0 < dim; ++d0)
+          for (unsigned int d1 = 0; d1 < dim; ++d1)
+            for (unsigned int d2 = 0; d2 < dim; ++d2)
+              for (unsigned int d3 = 0; d3 < dim; ++d3)
+                poly_fourth_derivatives[i][d0][d1][d2][d3] =
+                  poly_third_derivatives[i][d0][d1][d2].derivative(d3);
+      }
+  }
+
+
+
+  template <int dim>
+  void
+  BarycentricPolynomials<dim>::evaluate(
+    const Point<dim> &           unit_point,
+    std::vector<double> &        values,
+    std::vector<Tensor<1, dim>> &grads,
+    std::vector<Tensor<2, dim>> &grad_grads,
+    std::vector<Tensor<3, dim>> &third_derivatives,
+    std::vector<Tensor<4, dim>> &fourth_derivatives) const
+  {
+    Assert(values.size() == this->n() || values.size() == 0,
+           ExcDimensionMismatch2(values.size(), this->n(), 0));
+    Assert(grads.size() == this->n() || grads.size() == 0,
+           ExcDimensionMismatch2(grads.size(), this->n(), 0));
+    Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
+           ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
+    Assert(third_derivatives.size() == this->n() ||
+             third_derivatives.size() == 0,
+           ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
+    Assert(fourth_derivatives.size() == this->n() ||
+             fourth_derivatives.size() == 0,
+           ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
+
+    for (std::size_t i = 0; i < polys.size(); ++i)
+      {
+        if (values.size() == this->n())
+          values[i] = polys[i].value(unit_point);
+
+        // gradients
+        if (grads.size() == this->n())
+          for (unsigned int d = 0; d < dim; ++d)
+            grads[i][d] = poly_grads[i][d].value(unit_point);
+
+        // hessians
+        if (grad_grads.size() == this->n())
+          for (unsigned int d0 = 0; d0 < dim; ++d0)
+            for (unsigned int d1 = 0; d1 < dim; ++d1)
+              grad_grads[i][d0][d1] =
+                poly_hessians[i][d0][d1].value(unit_point);
+
+        // third derivatives
+        if (third_derivatives.size() == this->n())
+          for (unsigned int d0 = 0; d0 < dim; ++d0)
+            for (unsigned int d1 = 0; d1 < dim; ++d1)
+              for (unsigned int d2 = 0; d2 < dim; ++d2)
+                third_derivatives[i][d0][d1][d2] =
+                  poly_third_derivatives[i][d0][d1][d2].value(unit_point);
+
+        // fourth derivatives
+        if (fourth_derivatives.size() == this->n())
+          for (unsigned int d0 = 0; d0 < dim; ++d0)
+            for (unsigned int d1 = 0; d1 < dim; ++d1)
+              for (unsigned int d2 = 0; d2 < dim; ++d2)
+                for (unsigned int d3 = 0; d3 < dim; ++d3)
+                  fourth_derivatives[i][d0][d1][d2][d3] =
+                    poly_fourth_derivatives[i][d0][d1][d2][d3].value(
+                      unit_point);
+      }
+  }
+
+
+
+  template <int dim>
+  double
+  BarycentricPolynomials<dim>::compute_value(const unsigned int i,
+                                             const Point<dim> & p) const
+  {
+    AssertIndexRange(i, this->n());
+    return polys[i].value(p);
+  }
+
+
+
+  template <int dim>
+  Tensor<1, dim>
+  BarycentricPolynomials<dim>::compute_1st_derivative(const unsigned int i,
+                                                      const Point<dim> &p) const
+  {
+    Tensor<1, dim> result;
+    for (unsigned int d = 0; d < dim; ++d)
+      result[d] = poly_grads[i][d].value(p);
+    return result;
+  }
+
+
+
+  template <int dim>
+  Tensor<2, dim>
+  BarycentricPolynomials<dim>::compute_2nd_derivative(const unsigned int i,
+                                                      const Point<dim> &p) const
+  {
+    Tensor<2, dim> result;
+    for (unsigned int d0 = 0; d0 < dim; ++d0)
+      for (unsigned int d1 = 0; d1 < dim; ++d1)
+        result[d0][d1] = poly_hessians[i][d0][d1].value(p);
+
+    return result;
+  }
+
+
+
+  template <int dim>
+  Tensor<3, dim>
+  BarycentricPolynomials<dim>::compute_3rd_derivative(const unsigned int i,
+                                                      const Point<dim> &p) const
+  {
+    Tensor<3, dim> result;
+    for (unsigned int d0 = 0; d0 < dim; ++d0)
+      for (unsigned int d1 = 0; d1 < dim; ++d1)
+        for (unsigned int d2 = 0; d2 < dim; ++d2)
+          result[d0][d1][d2] = poly_third_derivatives[i][d0][d1][d2].value(p);
+
+    return result;
+  }
+
+
+
+  template <int dim>
+  Tensor<4, dim>
+  BarycentricPolynomials<dim>::compute_4th_derivative(const unsigned int i,
+                                                      const Point<dim> &p) const
+  {
+    Tensor<4, dim> result;
+    for (unsigned int d0 = 0; d0 < dim; ++d0)
+      for (unsigned int d1 = 0; d1 < dim; ++d1)
+        for (unsigned int d2 = 0; d2 < dim; ++d2)
+          for (unsigned int d3 = 0; d3 < dim; ++d3)
+            result[d0][d1][d2][d3] =
+              poly_fourth_derivatives[i][d0][d1][d2][d3].value(p);
+
+    return result;
+  }
+
+
+
+  template <int dim>
+  Tensor<1, dim>
+  BarycentricPolynomials<dim>::compute_grad(const unsigned int i,
+                                            const Point<dim> & p) const
+  {
+    return compute_1st_derivative(i, p);
+  }
+
+
+
+  template <int dim>
+  Tensor<2, dim>
+  BarycentricPolynomials<dim>::compute_grad_grad(const unsigned int i,
+                                                 const Point<dim> & p) const
+  {
+    return compute_2nd_derivative(i, p);
+  }
+
+
+
+  template <int dim>
+  std::unique_ptr<ScalarPolynomialsBase<dim>>
+  BarycentricPolynomials<dim>::clone() const
+  {
+    return std::make_unique<BarycentricPolynomials<dim>>(*this);
+  }
+
+
+
+  template <int dim>
+  std::string
+  BarycentricPolynomials<dim>::name() const
+  {
+    return "BarycentricPolynomials<" + std::to_string(dim) + ">";
+  }
+
+
+
+  template <int dim>
+  std::size_t
+  BarycentricPolynomials<dim>::memory_consumption() const
+  {
+    std::size_t poly_memory = 0;
+    for (const auto &poly : polys)
+      poly_memory += poly.memory_consumption();
+    return ScalarPolynomialsBase<dim>::memory_consumption() + poly_memory +
+           poly_grads.memory_consumption() +
+           poly_hessians.memory_consumption() +
+           poly_third_derivatives.memory_consumption() +
+           poly_fourth_derivatives.memory_consumption();
+  }
+
+  template class BarycentricPolynomials<1>;
+  template class BarycentricPolynomials<2>;
+  template class BarycentricPolynomials<3>;
+
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
index 83f965f66ceda771d124684c77498e85cc361642..d06188f983304f1f7ec5ffba8c28525f934a3600 100644 (file)
@@ -20,6 +20,7 @@
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_tools.h>
 
+#include <deal.II/simplex/barycentric_polynomials.h>
 #include <deal.II/simplex/fe_lib.h>
 
 DEAL_II_NAMESPACE_OPEN
diff --git a/tests/simplex/barycentric_01.cc b/tests/simplex/barycentric_01.cc
new file mode 100644 (file)
index 0000000..6bd73fb
--- /dev/null
@@ -0,0 +1,170 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+// Test Simplex::BarycentricPolynomial and Simplex::BarycentricPolynomials.
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/table.h>
+
+#include <deal.II/simplex/barycentric_polynomials.h>
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/polynomials.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+int
+main()
+{
+  initlog();
+
+  Simplex::BarycentricPolynomial<2> bp2({1, 0, 0}, 1.0);
+  deallog << bp2 << std::endl;
+
+  // test some basic algebra with barycentric polynomials
+  {
+    deallog << "1D:" << std::endl;
+    const auto bp1_0 = Simplex::BarycentricPolynomial<1>::monomial(0);
+    const auto bp1_1 = Simplex::BarycentricPolynomial<1>::monomial(1);
+
+    deallog << "bp1_0 = " << bp1_0 << std::endl;
+    deallog << "bp1_1 = " << bp1_1 << std::endl;
+    deallog << "bp1_0 * 2 * bp1_1 / 2 = " << bp1_0 * 2 * bp1_1 / 2 << std::endl
+            << std::endl;
+  }
+
+  {
+    deallog << std::endl << "2D:" << std::endl;
+    const auto bp2_0 = Simplex::BarycentricPolynomial<2>::monomial(0) * 2;
+    deallog << "bp2_0 = " << bp2_0 << std::endl;
+
+    const auto bp2_1 = 3.0 * Simplex::BarycentricPolynomial<2>::monomial(1);
+    deallog << "bp2_1 = " << bp2_1 << std::endl;
+
+    const auto bp2_2 = Simplex::BarycentricPolynomial<2>::monomial(2);
+    deallog << "bp2_2 = " << bp2_2 << std::endl;
+
+    const auto prod1 = bp2_0 + bp2_1;
+    deallog << "bp2_0 + bp2_1 = " << prod1 << std::endl;
+
+    const auto prod2 = prod1 * bp2_0;
+    deallog << "(bp2_0 + bp2_1) * bp2_0 = " << prod2 << std::endl;
+    deallog << "bp2_0 * bp2_0 + bp2_1 * bp2_0 = "
+            << bp2_0 * bp2_0 + bp2_1 * bp2_0 << std::endl;
+    deallog << "bp2_1 * bp2_0 + bp2_0 * bp2_0 = "
+            << bp2_1 * bp2_0 + bp2_0 * bp2_0 << std::endl;
+
+    // test derivatives
+    deallog << "d/dx bp2_0 = " << bp2_0.derivative(0) << std::endl;
+    deallog << "d/dy bp2_0 = " << bp2_0.derivative(1) << std::endl;
+
+    deallog << "d/dx bp2_2 = " << bp2_2.derivative(0) << std::endl;
+    deallog << "d/dy bp2_2 = " << bp2_2.derivative(1) << std::endl;
+  }
+
+  // test various finite element spaces
+  {
+    deallog << std::endl << "Test with TRI6" << std::endl;
+
+    const auto t1 = Simplex::BarycentricPolynomial<2>::monomial(0);
+    const auto t2 = Simplex::BarycentricPolynomial<2>::monomial(1);
+    const auto t3 = Simplex::BarycentricPolynomial<2>::monomial(2);
+
+    std::vector<Simplex::BarycentricPolynomial<2>> p2;
+    p2.push_back(t1 * (2 * t1 - 1));
+    p2.push_back(t2 * (2 * t2 - 1));
+    p2.push_back(t3 * (2 * t3 - 1));
+    p2.push_back(4 * t2 * t1);
+    p2.push_back(4 * t2 * t3);
+    p2.push_back(4 * t3 * t1);
+
+    Simplex::FE_P<2> fe(2);
+    for (unsigned int i = 0; i < 6; ++i)
+      {
+        deallog << "p = " << p2[i] << std::endl;
+        deallog << "p_x = " << p2[i].derivative(0) << std::endl;
+        deallog << "p_y = " << p2[i].derivative(1) << std::endl;
+        for (unsigned int j = 0; j < 6; ++j)
+          {
+            Assert(std::abs(p2[i].value(fe.get_unit_support_points()[j]) -
+                            double(i == j)) < 1e-12,
+                   ExcInternalError());
+          }
+        deallog << std::endl;
+      }
+  }
+
+  {
+    deallog << std::endl << "Test with TET4" << std::endl;
+    const auto tet4 = Simplex::BarycentricPolynomials<3>::get_fe_p_basis(1);
+
+    Simplex::FE_P<3> fe(1);
+    const auto &     points = fe.get_unit_support_points();
+    for (unsigned int i = 0; i < 4; ++i)
+      {
+        Assert(points.size() == 4, ExcInternalError());
+        for (unsigned int j = 0; j < 4; ++j)
+          {
+            Assert(std::abs(tet4.compute_value(i, points[j]) - double(i == j)) <
+                     1e-12,
+                   ExcInternalError());
+
+            // first derivatives should be constant
+            Assert((tet4.compute_grad(i, points[0]) -
+                    tet4.compute_grad(i, points[j]))
+                       .norm() == 0.0,
+                   ExcInternalError());
+            Assert(tet4.compute_2nd_derivative(i, points[j]).norm() == 0.0,
+                   ExcInternalError());
+            Assert(tet4.compute_3rd_derivative(i, points[j]).norm() == 0.0,
+                   ExcInternalError());
+            Assert(tet4.compute_4th_derivative(i, points[j]).norm() == 0.0,
+                   ExcInternalError());
+          }
+      }
+    deallog << "Test with TET4 - Success" << std::endl;
+  }
+
+  {
+    deallog << "Test with TET10" << std::endl;
+    const auto tet10 = Simplex::BarycentricPolynomials<3>::get_fe_p_basis(2);
+
+    Simplex::FE_P<3> fe(2);
+    const auto &     points = fe.get_unit_support_points();
+    for (unsigned int i = 0; i < 10; ++i)
+      {
+        Assert(points.size() == 10, ExcInternalError());
+        for (unsigned int j = 0; j < 10; ++j)
+          {
+            Assert(std::abs(tet10.compute_value(i, points[j]) -
+                            double(i == j)) < 1e-12,
+                   ExcInternalError());
+
+            // second derivatives should be constant
+            Assert((tet10.compute_2nd_derivative(i, points[0]) -
+                    tet10.compute_2nd_derivative(i, points[j]))
+                       .norm() == 0.0,
+                   ExcInternalError());
+
+            Assert(tet10.compute_3rd_derivative(i, points[j]).norm() == 0.0,
+                   ExcInternalError());
+            Assert(tet10.compute_4th_derivative(i, points[j]).norm() == 0.0,
+                   ExcInternalError());
+          }
+      }
+    deallog << "Test with TET10 - Success" << std::endl;
+  }
+}
diff --git a/tests/simplex/barycentric_01.output b/tests/simplex/barycentric_01.output
new file mode 100644 (file)
index 0000000..f083f6c
--- /dev/null
@@ -0,0 +1,51 @@
+
+DEAL::1.00000 * t0^1
+DEAL::1D:
+DEAL::bp1_0 = 1.00000 * t0^1
+DEAL::bp1_1 = 1.00000 * t1^1
+DEAL::bp1_0 * 2 * bp1_1 / 2 = 1.00000 * t0^1 * t1^1
+DEAL::
+DEAL::
+DEAL::2D:
+DEAL::bp2_0 = 2.00000 * t0^1
+DEAL::bp2_1 = 3.00000 * t1^1
+DEAL::bp2_2 = 1.00000 * t2^1
+DEAL::bp2_0 + bp2_1 = 3.00000 * t1^1 + 2.00000 * t0^1
+DEAL::(bp2_0 + bp2_1) * bp2_0 = 6.00000 * t0^1 * t1^1 + 4.00000 * t0^2
+DEAL::bp2_0 * bp2_0 + bp2_1 * bp2_0 = 6.00000 * t0^1 * t1^1 + 4.00000 * t0^2
+DEAL::bp2_1 * bp2_0 + bp2_0 * bp2_0 = 6.00000 * t0^1 * t1^1 + 4.00000 * t0^2
+DEAL::d/dx bp2_0 = -2.00000
+DEAL::d/dy bp2_0 = -2.00000
+DEAL::d/dx bp2_2 = 0.00000
+DEAL::d/dy bp2_2 = 1.00000
+DEAL::
+DEAL::Test with TRI6
+DEAL::p = -1.00000 * t0^1 + 2.00000 * t0^2
+DEAL::p_x = 1.00000 + -4.00000 * t0^1
+DEAL::p_y = 1.00000 + -4.00000 * t0^1
+DEAL::
+DEAL::p = -1.00000 * t1^1 + 2.00000 * t1^2
+DEAL::p_x = -1.00000 + 4.00000 * t1^1
+DEAL::p_y = 0.00000
+DEAL::
+DEAL::p = -1.00000 * t2^1 + 2.00000 * t2^2
+DEAL::p_x = 0.00000
+DEAL::p_y = -1.00000 + 4.00000 * t2^1
+DEAL::
+DEAL::p = 4.00000 * t0^1 * t1^1
+DEAL::p_x = -4.00000 * t1^1 + 4.00000 * t0^1
+DEAL::p_y = -4.00000 * t1^1
+DEAL::
+DEAL::p = 4.00000 * t1^1 * t2^1
+DEAL::p_x = 4.00000 * t2^1
+DEAL::p_y = 4.00000 * t1^1
+DEAL::
+DEAL::p = 4.00000 * t0^1 * t2^1
+DEAL::p_x = -4.00000 * t2^1
+DEAL::p_y = -4.00000 * t2^1 + 4.00000 * t0^1
+DEAL::
+DEAL::
+DEAL::Test with TET4
+DEAL::Test with TET4 - Success
+DEAL::Test with TET10
+DEAL::Test with TET10 - Success

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.