polynomial_space.cc

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2002 - 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.
//
// -----------------------------------------------------------------------------
 
#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/polynomial.h>
#include <deal.II/base/polynomial_space.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/types.h>
#include <deal.II/base/utilities.h>
 
#include <boost/container/small_vector.hpp>
 
#include <Kokkos_Macros.hpp>
 
#include <array>
#include <memory>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
 
template <int dim>
unsigned int
PolynomialSpace<dim>::n_polynomials(const unsigned int n)
{
  unsigned int n_pols = n;
  for (unsigned int i = 1; i < dim; ++i)
    {
      n_pols *= (n + i);
      n_pols /= (i + 1);
    }
  return n_pols;
}
 
 
template <>
unsigned int
PolynomialSpace<0>::n_polynomials(const unsigned int)
{
  return 0;
}
 
 
template <>
std::array<unsigned int, 1>
PolynomialSpace<1>::compute_index(const unsigned int i) const
{
  AssertIndexRange(i, index_map.size());
  return {{index_map[i]}};
}
 
 
 
template <>
std::array<unsigned int, 2>
PolynomialSpace<2>::compute_index(const unsigned int i) const
{
  AssertIndexRange(i, index_map.size());
  const unsigned int n = index_map[i];
  // there should be a better way to
  // write this function (not
  // linear in n_1d), someone
  // should think about this...
  const unsigned int n_1d = polynomials.size();
  unsigned int       k    = 0;
  for (unsigned int iy = 0; iy < n_1d; ++iy)
    if (n < k + n_1d - iy)
      {
        return {{n - k, iy}};
      }
    else
      k += n_1d - iy;
 
  DEAL_II_ASSERT_UNREACHABLE();
  return {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}};
}
 
 
 
template <>
std::array<unsigned int, 3>
PolynomialSpace<3>::compute_index(const unsigned int i) const
{
  AssertIndexRange(i, index_map.size());
  const unsigned int n = index_map[i];
  // there should be a better way to
  // write this function (not
  // quadratic in n_1d), someone
  // should think about this...
  //
  // (ah, and yes: the original
  // algorithm was even cubic!)
  const unsigned int n_1d = polynomials.size();
  unsigned int       k    = 0;
  for (unsigned int iz = 0; iz < n_1d; ++iz)
    for (unsigned int iy = 0; iy < n_1d - iz; ++iy)
      if (n < k + n_1d - iy - iz)
        {
          return {{n - k, iy, iz}};
        }
      else
        k += n_1d - iy - iz;
 
  DEAL_II_ASSERT_UNREACHABLE();
  return {{numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}};
}
 
 
template <int dim>
void
PolynomialSpace<dim>::set_numbering(const std::vector<unsigned int> &renumber)
{
  Assert(renumber.size() == index_map.size(),
         ExcDimensionMismatch(renumber.size(), index_map.size()));
 
  index_map = renumber;
  for (unsigned int i = 0; i < index_map.size(); ++i)
    index_map_inverse[index_map[i]] = i;
}
 
 
 
template <int dim>
double
PolynomialSpace<dim>::compute_value(const unsigned int i,
                                    const Point<dim>  &p) const
{
  const auto ix = compute_index(i);
  // take the product of the
  // polynomials in the various space
  // directions
  double result = 1.;
  for (unsigned int d = 0; d < dim; ++d)
    result *= polynomials[ix[d]].value(p[d]);
  return result;
}
 
 
 
template <int dim>
Tensor<1, dim>
PolynomialSpace<dim>::compute_grad(const unsigned int i,
                                   const Point<dim>  &p) const
{
  const auto ix = compute_index(i);
 
  Tensor<1, dim> result;
  for (unsigned int d = 0; d < dim; ++d)
    result[d] = 1.;
 
  // get value and first derivative
  std::vector<double> v(2);
  for (unsigned int d = 0; d < dim; ++d)
    {
      polynomials[ix[d]].value(p[d], v);
      result[d] *= v[1];
      for (unsigned int d1 = 0; d1 < dim; ++d1)
        if (d1 != d)
          result[d1] *= v[0];
    }
  return result;
}
 
 
template <int dim>
Tensor<2, dim>
PolynomialSpace<dim>::compute_grad_grad(const unsigned int i,
                                        const Point<dim>  &p) const
{
  const auto ix = compute_index(i);
 
  Tensor<2, dim> result;
  for (unsigned int d = 0; d < dim; ++d)
    for (unsigned int d1 = 0; d1 < dim; ++d1)
      result[d][d1] = 1.;
 
  // get value, first and second
  // derivatives
  std::vector<double> v(3);
  for (unsigned int d = 0; d < dim; ++d)
    {
      polynomials[ix[d]].value(p[d], v);
      result[d][d] *= v[2];
      for (unsigned int d1 = 0; d1 < dim; ++d1)
        {
          if (d1 != d)
            {
              result[d][d1] *= v[1];
              result[d1][d] *= v[1];
              for (unsigned int d2 = 0; d2 < dim; ++d2)
                if (d2 != d)
                  result[d1][d2] *= v[0];
            }
        }
    }
  return result;
}
 
 
template <int dim>
void
PolynomialSpace<dim>::evaluate(
  const Point<dim>            &p,
  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
{
  const unsigned int n_1d = polynomials.size();
 
  Assert(values.size() == this->n() || values.empty(),
         ExcDimensionMismatch2(values.size(), this->n(), 0));
  Assert(grads.size() == this->n() || grads.empty(),
         ExcDimensionMismatch2(grads.size(), this->n(), 0));
  Assert(grad_grads.size() == this->n() || grad_grads.empty(),
         ExcDimensionMismatch2(grad_grads.size(), this->n(), 0));
  Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
         ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0));
  Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
         ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0));
 
  unsigned int v_size = 0;
  bool update_values = false, update_grads = false, update_grad_grads = false;
  bool update_3rd_derivatives = false, update_4th_derivatives = false;
  if (values.size() == this->n())
    {
      update_values = true;
      v_size        = 1;
    }
  if (grads.size() == this->n())
    {
      update_grads = true;
      v_size       = 2;
    }
  if (grad_grads.size() == this->n())
    {
      update_grad_grads = true;
      v_size            = 3;
    }
  if (third_derivatives.size() == this->n())
    {
      update_3rd_derivatives = true;
      v_size                 = 4;
    }
  if (fourth_derivatives.size() == this->n())
    {
      update_4th_derivatives = true;
      v_size                 = 5;
    }
 
  // Store data in a single
  // object. Access is by
  // v[d][n][o]
  //  d: coordinate direction
  //  n: number of 1d polynomial
  //  o: order of derivative
  //
  // Our rule-of-thumb for stack arrays is 200 elements in a small_vector. Here,
  // we have (in 3d) 3 * 14 * 5 = 210 elements, in which 5 is the maximum number
  // derivatives (4) plus one value.
  std::array<boost::container::small_vector<std::array<double, 5>, 14>, dim> v;
  for (unsigned int d = 0; d < dim; ++d)
    {
      v[d].resize(n_1d);
      for (unsigned int i = 0; i < n_1d; ++i)
        {
          Assert(v_size > 0, ExcInternalError());
          Assert(v_size <= v[d][i].size(), ExcInternalError());
          if constexpr (running_in_debug_mode())
            v[d][i].fill(std::numeric_limits<double>::signaling_NaN());
          polynomials[i].value(p[d], v_size - 1, v[d][i].data());
        }
    }
 
  if (update_values)
    {
      unsigned int k = 0;
 
      for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
        for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
          for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
            values[index_map_inverse[k++]] = v[0][ix][0] *
                                             ((dim > 1) ? v[1][iy][0] : 1.) *
                                             ((dim > 2) ? v[2][iz][0] : 1.);
    }
 
  if (update_grads)
    {
      unsigned int k = 0;
 
      for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
        for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
          for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
            {
              const unsigned int k2 = index_map_inverse[k++];
              for (unsigned int d = 0; d < dim; ++d)
                grads[k2][d] = v[0][ix][(d == 0) ? 1 : 0] *
                               ((dim > 1) ? v[1][iy][(d == 1) ? 1 : 0] : 1.) *
                               ((dim > 2) ? v[2][iz][(d == 2) ? 1 : 0] : 1.);
            }
    }
 
  if (update_grad_grads)
    {
      unsigned int k = 0;
 
      for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
        for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
          for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
            {
              const unsigned int k2 = index_map_inverse[k++];
              for (unsigned int d1 = 0; d1 < dim; ++d1)
                for (unsigned int d2 = 0; d2 < dim; ++d2)
                  {
                    // Derivative
                    // order for each
                    // direction
                    const unsigned int j0 =
                      ((d1 == 0) ? 1 : 0) + ((d2 == 0) ? 1 : 0);
                    const unsigned int j1 =
                      ((d1 == 1) ? 1 : 0) + ((d2 == 1) ? 1 : 0);
                    const unsigned int j2 =
                      ((d1 == 2) ? 1 : 0) + ((d2 == 2) ? 1 : 0);
 
                    grad_grads[k2][d1][d2] = v[0][ix][j0] *
                                             ((dim > 1) ? v[1][iy][j1] : 1.) *
                                             ((dim > 2) ? v[2][iz][j2] : 1.);
                  }
            }
    }
 
  if (update_3rd_derivatives)
    {
      unsigned int k = 0;
 
      for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
        for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
          for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
            {
              const unsigned int k2 = index_map_inverse[k++];
              for (unsigned int d1 = 0; d1 < dim; ++d1)
                for (unsigned int d2 = 0; d2 < dim; ++d2)
                  for (unsigned int d3 = 0; d3 < dim; ++d3)
                    {
                      // Derivative
                      // order for each
                      // direction
                      std::array<unsigned int, dim> deriv_order{};
                      for (unsigned int x = 0; x < dim; ++x)
                        {
                          if (d1 == x)
                            ++deriv_order[x];
                          if (d2 == x)
                            ++deriv_order[x];
                          if (d3 == x)
                            ++deriv_order[x];
                        }
 
                      third_derivatives[k2][d1][d2][d3] =
                        v[0][ix][deriv_order[0]] *
                        ((dim > 1) ? v[1][iy][deriv_order[1]] : 1.) *
                        ((dim > 2) ? v[2][iz][deriv_order[2]] : 1.);
                    }
            }
    }
 
  if (update_4th_derivatives)
    {
      unsigned int k = 0;
 
      for (unsigned int iz = 0; iz < ((dim > 2) ? n_1d : 1); ++iz)
        for (unsigned int iy = 0; iy < ((dim > 1) ? n_1d - iz : 1); ++iy)
          for (unsigned int ix = 0; ix < n_1d - iy - iz; ++ix)
            {
              const unsigned int k2 = index_map_inverse[k++];
              for (unsigned int d1 = 0; d1 < dim; ++d1)
                for (unsigned int d2 = 0; d2 < dim; ++d2)
                  for (unsigned int d3 = 0; d3 < dim; ++d3)
                    for (unsigned int d4 = 0; d4 < dim; ++d4)
                      {
                        // Derivative
                        // order for each
                        // direction
                        std::array<unsigned int, dim> deriv_order{};
                        for (unsigned int x = 0; x < dim; ++x)
                          {
                            if (d1 == x)
                              ++deriv_order[x];
                            if (d2 == x)
                              ++deriv_order[x];
                            if (d3 == x)
                              ++deriv_order[x];
                            if (d4 == x)
                              ++deriv_order[x];
                          }
 
                        fourth_derivatives[k2][d1][d2][d3][d4] =
                          v[0][ix][deriv_order[0]] *
                          ((dim > 1) ? v[1][iy][deriv_order[1]] : 1.) *
                          ((dim > 2) ? v[2][iz][deriv_order[2]] : 1.);
                      }
            }
    }
}
 
 
 
template <int dim>
std::unique_ptr<ScalarPolynomialsBase<dim>>
PolynomialSpace<dim>::clone() const
{
  return std::make_unique<PolynomialSpace<dim>>(*this);
}
 
 
template class PolynomialSpace<1>;
template class PolynomialSpace<2>;
template class PolynomialSpace<3>;
 
DEAL_II_NAMESPACE_CLOSE