tree.hh

Go to the documentation of this file.

 
#ifndef clusterTree_hh
#define clusterTree_hh
 
#include <iostream>
 
#include "basics/output.hh"
#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 << concepts::typeOf(*this) << "(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::" << concepts::typeOf(*this) << "(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