numerics/splines.h

// -*- c++ -*-

// +--------------------------------------------------------------------+
// | This file is part of MathTL - the Mathematical Template Library    |
// |                                                                    |
// | Copyright (c) 2002-2009                                            |
// | Thorsten Raasch, Manuel Werner                                     |
// +--------------------------------------------------------------------+

#ifndef _MATHTL_SPLINES_H
#define _MATHTL_SPLINES_H

#include <cmath>
#include <iostream>
#include <utils/array1d.h>
#include <utils/function.h>
#include <algebra/matrix.h>

namespace MathTL
{
  class KnotSequence
  {
  public:
    virtual ~KnotSequence() = 0;
    
    /*
      index of the first knot
     */
    virtual int k0() const = 0;
    
    virtual double knot(const int k) const = 0;
  };

  KnotSequence::~KnotSequence() {}

  template <int d>
  double evaluate_Bspline(const KnotSequence* knots, const int j, const double x)
  {
    double r(0);

    // take care of multiple knots
    double diff = knots->knot(j+d-1) - knots->knot(j);
    if (diff > 0) r += (x - knots->knot(j)) * evaluate_Bspline<d-1>(knots, j, x) / diff;
    diff = knots->knot(j+d) - knots->knot(j+1);
    if (diff > 0) r += (knots->knot(j+d) - x) * evaluate_Bspline<d-1>(knots, j+1, x) / diff;
    
    return r;
  }
  
  template <>
  inline
  double evaluate_Bspline<1>(const KnotSequence* knots, const int j, const double x)
  {
    return (x >= knots->knot(j) && x < knots->knot(j+1) ? 1.0 : 0.0);
  }

  template <int d>
  void
  compute_Bspline_refinement_matrix(const KnotSequence* knots, Matrix<double>& M);


  template <int d>
  double evaluate_Bspline(const Array1D<double>& knots, const unsigned int j, const double x)
  {
    double r(0);

    // take care of multiple knots
    double diff = knots[j+d-1] - knots[j];
    if (diff > 0) r += (x - knots[j]) * evaluate_Bspline<d-1>(knots, j, x) / diff;
    diff = knots[j+d] - knots[j+1];
    if (diff > 0) r += (knots[j+d] - x) * evaluate_Bspline<d-1>(knots, j+1, x) / diff;
    
    return r;
  }

  template <>
  inline
  double evaluate_Bspline<1>(const Array1D<double>& knots, const unsigned int j, const double x)
  {
    return (x >= knots[j] && x < knots[j+1] ? 1.0 : 0.0);
  }

  template <unsigned int d>
  class Spline : public Function<1>
  {
  public:
    Spline()
      : Function<1>(1), knots_(d+1), coeffs_(1)
    {
      for (unsigned int i(0); i <= d; i++)
        knots_[i] = i;

      coeffs_[0] = 1.0;
    }

    Spline(const Array1D<double>& knots, const Array1D<double>& coeffs)
      : Function<1>(1), knots_(knots), coeffs_(coeffs)
    {
      assert(knots.size() == coeffs.size()+d);
    }

    virtual ~Spline() {}

    inline double value(const Point<1>& p,
                        const unsigned int component = 0) const
    {
      double r(0);
      
      for (unsigned int i(0); i < coeffs_.size(); i++) {
        if (coeffs_[i] != 0)
          r += coeffs_[i] * evaluate_Bspline<d>(knots_, i, p(0));
      }

      return r;
    }
  
    void vector_value(const Point<1> &p,
                      Vector<double>& values) const
    {
      values.resize(1, false);
      values[0] = value(p);
    }
    
  protected:
    Array1D<double> knots_;

    Array1D<double> coeffs_;
  };
}

#include <numerics/splines.cpp>

#endif