cluster/cebysevK.hh

Go to the documentation of this file.

/* Cebysev kernel expansion for arbitrary kernels
 */

#ifndef clusterCebysevK_hh
#define clusterCebysevK_hh

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

namespace cluster {

  // ******************************************************* CebysevKFTraits **

  template <class K>
  class CebysevKFTraits {
    public:
    typedef typename K::F Fkrnl;

    static Fkrnl eval(K& k, const concepts::Real3d& x,
          const concepts::Real3d& y) {return k(x, y);}
  };

  // ************************************************************* CebysevKF **

  template <class K = bem::Laplace, class Fspc = concepts::Real>
  class CebysevKF : public ExpansionF<Fspc> {
  public:
    typedef CebysevKFTraits<K>     Traits;
    typedef typename Traits::Fkrnl Fkrnl;

  private:
    uint m_;
    concepts::Real eta_;
    concepts::Real** cebysev_;
    concepts::Real* x_;
    concepts::Real* y_;
    concepts::Real* z_;
    Fkrnl* kernel_;
    K& krnl_;
    Fkrnl* d_;
    Fkrnl* dd_;

    inline int index_(int ix, int iy, int iz) const;

  public:
    CebysevKF(K& krnl, uint m, concepts::Real eta);

    virtual ~CebysevKF();

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

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

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

  template <class K, class Fspc>
  inline int CebysevKF<K, 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 K, class Fspc>
  inline void CebysevKF<K, Fspc>::evaluate(uint m, const concepts::Real3d& z,
             FColExp* Fexp) const {
    FColF<Fkrnl>* fexp = dynamic_cast<FColF<Fkrnl>*>(Fexp);
    if (fexp) {
      Fkrnl* F = fexp->value();
      fit(m, z, F);  ceby2poly(m, z, F);
      return;
    }

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

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

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

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

} // namespace cluster

#endif // clusterCebysevK_hh