space/domainDecomp.hh

Go to the documentation of this file.

#ifndef spcDomainDecomp_hh
#define spcDomainDecomp_hh

#include "basics/debug.hh"
#include "geometry/boundaryConditions.hh"
#include "geometry/cellConditions.hh"
#include "geometry/mesh.hh"
#include "space/space.hh"
#include "space/spaceSet.hh"
#include "toolbox/hashMap.hh"
#include "toolbox/scannerConnectors.hh"
#include "toolbox/sequence.hh"
#include "toolbox/set.hh"

#define DDSpaceConstr_D 0
#define DDSpaceDestr_D 0
#define DDSpaceRebuild_D 0
#define DDSpaceGetIndices_D 0
#define DDSpaceInfoIndices_D 0

namespace concepts {

  // *************************************************************** DDSpace **

  template<class F>
  class DDSpace : public Space<F> {
  public:
    DDSpace(uint domains = 0) : domains_(domains) {}
    virtual ~DDSpace() {}
    const uint domains() const { return domains_; }
    virtual const Space<F>& space(uint i) const = 0;
    virtual const Set<IndexRange> indicesI(uint i) const = 0;
//     virtual const Sequence<Set<IndexRange> >& indicesI() const {
//       return indicesI_; }
    virtual const Set<IndexRange> indicesB(uint i) const = 0;
//     virtual const Sequence<Set<IndexRange> >& indicesB() const { 
//       return indicesB_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;

    uint domains_;
    Sequence<Set<IndexRange> > indicesI_, indicesB_;
  };

  template<class F>
  std::ostream& DDSpace<F>::info(std::ostream& os) const {
    return os << "DDSpace(" << domains_ << " domains)";
  }

  // ********************************************************** DomainDecomp **

  template<class F>
  class DomainDecomp : public DDSpace<typename F::t_type>, public Subspace {
  public:
    typedef concepts::Element<typename F::t_type> Element;
    typedef concepts::Scan<Element> Scan;

    template<class G>
    DomainDecomp(G& prebuild, Sequence<Set<uint> > domains,
                 BoundaryConditions* bc = 0, CellConditions* cc = 0,
                 uint spcNo = 0, uint* offset = 0);
    virtual ~DomainDecomp();

    void rebuild();
    bool available() const;

    inline virtual uint dim() const;
    inline virtual uint nelm() const;
    inline virtual Scan* scan() const;
//     inline Scan* scan();

    virtual uint& lastIdx();
    virtual uint offset() const;
    virtual const F& space(uint i) const;
    F& space(uint i);
    virtual const Set<IndexRange> indicesI(uint i) const;
//     const Sequence<Set<IndexRange> >& indicesI() const { return indicesI_; }
    virtual const Set<IndexRange> indicesB(uint i) const;
//     const Sequence<Set<IndexRange> >& indicesB() const { return indicesB_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    uint spcNo_;
    Sequence<F*> spaces_;
    Array<uint> spcBuild_;
    Sequence<CellConditions> cc_;

    HashMap<uint> attrToDomain_;

    concepts::Joiner<Element*, 1>* elm_;
    uint nelm_;

    uint idx_;

    template<class G, class H>
    void getIndices_(const G& prebuild, const H& cntr,
                     Set<IndexRange>& indices);
  };

  template<class F>
  template<class G>
  DomainDecomp<F>::DomainDecomp(G& prebuild, Sequence<Set<uint> > domains,
                                BoundaryConditions* bc, CellConditions* cc,
                                uint spcNo, uint* offset) : 
    DDSpace<typename F::t_type>(), spcNo_(spcNo), spcBuild_(0),
    elm_(0), nelm_(0) {
    // collect all attributes of coarse mesh
    Set<uint> attributes, tmp;
    std::auto_ptr<concepts::Scan2> sc(prebuild.mesh().scan());
    while (*sc) {
      const Attribute& attr = (*sc)++.connector().attrib();
      if (!cc || ((*cc)(attr).type() != CellCondition::INACTIVE))
        attributes.insert(attr);
    }
    tmp = attributes;
    DEBUGL(DDSpaceConstr_D, "attributes in mesh = " << attributes);


    // collect all attributes and delete empty domains
    Sequence<Set<uint> >::iterator i = domains.begin();
    for (; i != domains.end(); ) {
      // set difference with valid attributes
      Set<uint> nonvalid = *i - tmp;
      // there are only valid attributes given
      conceptsAssert3(nonvalid.empty(), Assertion(), "non valid attributes " <<
                      nonvalid << " in given domain " << *i);

      if ((*i).empty())
        i = domains.erase(i);             // delete empty domain
      else {
        // contribute to mapping from attribute to domain
        for(Set<uint>::const_iterator j = (*i).begin(); j != (*i).end(); ++j)
          attrToDomain_[*j] = this->domains_;
        tmp = tmp - *i++;                 // erase attributes from set
        ++this->domains_;                 // count non-empty domains
      }
    }

    // rest of attributes to a new domain
    if (!tmp.empty()) {
      domains.push_back(tmp);
      // contribute to mapping from attribute to domain
      for(Set<uint>::const_iterator j = tmp.begin(); j != tmp.end(); ++j)
        attrToDomain_[*j] = this->domains_;
      ++this->domains_;
    }
    
    DEBUGL(DDSpaceConstr_D, "domains = " << domains);

    // set the number of builds to zero
    spcBuild_.resize(this->domains_);
    spcBuild_.zeros();

    i = domains.begin();
    for (; i != domains.end(); ) {

      // create for each domain a entrance
      cc_.insert(cc_.end(), this->domains_, CellConditions());

      i = domains.begin();
      for (Sequence<CellConditions>::iterator j = cc_.begin(); j != cc_.end();
           ++j, ++i) {
        // set difference with valid attributes
        const Set<uint> nonactiv = attributes - *i;
        for (Set<uint>::const_iterator k = nonactiv.begin();
             k != nonactiv.end(); ++k)
          (*j).add(*k, CellCondition::INACTIVE);
      }
    }
    DEBUGL(DDSpaceConstr_D, "cell conditions of domains = " << cc_);

    // create spaces
    uint* idx = 0;
    F* spc;
    for (Sequence<CellConditions>::iterator j = cc_.begin(); j != cc_.end(); ++j)
      {
        spaces_.push_back(spc = new F(prebuild, bc, &*j, spcNo_, offset, idx));
        idx = &spc->lastIdx();
      }
    DEBUGL(DDSpaceConstr_D, "spaces of domains = " << spaces_);
    
    // reserve space for index sets for each sub domains
    this->indicesI_.resize(this->domains_, Set<IndexRange>());
    this->indicesB_.resize(this->domains_, Set<IndexRange>());
  }

  template<class F>
  DomainDecomp<F>::~DomainDecomp() {
    // deletes in reverse order, because the first space holds the
    // pointer to the index
    uint k = this->domains();
    typename Sequence<F*>::reverse_iterator i = spaces_.rbegin();
    for (; i != spaces_.rend(); ++i) {
      DEBUGL(DDSpaceDestr_D, "Try to delete Space(" << k-- << ") = " << **i);
      delete *i;
    }
    concepts::Joiner<Element*, 1>::destructor(elm_, false);
    DEBUGL(DDSpaceDestr_D, "done");
  }

  template<class F>
  std::ostream& DomainDecomp<F>::info(std::ostream& os) const {
    os << "DomainDecomp(" << this->domains_ << " domains, ";
    if (!available())
      os << "not built";
    else {
      os << "dim = " << dim() << ", nelm = " << nelm() << ", ";
      for (uint i = 0; i < spaces_.size();) {
        os << "space(" << i << ") = " << space(i);
#if DDSpaceInfoIndices_D
        os << ", I(" << i << ") = " << indicesI(i)
           << ", B(" << i << ") = " << indicesB(i);
#else
        os << ", I(" << i << ").dim = " << indicesI(i).dim()
           << ", B(" << i << ").dim = " << indicesB(i).dim();
#endif
        if (++i < spaces_.size()) os << ", ";
      }
    }
    return os << ")";
  }

  template<class F>
  void DomainDecomp<F>::rebuild() {
    DEBUGL(DDSpaceRebuild_D, "Start rebuilding");
    // If nothing to do, do nothing
    if (available()) {
      DEBUGL(DDSpaceRebuild_D, "space already built - nothing to do");
      return;
    }

    // Remove old list, but not elements. They are removed in spaces
    // itself.
    concepts::Joiner<Element*, 1>::destructor(elm_, false);
    nelm_ = 0;

    // Rebuild all spaces, but don't refresh own array of all
    // elements. That is done in method scan(), because the
    // rebuild()-method of the spaces could be called also directly.
    uint k = 0;
    typename Sequence<F*>::iterator i = spaces_.begin();
    // reset the index counter (identical for all spaces/domains)
    (*i)->lastIdx() = (*i)->offset();
    // clear all indices
    (*i)->prebuild().clearAllIndices(spcNo_);
    uint* b = spcBuild_;
    for (; i != spaces_.end(); ++i) {
      (*i)->rebuild(true);
      DEBUGL(DDSpaceRebuild_D, "Space(" << k++ << ") = " << **i);
      // space should be available now
      conceptsAssert((*i)->available().first, Assertion());
      // memorize number of the build
      *b++ = (*i)->available().second;
      conceptsAssert(nelm_ == 0 || (*i)->nelm() == nelm_, Assertion());
      nelm_ = (*i)->nelm();
      DEBUGL(DDSpaceRebuild_D, "nelm = " << nelm_);
    }

    // fill list of elements
    uint j = 0;
    std::auto_ptr<typename Space<typename F::t_type>::Scanner> sc(0);
    for (i = spaces_.begin(); i != spaces_.end(); ++i) {
      sc.reset((*i)->scan());
      while (*sc) {
        Element* elm = &(*sc)++;
        DEBUGL(DDSpaceRebuild_D, "Element = " << elm);
        DEBUGL(DDSpaceRebuild_D, "Element = " << *elm);
        if (elm->T().n() > 0) {
          DEBUGL(DDSpaceRebuild_D, "Element(" << j << ") = " << *elm);
          elm_ = new concepts::Joiner<Element*, 1>(elm, elm_);
          ++j;
        }
      }
    }
    conceptsAssert(j == nelm_, Assertion());

    for(uint j = 0; j < this->domains_; ++j) {
      this->indicesI_[j].clear();
      this->indicesB_[j].clear();
    }

    // collect all index ranges into this set to determine if a index
    // range occurs at least two times
    Set<IndexRange> indices;
    std::auto_ptr<concepts::Scan2> se(spaces_[0]->prebuild().mesh().scan());
    while (*se)
      getIndices_(spaces_[0]->prebuild(), (*se)++.connector(), indices);

#if DDSpaceRebuild_D
    for(uint j = 0; j < this->domains_; ++j) {
      DEBUGL(1, "indicesI_[" << j << "] = " << this->indicesI_[j]);
      DEBUGL(1, "indicesB_[" << j << "] = " << this->indicesB_[j]);
    }

    i = spaces_.begin(); k = 0;
    for (; i != spaces_.end(); ++i)
      DEBUGL(1, "Space(" << k++ << ") = " << **i);
#endif
  }

  template<class F>
  bool DomainDecomp<F>::available() const {
//     uint build;
    const uint* b = spcBuild_;
    typename Sequence<F*>::const_iterator i = spaces_.begin();
    while(i != spaces_.end() && (*i)->available().first && 
    (*i)->available().second == *b++) ++i;
    return i == spaces_.end();
  }

  template<class F>
  uint DomainDecomp<F>::dim() const {
    if (!available())
      throw conceptsException(concepts::SpaceNotBuilt());
    // dimension is for the spaces of all domains the same, they live
    // on the same index range
    return (*spaces_.begin())->dim();
  }

  template<class F>
  uint DomainDecomp<F>::nelm() const {
    if (!available())
      throw conceptsException(concepts::SpaceNotBuilt());
    return nelm_;
  }

  template<class F>
  typename DomainDecomp<F>::Scan* DomainDecomp<F>::scan() const {
    if (!available())
      throw conceptsException(concepts::SpaceNotBuilt());
    return new concepts::PListScan<Element>(*elm_);
  }

  template<class F>
  uint& DomainDecomp<F>::lastIdx() {
    return (*spaces_.begin())->lastIdx();
  }

  template<class F>
  uint DomainDecomp<F>::offset() const {
    return (*spaces_.begin())->offset();
  }

  template<class F>
  const F& DomainDecomp<F>::space(uint i) const {
    conceptsAssert(i < this->domains_, Assertion());
    return *spaces_[i];
  }

  template<class F>
  F& DomainDecomp<F>::space(uint i) {
    conceptsAssert(i < this->domains_, Assertion());
    return *spaces_[i];
  }

  template<class F>
  const Set<IndexRange> DomainDecomp<F>::indicesI(uint i) const {
    conceptsAssert(i < this->domains_, Assertion());
    if (!available())
      throw conceptsException(concepts::SpaceNotBuilt());
    return this->indicesI_[i];
  }

  template<class F>
  const Set<IndexRange> DomainDecomp<F>::indicesB(uint i) const {
    conceptsAssert(i < this->domains_, Assertion());
    if (!available())
      throw conceptsException(concepts::SpaceNotBuilt());
    return this->indicesB_[i];
  }

  template<class F>
  template<class G, class H>
  void DomainDecomp<F>::getIndices_(const G& prebuild, const H& cntr,
                                    Set<IndexRange>& indices)
  {
    DEBUGL(DDSpaceGetIndices_D, "cntr = " << cntr);
    HashMap<uint>::const_iterator i = attrToDomain_.find(cntr.attrib());
    conceptsAssert(i != attrToDomain_.end(), Assertion());
    // domain
    uint d = i->second;
    DEBUGL(DDSpaceGetIndices_D, "domain = " << d);

    Set<IndexRange> &indicesB = this->indicesB_[d],
      &indicesI = this->indicesI_[d];

    uint dim = 0;
    try {
      while(1) {
        // Index ranges of this cell of entities of dimension dim. 
        // Throws an exception if dimension is to high and leaves the
        // while loop.
        const Set<IndexRange> idx = prebuild.indices(dim, cntr, spcNo_);
        if (!idx.empty()) {
          DEBUGL(DDSpaceGetIndices_D,
                 cntr.key() << " - dim = " << dim << ", idx = " << idx);
          DEBUGL(DDSpaceGetIndices_D, "indices = " << indices);
          // Indices, which occur first time or lie the interior this domain.
          Set<IndexRange> recent = idx - (indices - indicesI);
          // already occured indices
          Set<IndexRange> occured = idx - recent;
          DEBUGL(DDSpaceGetIndices_D,
                 "recent = " << recent << ", occured = " << occured);
          if (!occured.empty()) {
            // its on the boundary of this domain
            indicesB = indicesB || occured;
            DEBUGL(DDSpaceGetIndices_D, "indicesB = " << indicesB);
            // Delete it from index set inside other domains, could be
            // counted there
            for(uint j = 0; j < this->domains_; ++j)
              if (j != d) {
                Set<IndexRange> occuredDomain = occured && this->indicesI_[j];
                this->indicesI_[j] = this->indicesI_[j] - occuredDomain;
                this->indicesB_[j] = this->indicesB_[j] || occuredDomain;
              }
          }
          if (!recent.empty()) {
            indicesI = indicesI || recent;
            DEBUGL(DDSpaceGetIndices_D, "indicesI = " << indicesI);
          }
          indices = indices || idx;
          // set difference with index ranges at boundary are inside domain
          DEBUGL(DDSpaceGetIndices_D,
                 "indicesI_[" << d << "] = " << this->indicesI_[d]);
          DEBUGL(DDSpaceGetIndices_D,
                 "indicesB_[" << d << "] = " << this->indicesB_[d]);
        } // idx not empty
        ++dim;
        DEBUGL(DDSpaceGetIndices_D, "dim = " << dim);
      } // while (1)
    }
    catch (concepts::MissingFeature& e) {}

    DEBUGL(DDSpaceGetIndices_D, "done for " << cntr);
    const H* chld = 0;
    for(uint j = 0; (chld = cntr.child(j)) != 0; ++j)
      getIndices_(prebuild, *chld, indices);
  }

} // namespace concepts

#endif