cebysev.hh

Go to the documentation of this file.

 
#ifndef cebysev_hh
#define cebysev_hh
 
#include "basics/output.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