aglowav2/xfy.hh

Go to the documentation of this file.

/* Far field matrices X, F and Y.
 */

#ifndef aglowav2xy_hh
#define aglowav2xy_hh

#include "cluster/expansion.hh"
#include "aglowav2/space.hh"

namespace aglowav2 {

  // ****************************************************************** XY00 **

  template<class F = concepts::Real>
  class XY00 {
  public:
    XY00(const Haar3d0<F>& spc, const cluster::ExpansionXY<F>& exp,
   const cluster::ExpansionXY<F>* shftExp = 0);
    virtual ~XY00();

    inline const Haar3d0<F>& space() const {return spc_;}
    cluster::XYColExpPtr* operator[](uint i) const {return elmptr_[i];}
    inline uint memory() const;

  private:
    void constructor_(cluster::XYColReal* xy,
          const cluster::ExpansionXY<F>* shftExp);

  protected:
    const cluster::ExpansionXY<F>& exp_;

  private:
    const Haar3d0<F>& spc_;
    cluster::XYColExpPtr** elmptr_;
    cluster::XYColExp* xy_;
  };

  template<class F>
  XY00<F>::~XY00() {
    for(uint i = 0; i < spc_.nelm(); i++)  delete elmptr_[i];
    delete[] elmptr_;
    delete xy_;
  }

  template<class F>
  uint XY00<F>::memory() const {
    return sizeof(XY00<F>) + spc_.nelm()*sizeof(elmptr_[0]) + xy_->memory();
  }

  // ******************************************************************* Y00 **

  template<class F = concepts::Real>
  class Y00 : private XY00<F> {
  public:
    inline Y00(const cluster::ExpansionXY<F>& exp, const Haar3d0<F>& spc);

    inline const Haar3d0<F>& space() const {return XY00<F>::space();}
    inline cluster::XYColExpPtr* operator[](uint i) const;
    inline uint memory() const {return XY00<F>::memory();}
  };

  template<class F>
  Y00<F>::Y00(const cluster::ExpansionXY<F>& exp, const Haar3d0<F>& spc)
    : XY00<F>(spc, exp, 0) {}

  template<class F>
  cluster::XYColExpPtr* Y00<F>::operator[](uint i) const {
    return XY00<F>::operator[](i);
  }

  // ******************************************************************* X00 **

  template<class F = concepts::Real>
  class X00 : protected XY00<F> {
  public:
    inline X00(const cluster::ExpansionXY<F>& exp,
         const cluster::ExpansionXY<F>& shftExp, const Haar3d0<F>& spc);

    inline const Haar3d0<F>& space() const {return XY00<F>::space();}
    inline cluster::XYColExpPtr* operator[](uint i) const;
    inline uint memory() const {return XY00<F>::memory();}
    inline void apply(cluster::XYColExpPtr* x, const F* src, F* dst) const;
  };

  template<class F>
  X00<F>::X00(const cluster::ExpansionXY<F>& exp,
        const cluster::ExpansionXY<F>& shftExp, const Haar3d0<F>& spc)
    : XY00<F>(spc, exp, &shftExp) {}

  template<class F>
  cluster::XYColExpPtr* X00<F>::operator[](uint i) const {
    return XY00<F>::operator[](i);
  }

  template<class F>
  void X00<F>::apply(cluster::XYColExpPtr* x, const F* src, F* dst) const {
    this->exp_.apply(x, src, dst);
  }

  // ******************************************************************* F00 **

  template<class F = concepts::Real>
  class F00 {
  public:
    F00(const X00<F>& X, const Y00<F>& Y, const cluster::ExpansionF<F>& exp,
  concepts::Real eta);
    ~F00() {delete foo_;  delete[] fee_;}

    inline const Haar3d0<F>& space() const {return X_.space();}
    inline bool admissible(const Haar3dXXX<F>& elmX,
         const Haar3dXXX<F>& elmY);
    void operator()(uint idxX, uint idxY, F* valXY,
        uint dx, uint dy, uint nx, uint ny);
    void operator()(uint idxX, uint idxY, F* valXY, F* valYX, uint dx, uint dy,
        uint nx, uint ny);

  private:
    const X00<F>& X_;
    const Y00<F>& Y_;
    const cluster::ExpansionF<F>& exp_;
    uint m_;
    concepts::Real eta_;

    F* fee_;
    cluster::FColExp* foo_;
    concepts::Real3d z_;
  };

  template<class F>
  F00<F>::F00(const X00<F>& X, const Y00<F>& Y,
        const cluster::ExpansionF<F>& exp, concepts::Real eta)
    : X_(X), Y_(Y), exp_(exp), m_(exp.m()), eta_(eta) {

    if (&X_.space() != &Y_.space())
      throw conceptsException(concepts::
                  MissingFeature("Different spaces not supported"));

    fee_ = new F[exp.blksz(m_)];
    foo_ = exp.getCol(exp.blksz(m_));
  }

  template<class F>
  bool F00<F>::admissible(const Haar3dXXX<F>& elmX,
        const Haar3dXXX<F>& elmY) {
    z_  = elmY.center();
    z_ -= elmX.center();
    return (elmX.radius() + elmY.radius() < eta_*z_.l2());
  }

} // namespace aglowav2

#endif // aglowav2xy_hh