cluster/cebysev.hh

Go to the documentation of this file.

/* Cebysev kernel expansion
 */

#ifndef cebysev_hh
#define cebysev_hh

#include "basics/exceptions.hh"
#include "cluster/expansion.hh"

namespace cluster {

  // ******************************************************* CebysevLaplaceF **

  template <class Fspc = concepts::Real>
  class CebysevLaplaceF : public ExpansionF<Fspc> {

    class FColCebyLpl : public FColExp {
    public:
      concepts::Real* val;

      inline FColCebyLpl(uint blksz) {val = new concepts::Real[blksz];}
      inline ~FColCebyLpl() {delete[] val;}

      uint memory(uint blksz) const {
        return sizeof(FColCebyLpl) + blksz * sizeof(concepts::Real);
      }
      void info(std::ostream& os, uint blksz) const {
        for(uint i = 0; i < blksz-1; i++)  os << val[i] << ", ";
        os << val[blksz-1];
      }
    };

    uint m_;
    concepts::Real eta_;
    concepts::Real** cebysev_;
    concepts::Real* x_;
    concepts::Real* y_;
    concepts::Real* z_;
    concepts::Real* kernel_;
    concepts::Real* d_;
    concepts::Real* dd_;
  
    inline int index_(int ix, int iy, int iz) const;
  
  public:
    CebysevLaplaceF(uint m, concepts::Real eta);
    
    virtual ~CebysevLaplaceF();

    uint blksz(uint m) const {return (m * (m * (m + 3) + 2)) / 6;}
    uint m() const {return m_;}
    virtual FColReal* getCol(uint blksz) const {
      return new FColReal(blksz);
    }

    void fit(uint m, const concepts::Real3d& z, concepts::Real F[]) const;
    void ceby2poly(uint m, const concepts::Real3d& z,
       concepts::Real F[]) const;

    inline void evaluate(uint m, const concepts::Real3d& z,
       FColExp* Fexp) const;
    inline void evaluate(uint m, const concepts::Real3d& z,
       FColReal Fexp[]) const;
    inline void apply(uint m, const FColExp* Fexp, const Fspc src[],
          Fspc dst[]) const;
    void apply(uint m, const FColReal Fexp[], const Fspc src[],
         Fspc dst[]) const;
  };

  template <class Fspc>
  inline int CebysevLaplaceF<Fspc>::index_(int ix, int iy, int iz) const {
    int m = ix + iy + iz;
    return (((m+3)*m+2)*m + ((3-ix)*3*ix)) / 6 + m*ix + iy;
  }

  template <class Fspc>
  inline void CebysevLaplaceF<Fspc>::evaluate(uint m,
                const concepts::Real3d& z,
                FColExp* Fexp) const {
    FColReal* fexp = dynamic_cast<FColReal*>(Fexp);
    if (fexp) {
      concepts::Real* F = fexp->value();
      fit(m, z, F); ceby2poly(m, z, F);
      return;
    }

    throw
      conceptsException(concepts::MissingFeature("FColExp not supported"));
  }

  template <class Fspc>
  inline void CebysevLaplaceF<Fspc>::evaluate(uint m,
                const concepts::Real3d& z,
                FColReal Fexp[]) const {
    concepts::Real* F = Fexp->value();
    fit(m, z, F);  ceby2poly(m, z, F);
    return;
  }

  template <class Fspc>
  inline void CebysevLaplaceF<Fspc>::apply(uint m, const FColExp* Fexp,
             const Fspc src[], Fspc dst[]) const {
    const FColReal* fexp = dynamic_cast<const FColReal*>(Fexp);
    if (fexp) {apply(m, fexp, src, dst);  return;}

    throw
      conceptsException(concepts::MissingFeature("FColExp not supported"));
  }

} // namespace cluster

#endif // cebysev_hh