aglowav/element.hh

Go to the documentation of this file.

/* Agglomerated wavelet elements
 */

#ifndef aglowavElement_hh
#define aglowavElement_hh

#ifdef __GUNG__
  #pragma interface
#endif

#include "basics/outputOperator.hh"
#include "basics/vectorsMatricesForward.hh"
#include "geometry/cell2D.hh"
#include "space/tmatrix.hh"
#include "space/element.hh"
#include "aglowav/tree.hh"

namespace aglowav {

  // ********************************************************************* M **

  template <uint d = 2>
  class M : public concepts::OutputOperator {
  public:
    //friend std::ostream& operator<< <>(std::ostream& os, const M<d>& m);

    M(const concepts::Real* m);
    virtual ~M() {delete[] m_;}

    template<class F>
    void mult(const F* src, F* dst) const;
    template<class F>
    void mult_T(const F* src, F* dst) const;
    uint n() const {return d;}
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    concepts::Real* m_;
  };

  template<uint d>
  std::ostream& M<d>::info(std::ostream& os) const {
    os << "aglowav::M<" << d << ">=[";
    for(uint i = 0; i < d-1; i++)  os << m_[i] << ", ";
    if (d)  os << m_[d-1];
    return os << "]";
  }

  template <uint d>
  inline M<d>::M(const concepts::Real* m) {
    m_ = new concepts::Real[d];
    concepts::Real A = 0;
    uint i;
    for(i = 0; i < d; i++)  A += m[i];
    for(i = 0; i < d; i++)  m_[i] = sqrt(m[i]/A);
  }

  template <>
  template <>
  inline void M<2>::mult<concepts::Real>(const concepts::Real* src,
           concepts::Real* dst) const {
    dst[0] = m_[0]*src[0] - m_[1]*src[1];
    dst[1] = m_[1]*src[0] + m_[0]*src[1];
  }

  template <>
  template <>
  inline void M<2>::mult_T<concepts::Real>(const concepts::Real* src,
             concepts::Real* dst) const {
    dst[0] =  m_[0]*src[0] + m_[1]*src[1];
    dst[1] = -m_[1]*src[0] + m_[0]*src[1];
  }

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

  template <class F = concepts::Real, uint d = 2>
  class Haar3dXXX : public concepts::Element<F> {
  public:
    class Key {
    public:
      Key(uint l, uint j) : l_(l), j_(j) {}

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

    Haar3dXXX(const Haar3dXXX<F, d>::Key& key,
        const concepts::Real* m) : key_(key), M_(m) {}
    virtual ~Haar3dXXX() {};

    virtual Haar3dXXX<F,d>* child(uint j) const = 0;
    uint nchild() const {return d;}

    virtual const bem::Constant3d001<F>* element(uint j) const = 0;
    Key key() const {return key_;}
    const M<d>& trafoM() const {return M_;}
    virtual const concepts::Real3d& center() const = 0;
    virtual concepts::Real radius() const = 0;

  private:
    Key key_;
    M<d> M_;
  };

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

  template <class F = concepts::Real>
  class Haar3d000 : public Haar3dXXX<F, 2> {
  public:
    Haar3d000(const typename Haar3dXXX<F, 2>::Key& key,
        const concepts::Real* m, const BiClNode00<F>* nd,
        uint idx[], Haar3d000<F>* lft = 0, Haar3d000<F>* rght = 0)
      : Haar3dXXX<F, 2>(key, m), nd_(nd),
        T_((key.l()==0) ? 2 : 1, (key.l()==0) ? 2 : 1, idx), lft_(lft),
        rght_(rght) {}

    inline Haar3d000<F>* child(uint j) const {
      return (!j) ? lft_ : (j == 1) ? rght_ : (Haar3d000<F>*)0;
    }
    inline Haar3d000<F>* left() const {return lft_;}
    inline Haar3d000<F>* right() const {return rght_;}

    inline const concepts::TMatrixBase<F>& T() const {return T_;}

    inline const bem::Constant3d001<F>* element(uint j) const;
    inline const bem::Constant3d001<F>* elmleft() const;
    inline const bem::Constant3d001<F>* elmright() const;

    const concepts::Real3d& center() const {return nd_->center();}
    concepts::Real radius() const {return nd_->radius();}

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

  private:
    const BiClNode00<F>* nd_;
    concepts::TIndex<F> T_;
    Haar3d000<F>* lft_;
    Haar3d000<F>* rght_;
  };

  template<class F>
  inline const bem::Constant3d001<F>* Haar3d000<F>::element(uint j) const {
    return (!j && !lft_) ? nd_->child(0)->element() :
      (j==1 && !rght_) ? nd_->child(1)->element() :
      (const bem::Constant3d001<F>*)0;
  }

  template<class F>
  inline const bem::Constant3d001<F>* Haar3d000<F>::elmleft() const {
    return (lft_) ? 0 : nd_->child(0)->element();
  }

  template<class F>
  inline const bem::Constant3d001<F>* Haar3d000<F>::elmright() const {
    return (rght_) ? 0 : nd_->child(1)->element();
  }

  template<class F>
  std::ostream& Haar3d000<F>::info(std::ostream& os) const {
    os << "aglowav::Haar3d000(" << *nd_ << ", key = (" << this->key().l();
    os << ',' << this->key().j() << "), " << T_ << ", ";
    return os << this->trafoM() << ")";
  }

} // namespace aglowav

#endif // aglowavElement_hh