meshGraph.hh

#ifndef meshGraph_hh
#define meshGraph_hh
 
#include "hp2D.hh"
#include "toolbox.hh"
 
//evtl hp2D
namespace concepts{
 
 
 
 
class MacroElementNode;
 
class MeshGraph2_Edge : public OutputOperator{
 
public:
  MeshGraph2_Edge(MacroElementNode *node1, MacroElementNode *node2){
    node1_ = node1;
    node2_ = node2;
    }
 
  //returns the requested node
  // i = 0 first, else the second
  MacroElementNode* node(ushort i){
    return i ? node2_ : node1_;
  }
 
//    Vertex* getOrigin() {return origin;}
//    Vertex* getDestination() {return destination;}
 
private:
    MacroElementNode* node1_;
    MacroElementNode* node2_;
};
 
//TODO: Ein macroElementNode kennt seine mögliche nachbaranzahl
// viereck = max 4, dreieck = max 3, usw
 
class MacroElementNode : public OutputOperator{
 
public:
 
  MacroElementNode():content_(0), south_(0), east_(0), north_(0), west_(0){};
 
  MacroElementNode(const hp2D::Quad<Real>* quad): content_(0), south_(0), east_(0), north_(0), west_(0){
    content_.push_back(quad);
  };
 
  //sets the pointer to edge in given direction :
  // direction  = 0 south
  //        = 1 east
    //        = 2 north
  //        = 3 west
  //
  // position of the node in the edge
  void setEdge(uint direction, MeshGraph2_Edge* edge, uint pos){
    switch(direction){
    case 0 : south_ = edge; s_p = pos; break;
    case 1 :  east_ = edge; e_p = pos; break;
    case 2 : north_ = edge; n_p = pos; break;
    case 3 :  west_ = edge; w_p = pos; break;
    default:
      throw conceptsException(concepts::MissingFeature("Only 4 direction supported."));
      break;
    }
  }
 
  void addContent(const hp2D::Quad<Real>* quad){
    content_.push_back(quad);
  }
 
  const uint nEdges() const {
    return bool(south_) + bool(east_) + bool(north_) + bool(west_);
  }
 
 
  const hp2D::Quad<Real>* get() const{
    if(content_.size() == 1)
      return content_[0];
    else{
      DEBUGL(1, "content_ "<< content_);
      throw conceptsException(concepts::MissingFeature("Only one content supported."));
    }
  }
 
  //returns pointer to requested neighbour //TODO generalize with "sides of node element
  //direction :south = 0 east = 1 north = 2 west = 3
  const MacroElementNode* neighbour(ushort direction) const{
    switch (direction) {
    case 0:
      if (south_)
        return south_->node((s_p + 1) % 2);
      else
        return 0;
      break;
    case 1:
      if (east_)
        return east_->node((e_p + 1) % 2);
      else
        return 0;
      break;
    case 2:
      if (north_)
        return north_->node((n_p + 1) % 2);
      else
        return 0;
      break;
    case 3:
      if (west_)
        return west_->node((w_p + 1) % 2);
      else
        return 0;
      break;
    default:
      throw conceptsException(concepts::MissingFeature("Only 4 direction supported."));
      break;
    }
  } // neighbour
 
protected:
  virtual std::ostream& info(std::ostream& os) const {
    os << "MacroElementNode[ #content = " << content_.size() << std::endl;
    if (content_.size() == 1){
      os  << " has south : " << ((south_) ? "true" : "false" )<< std::endl
        << " has east  : " << ((east_) ? "true" : "false" )<< std::endl
        << " has north : " << ((north_) ? "true" : "false") << std::endl
        << " has west  : " << ((west_) ? "true" : "false") << std::endl;
    }
    if(content_.size()==1){
      os << "content key = "<< content_[0]->support().key() << std::endl;
      if(south_)
        os << "south key = " << neighbour(0)->get()->support().key() << std::endl;
      if(east_)
        os << "east  key = " << neighbour(1)->get()->support().key() << std::endl;
      if(north_)
        os << "north key = " << neighbour(2)->get()->support().key() << std::endl;
      if(west_)
        os << "west key = " << neighbour(3)->get()->support().key() << std::endl;
    }
    return os << "]";
  }
 
private:
 
  MeshGraph2_Edge *south_, *east_, *north_, *west_;
  //position of the node on the Edge
  uint s_p, e_p, n_p, w_p;
 
  //information holded in that MacroNode, i.e. all Quads in it
  concepts::Sequence<const hp2D::Quad<Real>* > content_;
 
 
};
 
 
 
//a bidirected graph to walk trough a mesh possible having hanging nodes
//therefore socallded MacroElementNodes are introduced
 
//      only half splitting
//      _____
//     |  |  |
//  ___|__|__|__ ...
// |   |     |__|..:   bisection in both direction
// |___|_____|__|..:
//     //////
//     BDN
//
 
// definition
//   regular graph  : only for same refinement level neighbours
//  irregular graph :
 
 
 
 
class MeshGraph2 : public OutputOperator{
 
public:
 
  typedef concepts::ElementAndFacette<hp2D::Element<Real> > UnderlyingElement;
 
  MeshGraph2(const concepts::SpaceOnCells<Real>& spc, bool regular = true) : regular_(regular){
 
    //maps from topology edge key to underlying quads
    concepts::HashMap<concepts::Sequence<UnderlyingElement> > uelm;
    // O(1) access to quads later to avoid dynamic cast
    concepts::HashMap<const hp2D::Quad<Real>*> quads;
 
    //build edge topology map;
    std::auto_ptr<concepts::SpaceOnCells<Real>::Scanner> sc(spc.scan());
    while (*sc) {
      concepts::ElementWithCell<Real>& elm = (*sc)++;
      //TODO: what quad?
      const hp2D::Quad<Real>* quad = dynamic_cast<const hp2D::Quad<Real>*> (&elm);
      //element must be q quad
      conceptsAssert(quad, concepts::Assertion());
 
      quads[quad->support().key()] = quad;
 
      //Build information about underlying elements
      const concepts::Quad& cntr = quad->support();
      for (uint i = 0; i < 4; ++i) {
        const concepts::Edge& edge = *cntr.edge(i);
        uelm[edge.key()].push_back(UnderlyingElement(quad, i));
      }
    }//while
 
 
    //TODO: USE THE FOLLOWING CONCEPT OF BOUNDARY CONDITIONS
    //Space with a boundary condition
    const hp2D::hpAdaptiveSpace<Real>* bcSpace = dynamic_cast<const hp2D::hpAdaptiveSpace<Real>* >(&spc);
    if(!bcSpace)
        throw conceptsException(concepts::MissingFeature("Just Spaces with boundary conditions are allowed"));
 
    const BoundaryConditions* bc = bcSpace->helper().bc();
 
 
 
    //Build the Graph
 
    //build a regular graph, i.e.
    //  - neighboured MacroElementNodes have only one content size one
    //  - its content has same level in h-refinement sense
    if(regular){
 
      //iterate over all Quads and Build the Node contents
      concepts::HashMap<const hp2D::Quad<Real>*>::const_iterator qIter = quads.begin();
      for( ; qIter != quads.end(); ++qIter){
        uint K = qIter->first;
        nodes_[K] = new MacroElementNode(qIter->second);
      }
 
      concepts::HashMap<concepts::Sequence<UnderlyingElement> >::const_iterator
          iter = uelm.begin();
      for( ; iter != uelm.end(); ++iter){
        if(iter->second.size() == 2){
          uint E = iter->first;
          const hp2D::Quad<Real>* quad_1 =
              dynamic_cast<const hp2D::Quad<Real>*> (iter->second[0].elm);
          const hp2D::Quad<Real>* quad_2 =
              dynamic_cast<const hp2D::Quad<Real>*> (iter->second[1].elm);
          uint K1 = quad_1->support().key();
          uint K2 = quad_2->support().key();
          //build a graph Edge between both Nodes
          edges_[E] = new MeshGraph2_Edge(nodes_[K1],nodes_[K2]);
 
          //add the Edge its Nodes
          nodes_[K1]->setEdge(iter->second[0].k, edges_[E], 0);
          nodes_[K2]->setEdge(iter->second[1].k, edges_[E], 1);
        }
      }
 
 
 
    }
    //TODO : build irregular graph, i.e.
    //  - neighboured MacroElementNodes are MacroElements consisting of
    //    quads forming the Edge Support
    //  - neighboured refinement level may differs
    else{
      throw conceptsException(concepts::MissingFeature("Irregular Mapping not implemented"));
    }
  }
 
 
  virtual ~MeshGraph2(){
 
    HashMap<MacroElementNode* >::iterator nodeIter = nodes_.begin();
    HashMap<MeshGraph2_Edge* >::iterator edgeIter = edges_.begin();
    //delete allocated nodes
    for( ; nodeIter != nodes_.end(); ++nodeIter){
      if( nodeIter->second ){
        delete nodeIter->second;
        nodeIter->second = 0;
      }
    }
    //delete allocated edges
    for( ; edgeIter != edges_.end(); ++edgeIter){
      if(edgeIter->second){
        delete edgeIter->second;
        edgeIter->second = 0;
      }
    }
  }
 
  //returns the requested node (representation of Cell with key K
  const MacroElementNode* getNode(const uint K) const{
    HashMap<MacroElementNode* >::const_iterator iter;
    if( (iter = nodes_.find(K)) != nodes_.end() )
      return iter->second;
    return 0;
  }
 
 
 
 
 
protected:
  virtual std::ostream& info(std::ostream& os) const {
    os << "MeshGraph2[ #Nodes = " << nodes_.size() << " #Edges = " << edges_.size()
       << "  regular = "<< regular_ << std::endl;
 
    HashMap<MacroElementNode* >::const_iterator nodeIter =  nodes_.begin();
    for( ; nodeIter != nodes_.end(); ++nodeIter)
      os << *(nodeIter->second) << std::endl;
 
    return os;
 
  }
 
private:
  bool regular_;
 
  //HashMap<MacroElementNode> map_;
 
  //the edge number here corresponds with the coarsest adjacent topological
  // edge that is in the space, i.e. hanging edges
  //TODO: change this to a PList, as no direct contact as for nodes will be needed?!
  HashMap<MeshGraph2_Edge* > edges_;
  //the node corresponds numerical with its topological key
  HashMap<MacroElementNode* > nodes_;
 
};
 
 
} //namespace
 
 
#endif // meshGraph_hh