polynomials_piecewise.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) 2013 - 2024 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_polynomials_piecewise_h
#define dealii_polynomials_piecewise_h
 
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/point.h>
#include <deal.II/base/polynomial.h>
 
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
namespace Polynomials
{
  template <typename number>
  class PiecewisePolynomial : public EnableObserverPointer
  {
  public:
    PiecewisePolynomial(const Polynomial<number> &coefficients_on_interval,
                        const unsigned int        n_intervals,
                        const unsigned int        interval,
                        const bool                spans_next_interval);
 
    PiecewisePolynomial(const std::vector<Point<1, number>> &points,
                        const unsigned int                   index);
 
    number
    value(const number x) const;
 
    void
    value(const number x, std::vector<number> &values) const;
 
    void
    value(const number       x,
          const unsigned int n_derivatives,
          number            *values) const;
 
    unsigned int
    degree() const;
 
    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int version);
 
    virtual std::size_t
    memory_consumption() const;
 
  protected:
    Polynomial<number> polynomial;
 
    unsigned int n_intervals;
 
    unsigned int interval;
 
    bool spans_two_intervals;
 
    std::vector<number> points;
 
    std::vector<number> one_over_lengths;
 
    unsigned int index;
  };
 
 
 
  std::vector<PiecewisePolynomial<double>>
  generate_complete_Lagrange_basis_on_subdivisions(
    const unsigned int n_subdivisions,
    const unsigned int base_degree);
 
  std::vector<PiecewisePolynomial<double>>
  generate_complete_linear_basis_on_subdivisions(
    const std::vector<Point<1>> &points);
 
} // namespace Polynomials
 
 
/* -------------------------- inline functions --------------------- */
 
namespace Polynomials
{
  template <typename number>
  inline unsigned int
  PiecewisePolynomial<number>::degree() const
  {
    if (points.size() > 0)
      return 1;
    return polynomial.degree();
  }
 
 
 
  template <typename number>
  inline number
  PiecewisePolynomial<number>::value(const number x) const
  {
    if (points.size() > 0)
      {
        if (x > points[index])
          return std::max<number>(0.0,
                                  1.0 - (x - points[index]) *
                                          one_over_lengths[index]);
        else if (x < points[index])
          return std::max<number>(0.0,
                                  0.0 + (x - points[index - 1]) *
                                          one_over_lengths[index - 1]);
        else
          return 1.0;
      }
 
    AssertIndexRange(interval, n_intervals);
    number y = x;
    // shift polynomial if necessary
    if (n_intervals > 1)
      {
        const number step = 1. / n_intervals;
 
        // polynomial spans over two intervals
        if (spans_two_intervals == true)
          {
            const number offset = step * interval;
            if (x < offset)
              return 0;
            else if (x > offset + step + step)
              return 0;
            else if (x < offset + step)
              y = x - offset;
            else
              y = offset + step + step - x;
          }
        else
          {
            const number offset = step * interval;
            if (x < offset || x > offset + step)
              return 0;
            else
              y = x - offset;
          }
 
        return polynomial.value(y);
      }
    else
      return polynomial.value(x);
  }
 
 
 
  template <typename number>
  template <class Archive>
  inline void
  PiecewisePolynomial<number>::serialize(Archive &ar, const unsigned int)
  {
    // forward to serialization function in the base class.
    ar &static_cast<EnableObserverPointer &>(*this);
    ar &polynomial;
    ar &n_intervals;
    ar &interval;
    ar &spans_two_intervals;
    ar &points;
    ar &one_over_lengths;
    ar &index;
  }
} // namespace Polynomials
 
DEAL_II_NAMESPACE_CLOSE
 
#endif