cluster/tree.hh

Go to the documentation of this file.

/* Cluster tree
 */

#ifndef clusterTree_hh
#define clusterTree_hh

#include <iostream>

#include "basics/typedefs.hh"
#include "basics/vectorsMatricesForward.hh"
#include "basics/exceptions.hh"
#include "space/space.hh"
#include "bem/element.hh"

namespace cluster {

  // *************************************************************** cebysev **

  void cebysev(const concepts::Real3d& x0, const concepts::Real3d& x1,
         concepts::Real3d& c, concepts::Real& r);
  void cebysev(const concepts::Real3d& x0, const concepts::Real3d& x1,
         const concepts::Real3d& x2, concepts::Real3d& c,
         concepts::Real& r);

  // ***************************************************************** BBall **

  template<class F>
  class BBall {
  public:
    virtual void operator()(const concepts::Element<F>& elm,
          concepts::Real3d& c, concepts::Real& r) const = 0;
  };

  // ************************************************************** BBall000 **

  template<class F>
  class BBall000 : public BBall<F> {
  public:
    void operator()(const concepts::Element<F>& elm, concepts::Real3d& c,
        concepts::Real& r) const;
  };

  template<class F>
  void BBall000<F>::operator()(const concepts::Element<F>& elm,
             concepts::Real3d& c, concepts::Real& r) const {

    const bem::Linear3d000<F>* elmA =
      dynamic_cast<const bem::Linear3d000<F>*>(&elm);
    if (elmA != 0) {
      cebysev(elmA->vertex(0), elmA->vertex(1), elmA->vertex(2), c, r);
      return; 
    }
  
    const bem::Constant3d000<F>* elmB =
      dynamic_cast<const bem::Constant3d000<F>*>(&elm);
    if (elmB != 0) {
      cebysev(elmB->vertex(0), elmB->vertex(1), elmB->vertex(2), c, r);
      return;
    }
  
    const bem::Constant3d001<F>* elmC =
      dynamic_cast<const bem::Constant3d001<F>*>(&elm);
    if (elmC != 0) {
      cebysev(elmC->vertex(0), elmC->vertex(1), elmC->vertex(2), c, r);
      return;
    }

    const bem::Constant3d002<F>* elmD =
      dynamic_cast<const bem::Constant3d002<F>*>(&elm);
    if (elmD != 0) {
      cebysev(elmD->vertex(0), elmD->vertex(1), elmD->vertex(2), c, r);
      return;
    }
  
    const bem::Dirac3d000<F>* elmE =
      dynamic_cast<const bem::Dirac3d000<F>*>(&elm);
    if (elmE != 0) {
      switch (elmE->T().n()) {
      case 1:
  c = elmE->vertex(0); r = 0.0;
  break;
      case 2: 
  cebysev(elmE->vertex(0), elmE->vertex(1), c, r);
  break;
      case 3: 
  cebysev(elmE->vertex(0), elmE->vertex(1), elmE->vertex(2), c, r);
  break;
      }
      return;
    }

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

  // *************************************************************** Cluster **

  template<typename F>
  class Cluster;

  template<typename F>
  std::ostream& operator<<(std::ostream& os, const Cluster<F>& clst);

  template<class F>
  class Cluster {
  public:
    typedef F CF;

    Cluster(const concepts::Element<F>& elm, const BBall<F>& ball,
      Cluster<F>* lnk);
    Cluster(uint idx, Cluster<F>& chld, const concepts::Real3d& c,
      concepts::Real r, uint nlf = 0);

    inline Cluster<F>* child(uint j) const;
    inline Cluster<F>* child() const;
    inline Cluster<F>* link() const {return lnk_;}
    inline Cluster<F>*& link() {return lnk_;}
    inline const concepts::Element<F>* element() const;
    inline const concepts::Real3d& center() const {return c_;}
    inline concepts::Real radius() const {return r_;}
    inline uint nleaf() const {return nlf_;}
    inline uint index() const {return idx_;}
    std::ostream& info(std::ostream& os) const;

  private:
    friend std::ostream& operator<< <>(std::ostream& os,
               const Cluster<F>& clst);
    union CE {
      Cluster<F>*                 chld;
      const concepts::Element<F>* elm;
    };

    Cluster<F>* lnk_;
    CE          ce_;

    uint idx_;
    uint nlf_;

    concepts::Real3d c_;
    concepts::Real   r_;
  };

  template<class F>
  inline Cluster<F>* Cluster<F>::child(uint j) const {
    Cluster<F>* c = child();
    while (c && j)  {c = c->link();  j--;}
    return (j) ? (Cluster<F>*)0 : c;
  }

  template<class F>
  inline Cluster<F>* Cluster<F>::child() const {
    return (nlf_) ? ce_.chld : (Cluster<F>*)0;
  }

  template<class F>
  inline const concepts::Element<F>* Cluster<F>::element() const {
    return (nlf_) ? (concepts::Element<F>*)0 : ce_.elm;
  }

  template<class F>
  inline std::ostream& operator<<(std::ostream& os, const Cluster<F>& clst) {
    return clst.info(os);
  }

  template<class F>
  Cluster<F>::Cluster(uint idx, Cluster<F>& chld, const concepts::Real3d& c,
          concepts::Real r, uint nlf)
    : lnk_(0), idx_(idx), nlf_(nlf), c_(c), r_(r) {

    ce_.chld = &chld;

    if (nlf_ == 0)
      for(const Cluster<F>* c = &chld; c != 0; c = c->lnk_)
  nlf_ += c->nlf_ ? c->nlf_ : 1;
  }

  template<class F>
  Cluster<F>::Cluster(const concepts::Element<F>& elm,
          const BBall<F>& ball, Cluster<F>* lnk)
    : lnk_(lnk), idx_(0), nlf_(0) {

    ce_.elm = &elm;
    ball(elm, c_, r_);
  }

  template<class F>
  std::ostream& Cluster<F>::info(std::ostream& os) const {
    os << "Cluster(idx = " << idx_ << ", r = " << r_ << ", c = " << c_;
    os << ", nleaf = " << nlf_ << ", lnk = " << lnk_;
    os << ", chld = " << child() << ')';
    return os;
  }

  // ************************************************************ TreeTraits **

  template<class Node>
  class TreeTraits {
  public:
    typedef typename Node::CF F;

    static Node* child(const Node* nd, uint j) {return nd->child(j);}
    static const concepts::Element<F>* element(const Node* nd) {
      return nd->element();
    }
    static const concepts::Real3d& center(const Node* nd) {
      return nd->center();}
    
    static concepts::Real radius(const Node* nd) {
      return nd->radius();
    }
    static uint nleaf(const Node* nd) {return nd->nleaf();}
    static uint index(const Node* nd) {return nd->index();}
  };

  // ****************************************************************** Tree **
  
  template<class CNode>
  class Tree;
  
  template<class CNode>
  inline std::ostream& operator<<(std::ostream& os, const Tree<CNode>& t) {
    return t.info(os);
  }

  template<class CNode>
  class Tree {
    friend std::ostream& operator<< <>(std::ostream& os, const Tree<CNode>& t);

  public:
    typedef TreeTraits<CNode> Traits;
    typedef typename Traits::F F;
    typedef CNode Node;

    inline Tree(const concepts::Space<F>& spc) : spc_(spc) {}
    virtual ~Tree() {}

    virtual std::ostream& info(std::ostream& os) const;
    virtual uint nclst() const = 0;
    virtual uint nleaf() const = 0;
    virtual const Node* root() const = 0;
    inline const concepts::Space<F>& space() const {return spc_;}

 private:
    const concepts::Space<F>& spc_;

    void info_(const Node* clst, concepts::Real& rhomin,
         concepts::Real& rhomax, concepts::Real& taumin,
         concepts::Real& taumax, uint& h) const;
  };

  template<class CNode>
  void Tree<CNode>::info_(const CNode* clst, concepts::Real& rhomin,
        concepts::Real& rhomax, concepts::Real& taumin,
        concepts::Real& taumax, uint& h) const {
    uint j = 0;
    const CNode* chld = Traits::child(clst, j);

    while (chld) {
      if (Traits::nleaf(chld)) {
  uint hh = 0;
  info_(chld, rhomin, rhomax, taumin, taumax, hh);

  if (hh > h) h = hh;

  concepts::Real rho = Traits::radius(chld) / Traits::radius(clst);
  if (rho < rhomin) rhomin = rho;
  if (rho > rhomax) rhomax = rho;

  concepts::Real tau =
    (concepts::Real)Traits::nleaf(chld) / Traits::nleaf(clst);
  if (tau < taumin) taumin = tau;
  if (tau > taumax) taumax = tau;

  if (rho > 2) {
    std::cout << *clst << std::endl;
    std::cout << *chld << std::endl;
  }

      }
      chld = Traits::child(clst, ++j);
    }
    ++h;
  }

  template<class CNode>
  std::ostream& Tree<CNode>::info(std::ostream& os) const {

    concepts::Real rhomin = 1, rhomax = 0;
    concepts::Real taumin = 1, taumax = 0;
    uint h = 0;

    info_(root(), rhomin, rhomax, taumin, taumax, h);

    os << "cluster::Tree(nclst = " << nclst() << ", nleaf = " << nleaf();
    os << ", h = " << h;

    if (h > 1) {
      os << ", rho in [" << rhomin << ", " << rhomax;
      os << "], tau in [" << taumin << ", " << taumax << "])";
    } else {
      os << ')';
    }

    return os;
  }

} // namespace cluster

#endif // clusterTree_hh