Topology.hpp

#ifndef __FIBER_GRID_TOPOLOGY_HPP
#define __FIBER_GRID_TOPOLOGY_HPP

#include "GridAPI.h"
#include <vector/MultiIndex.hpp>


namespace Fiber
{

//      GRID_API BaseSpace : public MemCore::ReferenceBase<BaseSpace> , public MemCore::Intercube

namespace RegularTopology
{

template <int Dims>
inline index_t NumberOfEdges(const MultiIndex<Dims>&NumberOfVertices) throw()
{
index_t Edges = 0; 
        for(int i=0; i<Dims; i++)
        {
                Edges += ( NumberOfVertices - MultiIndex<Dims>::BitIndex(1<<i) ).size();
        }

        return Edges;
}

template <int Dims>
inline int      EdgeOrientation(index_t EdgeIndex, const MultiIndex<Dims>&NumberOfVertices) throw(int)
{
        for(int Orientation=0; Orientation<Dims; Orientation++)
        {
        index_t NumberOfOrientedEdges = (NumberOfVertices - MultiIndex<Dims>::Axis(Orientation )).size(); 
                if (EdgeIndex < NumberOfOrientedEdges)
                {
                        return Orientation;
                } 
                EdgeIndex -= NumberOfOrientedEdges;
        } 
        throw -1;
}

template <int Dims>
void ComputeFirstEdgeVertex(MultiIndex<Dims>&FirstVertex, int&Orientation,
                            index_t EdgeIndex, const MultiIndex<Dims>&NumberOfVertices) throw(int)
{
//      assert( EdgeIndex < NumberOfEdges( NumberOfVertices ) ); 

        for(Orientation=0; Orientation<Dims; Orientation++)
        {
        MultiIndex<Dims> NumberOfOrientedEdges = NumberOfVertices - MultiIndex<Dims>::Axis(Orientation );
        index_t NumberOfOrientedEdgesCount = NumberOfOrientedEdges.size(); 
                if (EdgeIndex < NumberOfOrientedEdgesCount)
                {
                        FirstVertex = MultiIndex<Dims>(EdgeIndex, NumberOfOrientedEdges, CreateFromLinear() );
                        return;
                } 
                EdgeIndex -= NumberOfOrientedEdgesCount;
        } 
        throw -1;
}

template <int Dims>
inline MultiIndex<Dims> ComputeSecondEdgeVertex(const MultiIndex<Dims>&Vertex, int Orientation) throw()
{
        return Vertex + MultiIndex<Dims>::BitIndex( 1 << Orientation );
}

template <int Dims>
inline std::pair<MultiIndex<Dims>, MultiIndex<Dims> >
ComputeEdgeVertices( index_t EdgeIndex, const MultiIndex<Dims>&NumberOfVertices) throw(int)
{
//      assert( EdgeIndex < NumberOfEdges( NumberOfVertices) ); 

MultiIndex<Dims> FirstVertex;
int     Orientation;
        ComputeFirstEdgeVertex(FirstVertex, Orientation,
                               EdgeIndex, NumberOfVertices); 

        return std::pair<MultiIndex<Dims>, MultiIndex<Dims> >(FirstVertex,
                                                         ComputeSecondEdgeVertex(FirstVertex, Orientation)
                );
}

template <int Dims>
inline index_t ComputeEdgeIDfromVertexAndOrientation(const MultiIndex<Dims>&Vertex, int Orientation,
                                                     const MultiIndex<Dims>&NumberOfVertices) throw(int)
{
        assert( Vertex < NumberOfVertices ); 
        assert( Orientation >= 0);

index_t EdgeNumberBias = 0;
        for(int O=0; O<Orientation+1; O++)
        {
        MultiIndex<Dims> NumberOfOrientedEdges = NumberOfVertices - MultiIndex<Dims>::Axis(O);
        index_t NumberOfOrientedEdgesCount = NumberOfOrientedEdges.size(); 

/*
                cout << "      ComputeEdgeIDfromVertexAndOrientation("<<Vertex<<","<<Orientation<<","<<NumberOfVertices<<")" 
                     << "Oriention Axis: " << MultiIndex<3>::Axis(O)
                     << "  O="<<O<<" NumberOfOrientedEdges="<<NumberOfOrientedEdges<<"   NumberOfOrientedEdgesCount="<<NumberOfOrientedEdgesCount
                     << "  EdgeNumberBias="<<EdgeNumberBias<<endl;
*/
                if (O<Orientation)
                {
                        EdgeNumberBias += NumberOfOrientedEdgesCount;
                }
                else
                {
                        return EdgeNumberBias + Vertex.linear(NumberOfOrientedEdges);
                }
        }
        throw EdgeNumberBias;
}

} // namespace RegularTopology

using namespace RegularTopology;

} /* namespace Fiber */ 

#endif /* __FIBER_GRID_TOPOLOGY_HPP */