polynomials_piecewise.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 2017 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_piecewise.h>


DEAL_II_NAMESPACE_OPEN



namespace Polynomials
{
  template <typename number>
  PiecewisePolynomial<number>::PiecewisePolynomial(
    const Polynomial<number> &coefficients_on_interval,
    const unsigned int        n_intervals,
    const unsigned int        interval,
    const bool                spans_next_interval)
    : polynomial(coefficients_on_interval)
    , n_intervals(n_intervals)
    , interval(interval)
    , spans_two_intervals(spans_next_interval)
  {
    Assert(n_intervals > 0, ExcMessage("No intervals given"));
    AssertIndexRange(interval, n_intervals);
  }



  template <typename number>
  void
  PiecewisePolynomial<number>::value(const number         x,
                                     std::vector<number> &values) const
  {
    Assert(values.size() > 0, ExcZero());

    value(x, values.size() - 1, values.data());
  }



  template <typename number>
  void
  PiecewisePolynomial<number>::value(const number       x,
                                     const unsigned int n_derivatives,
                                     number *           values) const
  {
    // shift polynomial if necessary
    number y                      = x;
    double derivative_change_sign = 1.;
    if (n_intervals > 0)
      {
        const number step = 1. / n_intervals;
        // polynomial spans over two intervals
        if (spans_two_intervals)
          {
            const double offset = step * interval;
            if (x < offset || x > offset + step + step)
              {
                for (unsigned int k = 0; k <= n_derivatives; ++k)
                  values[k] = 0;
                return;
              }
            else if (x < offset + step)
              y = x - offset;
            else
              {
                derivative_change_sign = -1.;
                y                      = offset + step + step - x;
              }
          }
        else
          {
            const double offset = step * interval;
            if (x < offset || x > offset + step)
              {
                for (unsigned int k = 0; k <= n_derivatives; ++k)
                  values[k] = 0;
                return;
              }
            else
              y = x - offset;
          }

        // on subinterval boundaries, cannot evaluate derivatives properly, so
        // set them to zero
        if ((std::abs(y) < 1e-14 &&
             (interval > 0 || derivative_change_sign == -1.)) ||
            (std::abs(y - step) < 1e-14 &&
             (interval < n_intervals - 1 || derivative_change_sign == -1.)))
          {
            values[0] = value(x);
            for (unsigned int d = 1; d <= n_derivatives; ++d)
              values[d] = 0;
            return;
          }
      }

    polynomial.value(y, n_derivatives, values);

    // change sign if necessary
    for (unsigned int j = 1; j <= n_derivatives; j += 2)
      values[j] *= derivative_change_sign;
  }



  std::vector<PiecewisePolynomial<double>>
  generate_complete_Lagrange_basis_on_subdivisions(
    const unsigned int n_subdivisions,
    const unsigned int base_degree)
  {
    std::vector<Polynomial<double>> p_base =
      LagrangeEquidistant::generate_complete_basis(base_degree);
    for (unsigned int i = 0; i < p_base.size(); ++i)
      p_base[i].scale(n_subdivisions);

    std::vector<PiecewisePolynomial<double>> p;
    p.reserve(n_subdivisions * base_degree + 1);

    p.emplace_back(p_base[0], n_subdivisions, 0, false);
    for (unsigned int s = 0; s < n_subdivisions; ++s)
      for (unsigned int i = 0; i < base_degree; ++i)
        p.emplace_back(p_base[i + 1],
                       n_subdivisions,
                       s,
                       i == (base_degree - 1) && s < n_subdivisions - 1);
    return p;
  }

} // namespace Polynomials

// ------------------ explicit instantiations --------------- //

namespace Polynomials
{
  template class PiecewisePolynomial<float>;
  template class PiecewisePolynomial<double>;
  template class PiecewisePolynomial<long double>;
} // namespace Polynomials

DEAL_II_NAMESPACE_CLOSE