polynomials_barycentric.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2021 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
 
#ifndef dealii_simplex_barycentric_polynomials_h
#define dealii_simplex_barycentric_polynomials_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exception_macros.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/point.h>
#include <deal.II/base/scalar_polynomials_base.h>
#include <deal.II/base/table.h>
#include <deal.II/base/table_indices.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
 
#include <Kokkos_Macros.hpp>
 
#include <algorithm>
#include <array>
#include <cstddef>
#include <limits>
#include <memory>
#include <ostream>
#include <string>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
template <int dim, typename Number = double>
class BarycentricPolynomial
{
public:
  BarycentricPolynomial();
 
  BarycentricPolynomial(const TableIndices<dim + 1> &powers,
                        const Number                 coefficient);
 
  static BarycentricPolynomial<dim, Number>
  monomial(const unsigned int d);
 
  void
  print(std::ostream &out) const;
 
  TableIndices<dim + 1>
  degrees() const;
 
  BarycentricPolynomial<dim, Number>
  operator-() const;
 
  template <typename Number2>
  BarycentricPolynomial<dim, Number>
  operator+(const Number2 &a) const;
 
  template <typename Number2>
  BarycentricPolynomial<dim, Number>
  operator-(const Number2 &a) const;
 
  template <typename Number2>
  BarycentricPolynomial<dim, Number>
  operator*(const Number2 &a) const;
 
  template <typename Number2>
  BarycentricPolynomial<dim, Number>
  operator/(const Number2 &a) const;
 
  BarycentricPolynomial<dim, Number>
  operator+(const BarycentricPolynomial<dim, Number> &augend) const;
 
  BarycentricPolynomial<dim, Number>
  operator-(const BarycentricPolynomial<dim, Number> &augend) const;
 
  BarycentricPolynomial<dim, Number>
  operator*(const BarycentricPolynomial<dim, Number> &multiplicand) const;
 
  BarycentricPolynomial<dim, Number>
  barycentric_derivative(const unsigned int coordinate) const;
 
  BarycentricPolynomial<dim, Number>
  derivative(const unsigned int coordinate) const;
 
  Number
  value(const Point<dim> &point) const;
 
  std::size_t
  memory_consumption() const;
 
protected:
  Table<dim + 1, Number> coefficients;
 
  static TableIndices<dim + 1>
  index_to_indices(const std::size_t           &index,
                   const TableIndices<dim + 1> &extents);
};
 
template <int dim>
class BarycentricPolynomials : public ScalarPolynomialsBase<dim>
{
public:
  using PolyType = BarycentricPolynomial<dim>;
 
  using GradType = std::array<PolyType, dim>;
 
  using HessianType = std::array<GradType, dim>;
 
  using ThirdDerivativesType = std::array<HessianType, dim>;
 
  using FourthDerivativesType = std::array<ThirdDerivativesType, dim>;
 
  static constexpr unsigned int dimension = dim;
 
  static BarycentricPolynomials<dim>
  get_fe_p_basis(const unsigned int degree);
 
  BarycentricPolynomials(
    const std::vector<BarycentricPolynomial<dim>> &polynomials);
 
  // We need to define the destructor to work around a compiler bug with hipcc
  // 6.4. Using default also triggers the error.
  // NOLINTNEXTLINE(modernize-use-equals-default)
  virtual ~BarycentricPolynomials(){};
 
  const BarycentricPolynomial<dim> &
  operator[](const std::size_t i) const;
 
  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;
 
  double
  compute_value(const unsigned int i, const Point<dim> &p) const override;
 
  Tensor<1, dim>
  compute_1st_derivative(const unsigned int i,
                         const Point<dim>  &p) const override;
 
  Tensor<2, dim>
  compute_2nd_derivative(const unsigned int i,
                         const Point<dim>  &p) const override;
 
  Tensor<3, dim>
  compute_3rd_derivative(const unsigned int i,
                         const Point<dim>  &p) const override;
 
  Tensor<4, dim>
  compute_4th_derivative(const unsigned int i,
                         const Point<dim>  &p) const override;
 
  Tensor<1, dim>
  compute_grad(const unsigned int i, const Point<dim> &p) const override;
 
  Tensor<2, dim>
  compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
 
  virtual std::size_t
  memory_consumption() const override;
 
  std::string
  name() const override;
 
  virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
  clone() const override;
 
protected:
  std::vector<PolyType>              polys;
  std::vector<GradType>              poly_grads;
  std::vector<HessianType>           poly_hessians;
  std::vector<ThirdDerivativesType>  poly_third_derivatives;
  std::vector<FourthDerivativesType> poly_fourth_derivatives;
};
 
// non-member template functions for algebra
 
template <int dim, typename Number1, typename Number2>
BarycentricPolynomial<dim, Number1>
operator*(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
{
  return bp * Number1(a);
}
 
template <int dim, typename Number1, typename Number2>
BarycentricPolynomial<dim, Number1>
operator+(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
{
  return bp + Number1(a);
}
 
template <int dim, typename Number1, typename Number2>
BarycentricPolynomial<dim, Number1>
operator-(const Number2 &a, const BarycentricPolynomial<dim, Number1> &bp)
{
  return bp - Number1(a);
}
 
template <int dim, typename Number>
std::ostream &
operator<<(std::ostream &out, const BarycentricPolynomial<dim, Number> &bp)
{
  bp.print(out);
  return out;
}
 
// Template function definitions
 
// 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 *= Utilities::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> &extents)
{
  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 *= extents[n2];
      result[n] = temp / slice_size;
      temp %= slice_size;
    }
  return result;
}
 
 
 
template <int dim>
const BarycentricPolynomial<dim> &
BarycentricPolynomials<dim>::operator[](const std::size_t i) const
{
  AssertIndexRange(i, polys.size());
  return polys[i];
}
 
DEAL_II_NAMESPACE_CLOSE
 
#endif