aglowav2/element.hh

Go to the documentation of this file.

/* Agglomerated wavelet elements
 */

#ifndef aglowav2Element_hh
#define aglowav2Element_hh

#ifdef __GUNG__
#pragma interface
#endif

#include "basics/vectorsMatricesForward.hh"
#include "space/tmatrix.hh"
#include "space/element.hh"
#include "bem/element.hh"

namespace aglowav2 {

  // ************************************************************* Haar3dXXX **

  template<class F = concepts::Real>
  class Haar3dXXX : public concepts::Element<F> {
  public:
    class Key {
    public:
      inline Key() : l_(0), j_(0) {}
      inline Key(uint l, uint j) : l_(l), j_(j) {}
      inline Key(const Key& k) : l_(k.l()), j_(k.j()) {}

      int operator==(const Key& k) const {return l_ == k.l() && j_ == k.j();}
      int operator<(const Key& k) const {
  return l_ < k.l() || (l_ == k.l() && j_ < k.j());
      }
      inline uint& l() {return l_;}
      inline uint l() const {return l_;}
      inline uint& j() {return j_;}
      inline uint j() const {return j_;}

    private:
      uint l_;
      uint j_;
    };

    inline Haar3dXXX(const concepts::Real3d& cntr, concepts::Real rd,
         concepts::Real sz, const Key& key = Key());
    virtual ~Haar3dXXX(){};

    inline const Key& key() const {return key_;}
    inline Key& key() {return key_;}
    inline const concepts::Real3d& center() const {return cntr_;}
    inline concepts::Real radius() const {return rd_;}
    inline concepts::Real size() const {return sz_;}

    virtual Haar3dXXX<F>* child() const = 0;
    virtual Haar3dXXX<F>* link() const = 0;

  private:
    concepts::Real3d cntr_;
    concepts::Real   rd_;
    Key key_;
    concepts::Real   sz_;
  };

  template<class F>
  Haar3dXXX<F>::Haar3dXXX(const concepts::Real3d& cntr, concepts::Real rd,
        concepts::Real sz, const Key& key)
    : cntr_(cntr), rd_(rd), key_(key), sz_(sz) {
  }

  // **************************************************************** Matrix **

  template<class F = concepts::Real>
  class Matrix {
  public:
    inline Matrix(uint n, uint m, F* val, bool row)
      : n_(n), m_(m), val_(val), row_(row) {}

    inline const F* source() const {return val_;}
    inline F* destination() const {return val_;}
    inline uint nrow() const {return n_;}
    inline uint ncol() const {return m_;}
    inline bool roworder() const {return row_;}
    inline void transpose() {
      row_ = row_ ? 0 : 1;  uint n = n_;  n_ = m_;  m_ = n;
    }

  private:
    uint n_;
    uint m_;
    F*   val_;
    bool row_;
  };


  // ****************************************************************** M000 **

  class M000 {
  public:
    friend std::ostream& operator<<(std::ostream& os, const M000& m);

    M000(const concepts::Real* m, uint n);
    ~M000() {if (m_) delete[] m_;}

    template<class F>
    void mult(const Matrix<F>& src, Matrix<F>& dst,
        uint n = 0, uint m = 0) const;

    template<class F>
    void mult_T(const Matrix<F>& src, Matrix<F>& dst,
    uint n = 0, uint m = 0) const;

    template<class F>
    void mult_T(const F* src, F* dst, const uint s, const uint t) const;
    inline const concepts::Real& operator()(uint i, uint j) const {
      return m_[i*n_ + j];
    }
    inline uint n() const {return n_;}

  private:
    concepts::Real* m_;
    uint n_;
  };

  template<class F>
  void M000::mult(const Matrix<F>& src, Matrix<F>& dst, uint n, uint m) const {
    const F* sv = src.source();
    F* dv = dst.destination();
    const concepts::Real* mat = &(*this)(n, m);

    if (src.roworder()) {
      if (dst.roworder()) {

  for(uint i = 0; i < dst.nrow(); i++) {
    const concepts::Real* matrow = mat + i*n_;
    for(uint j = 0; j < dst.ncol(); j++) {
      *dv = 0.0;
      const F* svj = sv + j;
      for(uint k = 0; k < src.nrow(); k++) {
        *dv += *svj * matrow[k];  svj += src.ncol();
      }
      dv++;
    }
  } // for i

      } // if dst row order
      else {

  for(uint i = 0; i < dst.ncol(); i++) {
    for(uint j = 0; j < dst.nrow(); j++) {
      *dv = 0.0;
      const concepts::Real* matrow = mat + j*n_;
      const F* svi = sv + i;
      for(uint k = 0; k < src.nrow(); k++) {
        *dv += *svi * matrow[k];  svi += src.ncol();
      }
      dv++;
    }
  } // for i

      } // if dst col order
    } // if src row order
    else {
      if (dst.roworder()) {

  for(uint i = 0; i < dst.nrow(); i++) {
    const concepts::Real* matrow = mat + i*n_;
    for(uint j = 0; j < dst.ncol(); j++) {
      *dv = 0.0;
      const F* svj = sv + j*src.nrow();
      for(uint k = 0; k < src.nrow(); k++) {
        *dv += *svj++ * matrow[k];
      }
      dv++;
    }
  } // for i

      } // if dst row order
      else {

  for(uint i = 0; i < dst.ncol(); i++) {
    for(uint j = 0; j < dst.nrow(); j++) {
      *dv = 0.0;
      const concepts::Real* matrow = mat + j*n_;
      const F* svi = sv + i*src.nrow();
      for(uint k = 0; k < src.nrow(); k++) {
        *dv += *svi++ * matrow[k];
      }
      dv++;
    }
  } // for i

      } // if dst col order
    } // if src col order
  }

  template<class F>
  void M000::mult_T(const Matrix<F>& src, Matrix<F>& dst,
        uint n, uint m) const {
    const F* sv = src.source();
    F* dv = dst.destination();
    const concepts::Real* mat = &(*this)(m, n);

    if (src.roworder()) {
      if (dst.roworder()) {

  for(uint i = 0; i < dst.nrow(); i++) {
    for(uint j = 0; j < dst.ncol(); j++) {
      *dv = 0.0;
      const concepts::Real* matrow = mat + i;
      const F* svj = sv + j;
      for(uint k = 0; k < src.nrow(); k++) {
        *dv += *svj * *matrow;  svj += src.ncol();  matrow += n_;
      }
      dv++;
    }
  } // for i

      } // if dst row order
      else {

  for(uint i = 0; i < dst.ncol(); i++) {
    for(uint j = 0; j < dst.nrow(); j++) {
      *dv = 0.0;
      const concepts::Real* matrow = mat + j;
      const F* svi = sv + i;
      for(uint k = 0; k < src.nrow(); k++) {
        *dv += *svi * *matrow;  svi += src.ncol();  matrow += n_;
      }
      dv++;
    }
  } // for i

      } // if dst col order
    } // if src row order
    else {
      if (dst.roworder()) {

  for(uint i = 0; i < dst.nrow(); i++) {
    for(uint j = 0; j < dst.ncol(); j++) {
      *dv = 0.0;
      const concepts::Real* matrow = mat + i;
      const F* svj = sv + j*src.nrow();
      for(uint k = 0; k < src.nrow(); k++) {
        *dv += *svj++ * *matrow;  matrow += n_;
      }
      dv++;
    }
  } // for i

      } // if dst row order
      else {

  for(uint i = 0; i < dst.ncol(); i++) {
    for(uint j = 0; j < dst.nrow(); j++) {
      *dv = 0.0;
      const concepts::Real* matrow = mat + j;
      const F* svi = sv + i*src.nrow();
      for(uint k = 0; k < src.nrow(); k++) {
        *dv += *svi++ * *matrow;  matrow += n_;
      }
      dv++;
    }
  } // for i

      } // if dst col order
    } // if src col order
  }

  template<class F>
  void M000::mult_T(const F* src, F* dst, const uint s, const uint t) const {
    for(uint i = 0; i < t; i++) {
      for(uint j = 0; j < s; j++) {
  concepts::Real* mi = m_ + i;
  const F* sj = src + j;
  *dst = 0.0;
  for(uint k = 0; k < n_; k++, mi += n_, sj += s)  *dst += *sj * *mi;
  dst++;
      }
    }
  }

  // ************************************************************* Haar3d000 **

  template<class F = concepts::Real>
  class Haar3d000 : public Haar3dXXX<F> {
  public:
    Haar3d000(const concepts::Real3d& cntr, concepts::Real rd,
        const concepts::Real* m, uint n, uint idx[], uint idxn,
        concepts::Real sz, Haar3d000<F>* chld,
        const bem::Constant3d002<F>** elm, uint nelm);
    ~Haar3d000() {delete T_;  if (nelm_)  delete[] elm_;}

    Haar3d000<F>* child() const {return chld_;}
    Haar3d000<F>* link() const {return lnk_;}
    inline Haar3d000<F>*& link() {return lnk_;}
    const bem::Constant3d002<F>* element(uint j) const;
    inline uint nelement() const {return nelm_;}
    inline uint gamma() const {return gamma_;}
    inline uint& gamma() {return gamma_;}

    inline const concepts::TMatrixBase<F>& T() const {return *T_;}
    inline void replaceT(uint idx[], uint idxn);
    inline const M000& trafoM() const {return M_;}
    inline uint index() const {return idx_;}

  protected:
    std::ostream& info(std::ostream& os) const;

  private:
    M000 M_;
    const concepts::TIndex<F>* T_;
    concepts::Real sz_;
    Haar3d000* lnk_;
    Haar3d000* chld_;
    const bem::Constant3d002<F>** elm_;
    uint nelm_;
    uint gamma_;
    uint idx_;

    static uint index_;
  };

  template<class F>
  Haar3d000<F>::Haar3d000(const concepts::Real3d& cntr, concepts::Real rd,
        const concepts::Real* m, uint n,
        uint idx[], uint idxn, concepts::Real sz,
        Haar3d000<F>* chld,
        const bem::Constant3d002<F>** elm, uint nelm)
    : Haar3dXXX<F>(cntr, rd, sz), M_(m, n),
      lnk_(0), chld_(chld), nelm_(nelm), gamma_(1), idx_(index_++) {

    T_ = new concepts::TIndex<F>(idxn, idxn, idx);

    if (nelm) {
      elm_ = new const bem::Constant3d002<F>*[nelm_];
      std::memcpy(elm_, elm, nelm_*sizeof(elm[0]));
    }
    else elm_ = 0;
  }

  template<class F>
  const bem::Constant3d002<F>* Haar3d000<F>::element(uint j) const {
    return (j < nelm_) ? elm_[j] : (bem::Constant3d002<F>*)0;
  }

  template<class F>
  void Haar3d000<F>::replaceT(uint idx[], uint idxn) {
    delete T_;
    T_ = new concepts::TIndex<F>(idxn, idxn, idx);
  }

} // namespace aglowav2

#endif //aglowav2Element_hh