#ifndef GEO_FAULT_FACTORY_HH
#define GEO_FAULT_FACTORY_HH

#include <math.h>

#include <dune/faultnetworks/facehierarchy.hh>
#include <dune/faultnetworks/utils/debugutils.hh>
#include <dune/faultnetworks/faultfactories/oscunitcube.hh>
#include <dune/faultnetworks/levelinterfacenetwork.hh>
#include <dune/faultnetworks/interfacenetwork.hh>
#include <dune/faultnetworks/faultinterface.hh>
#include <dune/faultnetworks/hierarchicleveliterator.hh>

#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>

#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/boundaryiterator.hh>

class FaultCorridor {
protected:
    int level_;

    int lowerY_;
    int upperY_;

    std::vector<int> faultToLevel_;
    int numberOfFaults_;


    std::vector<int> separatingYIDs_;
    std::vector<int> faultYs_;

private:
    void setup() {
        numberOfFaults_ = 0;
        for (size_t i=0; i<faultToLevel_.size(); i++) {
            if (faultToLevel_[i]==level_)
                numberOfFaults_++;
        }

        if (numberOfFaults_>0) {
            separatingYIDs_.resize(numberOfFaults_-1);
            faultYs_.resize(numberOfFaults_);

            int yStep = (upperY_-lowerY_)/numberOfFaults_;

            if (yStep<2) {
                std::cout << "FaultCorridor::setup() : FaultCorridor does not admit " + std::to_string(numberOfFaults_) + " faults on level " + std::to_string(level_) << std::endl;
                DUNE_THROW(Dune::Exception, "FaultCorridor::setup() : FaultCorridor does not admit " + std::to_string(numberOfFaults_) + " faults on level " + std::to_string(level_));
            }

            int offSet = yStep + (upperY_-lowerY_-yStep*numberOfFaults_)/2;
            int prevYID = lowerY_;

            for (size_t i=0; i<separatingYIDs_.size(); i++) {
                int newYID = lowerY_ + offSet + i*yStep;

                separatingYIDs_[i] = newYID;
                faultYs_[i] = (prevYID+newYID)/2;

                prevYID = newYID;
            }
            faultYs_[faultYs_.size()-1] = (upperY_+prevYID)/2;
        } else {
            separatingYIDs_.resize(0);
            faultYs_.resize(0);
        }
    }


public:
    FaultCorridor() {}

    FaultCorridor(const int level, const int lowerY, const int upperY, std::vector<int> faultToLevel) :
        level_(level),
        lowerY_(lowerY),
        upperY_(upperY),
        faultToLevel_(faultToLevel)
    {
        setup();
    }

    FaultCorridor(const FaultCorridor& faultCorridor) :
        level_(faultCorridor.level_),
        lowerY_(faultCorridor.lowerY_),
        upperY_(faultCorridor.upperY_),
        faultToLevel_(faultCorridor.faultToLevel_),
        numberOfFaults_(faultCorridor.numberOfFaults_),
        separatingYIDs_(faultCorridor.separatingYIDs_),
        faultYs_(faultCorridor.faultYs_)
    {}


    void nextLevel(std::vector<FaultCorridor>& newFaultCorridors, const int resolutionFactor = 2) {
        size_t idx = 0;
        int prevYID = resolutionFactor*lowerY_;

        for (size_t i=0; i<faultYs_.size(); i++) {
            std::vector<int> newFaultToLevel;
            for (size_t j=idx; j<faultToLevel_.size(); j++) {
                idx++;
                if (faultToLevel_[j]>level_)
                    newFaultToLevel.push_back(faultToLevel_[j]);
                else
                    break;
            }

            newFaultCorridors.push_back(FaultCorridor(level_+1, prevYID, resolutionFactor*faultYs_[i], newFaultToLevel));
            prevYID = resolutionFactor*faultYs_[i];
        }

        std::vector<int> newFaultToLevel;
        for (size_t j=idx; j<faultToLevel_.size(); j++) {
            idx++;
            if (faultToLevel_[j]>level_)
                newFaultToLevel.push_back(faultToLevel_[j]);
            else
                break;
        }
        newFaultCorridors.push_back(FaultCorridor(level_+1, prevYID, resolutionFactor*upperY_, newFaultToLevel));
    }

    bool isfaultY(int candidatY, int& idx) const {
        for (size_t i=0; i<faultYs_.size(); i++) {
            if (faultYs_[i]<candidatY)
                continue;

            if (faultYs_[i]==candidatY) {
                idx = i;
                return true;
            } else
                return false;
        }
        return false;
    }

    int upperY() const {
        return upperY_;
    }

    int lowerY() const {
        return lowerY_;
    }

    const std::vector<int>& faultToLevel() const {
        return faultToLevel_;
    }

    const std::vector<int>& separatingYIDs() const {
        return separatingYIDs_;
    }

    const std::vector<int>& faultYs() const {
        return faultYs_;
    }
};


template <class GridType>
class LevelGeoFaultFactory
{
    //! Parameter for mapper class
    template<int dim>
    struct FaceMapperLayout
    {
        bool contains (Dune::GeometryType gt)
        {
            return gt.dim() == dim-1;
        }
    };

    protected:
        static const int dimworld = GridType::dimensionworld;
        static const int dim = GridType::dimension;
        typedef typename GridType::ctype ctype;

        typedef typename GridType::LevelGridView GV;
        typedef typename GridType::LevelIntersection Intersection;

        typedef typename GridType::template Codim<0>::Entity Element;
        static const int dimElement = Element::dimension;

        typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;

    private:
        const std::shared_ptr<GridType> grid_;
        const int level_;
        const ctype resolution_;
        const bool allowLevelIntersections_;
        const GV gridView_;

        const typename GV::IndexSet& indexSet_;

        FaceMapper faceMapper_;

        std::vector<Intersection> faces_;
        std::vector<Dune::FieldVector<ctype, dimworld>> vertexPositions_;
        std::vector<std::vector<size_t>> vertexToFaces_;

    private:

        bool intersectionAllowed(const size_t faceID, const size_t vertexID, const double maxAngle, const std::set<size_t>& separatingDofs, const std::set<int>& separatingYIDs, const std::set<size_t>& faultDofs) const {
            const Intersection& intersection = faces_[faceID];

            //check if "back" edge, cannot deviate from desiredOrientation more than maxAngle degrees (radian)
            Dune::FieldVector<ctype, dimworld> desiredOrientation(1);
            desiredOrientation /= desiredOrientation.two_norm();

            Dune::FieldVector<ctype, dimworld> orientation(intersection.geometry().center());
            orientation -= vertexPositions_[vertexID];
            orientation /= orientation.two_norm();

            if (std::acos(desiredOrientation*orientation) > maxAngle)
                return false;

            // check if intersection has separating dofs or other fault dofs
            std::set<size_t> intersectionDofs;
            computeIntersectionDofs(indexSet_, intersection, intersectionDofs);
            intersectionDofs.erase(vertexID);

            std::set<size_t>::iterator iDofIt = intersectionDofs.begin();
            std::set<size_t>::iterator iDofItEnd = intersectionDofs.end();
            for (; iDofIt!=iDofItEnd; ++iDofIt) {
                int yID = computeYID(vertexPositions_[*iDofIt]);
                if (separatingDofs.count(*iDofIt) || faultDofs.count(*iDofIt) || separatingYIDs.count(yID)) {
                    return false;
                }
            }

            return true;
        }

        /*
        ctype distance(const Intersection& isec, const Dune::FieldVector<ctype, dimworld> & vertex) const {
            Dune::FieldVector<ctype, dimworld> vec(vertex);
            vec += isec.center();

            return isec.unitOuterNormal()*vec;
        }*/

        //works only in 2D, vertex(1) = vertex(0) - yid, intersection of line (given by vertex and derivative 1) with y axis
        int computeYID(Dune::FieldVector<ctype, dimworld> vertex) const {
            return (int) ((vertex[1]-vertex[0])*resolution_);
        }

        void computeIntersectionDofs(const typename GV::IndexSet& idxSet, const Intersection& intersection, std::set<size_t>& intersectionDofs) const {
            intersectionDofs.clear();

            // loop over all vertices of the intersection
            const Element& insideElement = intersection.inside();

            const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
            for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
                size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
                size_t globalIdx = idxSet.subIndex(insideElement, idxInElement, dimElement);

                intersectionDofs.insert(globalIdx);
            }
        }

        bool isTargetBoundaryVertex(const Dune::FieldVector<ctype, dimworld>& vertex) const {
            bool result = false;
            for (size_t i=0; i<vertex.size(); i++) {
                if (vertex[i]==1.0) {
                    result = true;
                    break;
                }
            }
            return result;
        }

        typedef std::vector<double>::const_iterator ConstVectorIterator;

        struct ordering {
            bool operator ()(std::pair<size_t, ConstVectorIterator> const& a, std::pair<size_t, ConstVectorIterator> const& b) {
                return *(a.second) < *(b.second);
            }
        };

        void distanceSort(std::vector<size_t>& admissibleFaces, std::vector<double>& distances) const {
            std::vector<std::pair<size_t, ConstVectorIterator>> order(distances.size());

            size_t j = 0;
            ConstVectorIterator it = distances.begin();
            ConstVectorIterator itEnd = distances.end();
            for (; it != itEnd; ++it, ++j)
                order[j] = std::make_pair(j, it);

            sort(order.begin(), order.end(), ordering());

            std::vector<size_t> initialAdmissibleFaces(admissibleFaces);
            for (size_t i=0; i<admissibleFaces.size(); i++)
                admissibleFaces[i] = initialAdmissibleFaces[order[i].first];
        }

        void generateFaultSeeds(std::vector<FaultCorridor>& faultCorridors, std::vector<std::vector<size_t>>& faultSeeds) {

            //std::cout << "LevelGeoFaultFactory::generateFaultSeeds() " << std::endl;
            //std::cout << "------------------------------------- " << std::endl << std::endl;

            faultSeeds.resize(faultCorridors.size());

            std::map<int, std::pair<size_t,size_t>> yIDtoFaultCorridor;
            for (size_t i=0; i<faultCorridors.size(); i++){
                FaultCorridor& faultCorridor = faultCorridors[i];
                const std::vector<int> & faultYs = faultCorridor.faultYs();

                faultSeeds[i].resize(faultYs.size());
                for (size_t j=0; j<faultYs.size(); j++) {
                    yIDtoFaultCorridor[faultYs[j]] = std::make_pair(i, j);
                }
            }

            BoundaryIterator<GV> bIt(gridView_, BoundaryIterator<GV>::begin);
            BoundaryIterator<GV> bEnd(gridView_, BoundaryIterator<GV>::end);
            for(; bIt!=bEnd; ++bIt) {

                const Element& insideElement = bIt->inside();
                const auto& geometry = insideElement.geometry();

                // neglect "upper" part of boundary, where (y==1 and x>0) or (y>0 and x==1)
                const auto& center = bIt->geometry().center();
                if (center[0] + center[1] >1)
                        continue;

                const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());

                for (int i=0; i<refElement.size(bIt->indexInInside(), 1, dimElement); i++) {
                    size_t idxInElement = refElement.subEntity(bIt->indexInInside(), 1, i, dimElement);
                    size_t globalIdx = indexSet_.subIndex(insideElement, idxInElement, dimElement);

                    const auto& vertex = geometry.corner(idxInElement);

                    const int yID = computeYID(vertex);

                    if (yIDtoFaultCorridor.count(yID)) {
                        const std::pair<size_t, size_t>& indices = yIDtoFaultCorridor[yID];
                        faultSeeds[indices.first][indices.second] = globalIdx;
                    }
                }
            }

            //std::cout << "------------------------------------- " << std::endl << std::endl;
        }


        void generateFault(std::shared_ptr<FaultInterface<GV>> fault, const size_t faultSeedIdx, const std::set<size_t>& separatingDofs, const std::set<int>& separatingYIDs, const double maxAngle) {
            // compute target at upper/right part of boundary

            //std::cout << "LevelGeoFaultFactory::generateFault() " << std::endl;
            //std::cout << "------------------------------------- " << std::endl << std::endl;

            bool success = false;

            const auto& vertexPosition = vertexPositions_[faultSeedIdx];
            ctype n = vertexPosition[1]-vertexPosition[0];

            Dune::FieldVector<ctype, dimworld> target;
            if (n>0) {
                target[0] = 1-n;
                target[1] = 1;
            } else {
                target[0] = 1;
                target[1] = 1+n;
            }

            std::set<size_t> faultDofs;
            faultDofs.insert(faultSeedIdx);

            std::vector<size_t> faultFaces;

            std::map<size_t, std::vector<size_t>> vertexToAdmissibleFaces;

            std::queue<size_t> vertexQueue;
            vertexQueue.push(faultSeedIdx);

            while (!vertexQueue.empty()) {
                const size_t vertexID = vertexQueue.front();
                vertexQueue.pop();

                if (isTargetBoundaryVertex(vertexPositions_[vertexID])) {
                    success = true;
                    break;
                }

                if (vertexToAdmissibleFaces.count(vertexID)==0) {
                    const std::vector<size_t>& faces = vertexToFaces_[vertexID];
                    std::vector<size_t> admissibleFaces;

                    for (size_t i=0; i<faces.size(); i++) {
                        if (intersectionAllowed(faces[i], vertexID, maxAngle, separatingDofs, separatingYIDs, faultDofs)) {
                            admissibleFaces.push_back(faces[i]);
                        }
                    }

                    std::vector<double> distances(admissibleFaces.size(), 0);
                    for (size_t i=0; i<admissibleFaces.size(); i++) {
                        Dune::FieldVector<ctype, dimworld> vec(target);
                        vec -= faces_[admissibleFaces[i]].geometry().center();
                        distances[i] = vec.two_norm();
                    }
                    distanceSort(admissibleFaces, distances);

                    vertexToAdmissibleFaces[vertexID] = admissibleFaces;
                }

                std::vector<size_t>& suitableFaces = vertexToAdmissibleFaces[vertexID];

                if (suitableFaces.size()==0) {
                    faultDofs.erase(vertexID);
                    faultFaces.pop_back();
                    vertexToAdmissibleFaces.erase(vertexID);
                } else {

                    // generate random number from (0,1)
                    double randomNumber = ((double) std::rand() / (RAND_MAX));

                    size_t randSelector = (size_t) std::abs(std::log2(randomNumber));
                    if (randSelector >= suitableFaces.size()) {
                        randSelector = suitableFaces.size() -1;
                    }

                    size_t nextFaceID = suitableFaces[randSelector];

                    suitableFaces.erase(suitableFaces.begin()+randSelector);

                    std::set<size_t> nextFaceDofs;
                    computeIntersectionDofs(indexSet_, faces_[nextFaceID], nextFaceDofs);
                    nextFaceDofs.erase(vertexID);

                    faultDofs.insert(nextFaceDofs.begin(), nextFaceDofs.end());
                    faultFaces.push_back(nextFaceID);

                    size_t nextVertexID = *(nextFaceDofs.begin());
                    vertexQueue.push(nextVertexID);
                }
            }

            if (!success) {
                std::cout << "Generating a fault failed: Unable to reach target boundary! This should not happend, FaultCorridors should prevent it!" << std::endl;
                DUNE_THROW(Dune::Exception, "Generating a fault failed: Unable to reach target boundary!");
            }

            for (size_t i=0; i<faultFaces.size(); i++) {
                fault->addFace(faces_[faultFaces[i]]);
            }

            //std::cout << "------------------------------------- " << std::endl << std::endl;
        }

    public:
        LevelGeoFaultFactory(const std::shared_ptr<GridType> grid, const int level, const ctype resolution, const bool allowLevelIntersections = false) :
            grid_(grid),
            level_(level),
            resolution_(resolution),
            allowLevelIntersections_(allowLevelIntersections),
            gridView_(grid_->levelGridView(level_)),
            indexSet_(gridView_.indexSet()),
            faceMapper_(gridView_)
            {
                // init faces_, vertexPositions_, vertexToFaces_
                faces_.resize(faceMapper_.size());
                vertexPositions_.resize(gridView_.size(dim));
                vertexToFaces_.resize(gridView_.size(dim));

                std::vector<bool> faceHandled(faceMapper_.size(), false);

                for (const auto& elem:elements(gridView_)) {
                    for (const auto& isect:intersections(gridView_, elem)) {
                        //const size_t faceID = faceMapper_.index(isect);

                        const int faceID = faceMapper_.subIndex(elem, isect.indexInInside(), 1);

                        if (isect.boundary())
                            continue;

                        if (faceHandled[faceID])
                            continue;

                        faceHandled[faceID] = true;
                        faces_[faceID] = isect;

                        const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(elem.type());

                        for (int i=0; i<refElement.size(isect.indexInInside(), 1, dimElement); i++) {
                            //size_t vertexID = idxSet.subIndex(elem.subEntity<1>(isect), i, dimElement);

                            size_t idxInElement = refElement.subEntity(isect.indexInInside(), 1, i, dimElement);
                            const size_t vertexID = indexSet_.subIndex(elem, idxInElement, dimElement);

                            vertexPositions_[vertexID] = elem.geometry().corner(idxInElement);
                            vertexToFaces_[vertexID].push_back(faceID);
                        }
                    }
                }
            }

        void build(std::vector<std::shared_ptr<FaultInterface<GV>>>& faults, std::vector<FaultCorridor>& faultCorridors, std::set<size_t> separatingDofs, double maxAngle = 2) {
            std::vector<std::vector<size_t>> faultSeeds;
            generateFaultSeeds(faultCorridors, faultSeeds);

            faults.resize(0);

            for (size_t i=0; i<faultCorridors.size(); i++) {
                const std::vector<int>& allSeparatingYIDs = faultCorridors[i].separatingYIDs();

                if (faultSeeds[i].size()==0)
                    continue;

                if (faultSeeds[i].size()==1) {
                    std::shared_ptr<FaultInterface<GV>> fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
                    std::set<int> separatingYIDs;

                    generateFault(fault, faultSeeds[i][0], separatingDofs, separatingYIDs, maxAngle);
                    faults.push_back(fault);
                } else {
                    // generate first fault whose limits are given by coarser fault and one separatingYID
                    std::shared_ptr<FaultInterface<GV>> fault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
                    std::set<int> separatingYIDs;
                    separatingYIDs.insert(allSeparatingYIDs[0]);

                    generateFault(fault, faultSeeds[i][0], separatingDofs, separatingYIDs, maxAngle);
                    faults.push_back(fault);

                    // generate all intermediate faults whose limits are given by two separatingYIDs
                    for (size_t j=1; j<faultSeeds[i].size()-1; j++) {
                        std::shared_ptr<FaultInterface<GV>> interface = std::make_shared<FaultInterface<GV>>(gridView_, level_);
                        std::set<int> separatingYIDs;
                        separatingYIDs.insert(allSeparatingYIDs[j-1]);
                        separatingYIDs.insert(allSeparatingYIDs[j]);

                        generateFault(interface, faultSeeds[i][j], separatingDofs, separatingYIDs, maxAngle);
                        faults.push_back(interface);
                    }

                    // generate last fault whose limits are given by coarser fault and one separatingYID
                    std::shared_ptr<FaultInterface<GV>> lastFault = std::make_shared<FaultInterface<GV>>(gridView_, level_);
                    std::set<int> lastSeparatingYIDs;
                    lastSeparatingYIDs.insert(allSeparatingYIDs[allSeparatingYIDs.size()-1]);

                    generateFault(lastFault, faultSeeds[i][faultSeeds[i].size()-1], separatingDofs, lastSeparatingYIDs, maxAngle);
                    faults.push_back(lastFault);
                }
            }
        }
};

template <class GridType>
class GeoFaultFactory
{
    //! Parameter for mapper class
    template<int dim>
    struct FaceMapperLayout
    {
        bool contains (Dune::GeometryType gt)
        {
            return gt.dim() == dim-1;
        }
    };

    protected:
        typedef OscUnitCube<GridType, 2> GridOb;
        static const int dimworld = GridType::dimensionworld ;
        static const int dim = GridType::dimension;
        typedef typename GridType::ctype ctype;

        typedef typename GridType::LevelGridView GV;
        typedef typename GridType::LevelIntersection Intersection;

        typedef typename GridType::template Codim<0>::Entity Element;
        static const int dimElement = Element::dimension;

        typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;


    private:
        const std::vector<int>& faultToLevel_;
        const int coarseResolution_;
        const size_t maxLevel_;
        const bool allowIntersections_;
        const bool allowLevelIntersections_;
        const int coarseGridN_;

        std::vector<std::set<size_t>> separatingDofs_;
        std::vector<int> faultPerLevel_;

        std::shared_ptr<GridType> grid_;
	
        std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_;

        std::vector<std::shared_ptr<LevelGeoFaultFactory<GridType>>> levelGeoFaultFactories_;

        std::vector<std::vector<FaultCorridor>> faultCorridorHierarchy_;

        std::vector<double> levelResolutions_;

private:

typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type
            almost_equal(ctype x, ctype y, int ulp) const
        {
            return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp
                   || std::abs(x-y) < std::numeric_limits<ctype>::min();
        }



bool isSeparatingIntersection(const Intersection& intersection, const std::set<int>& separatingYIDs) const {
    Dune::FieldVector<ctype, dimworld> faceCenter = intersection.geometry().center();

    ctype yID = faceCenter[1]*coarseGridN_;
    ctype intpart;

    if (almost_equal(std::modf(yID, &intpart), 0.0, 2)) {
        std::set<int>::iterator it = separatingYIDs.find((int) intpart);
        if (it!=separatingYIDs.end())
            return true;
        else
            return false;
    } else
        return false;
}


void computeIntersectionDofs(const typename GV::IndexSet& idxSet, const Intersection& intersection, std::set<size_t>& intersectionDofs) {
    intersectionDofs.clear();

    // loop over all vertices of the intersection
    const Element& insideElement = intersection.inside();

    const auto& refElement = Dune::ReferenceElements<double,dimElement>::general(insideElement.type());
    for (int i=0; i<refElement.size(intersection.indexInInside(), 1, dimElement); i++) {
        size_t idxInElement = refElement.subEntity(intersection.indexInInside(), 1, i, dimElement);
        size_t globalIdx = idxSet.subIndex(insideElement, idxInElement, dimElement);

        intersectionDofs.insert(globalIdx);
    }
}

public:
        //setup 
    GeoFaultFactory(const std::vector<int>& faultToLevel, const int coarseResolution, const size_t maxLevel, const bool allowIntersections = false, const bool allowLevelIntersections = false, const double maxAngle = 2) :
        faultToLevel_(faultToLevel),
        coarseResolution_(coarseResolution),
        maxLevel_(maxLevel),
        allowIntersections_(allowIntersections),
        allowLevelIntersections_(allowLevelIntersections),
        coarseGridN_(std::pow(2, coarseResolution_)),
        interfaceNetwork_(nullptr)
        {   
            // set faultPerLevel_:
            faultPerLevel_.resize(maxLevel_+1, 0);
            for (size_t i=0; i<faultToLevel_.size(); i++) {
                faultPerLevel_[faultToLevel_[i]]++;
            }

            print(faultPerLevel_, "faultPerLevel: ");

            Dune::UGGrid<dim>::setDefaultHeapSize(4000);
            GridOb unitCube(coarseGridN_);
            grid_ = unitCube.grid();

            grid_->globalRefine(maxLevel_);
            separatingDofs_.resize(maxLevel_+1);
            levelResolutions_.resize(maxLevel_+1);
            faultCorridorHierarchy_.resize(maxLevel_+1);
            levelGeoFaultFactories_.resize(maxLevel_+1);

            // init interface network
            interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);

            std::vector<FaultCorridor>& coarsefaultCorridors = faultCorridorHierarchy_[0];

            if (faultPerLevel_[0]==0) {
                FaultCorridor newFaultCorridor(0, -coarseGridN_, coarseGridN_, faultToLevel_);
                coarsefaultCorridors.push_back(newFaultCorridor);
            } else {
                int yStep = (2*coarseGridN_-2)/(faultPerLevel_[0]);
                int offSet = (2*coarseGridN_-yStep*faultPerLevel_[0])/2 - coarseGridN_;

                //std::cout << "yStep: " << yStep << " offSet: " << offSet << std::endl;

                FaultCorridor newFaultCorridor(0, offSet, faultPerLevel_[0]*yStep + offSet, faultToLevel_);
                coarsefaultCorridors.push_back(newFaultCorridor);
            }

            for (size_t i=0; i<=maxLevel_; i++) {
                std::vector<std::shared_ptr<FaultInterface<GV>>> faults;
                std::vector<FaultCorridor>& faultCorridors = faultCorridorHierarchy_[i];

                levelResolutions_[i] = std::pow(2, coarseResolution_+i);

                //generate faults on level
                levelGeoFaultFactories_[i] = std::make_shared<LevelGeoFaultFactory<GridType>>(grid_, i, levelResolutions_[i], allowLevelIntersections_);
                levelGeoFaultFactories_[i]->build(faults, faultCorridors, separatingDofs_[i], maxAngle);

                for (size_t j=0; j<faults.size(); j++) {
                    interfaceNetwork_->addInterface(faults[j]);
                }

                if (i<maxLevel_) {
                    const int newLevel = i+1;

                    //compute separating dofs for next level given by current levelinterfacenetwork
                    interfaceNetwork_->prolongLevel(i, newLevel);

                    if (!allowIntersections_) {
                        separatingDofs_[newLevel] = interfaceNetwork_->getInterfaceNetworkDofs(newLevel);
                    }

                    //compute faultCorridors for next level
                    std::vector<FaultCorridor>& newLevelFaultCorridors = faultCorridorHierarchy_[newLevel];
                    newLevelFaultCorridors.resize(0);

                    for (size_t j=0; j<faultCorridors.size(); j++) {
                        faultCorridors[j].nextLevel(newLevelFaultCorridors);
                    }
                }
            }
        }

    /*
    void prolongToAll() {
        // prolong all faults to all subsequent levels
        for (int i=maxLevel_-1; i>=0; i--) {
            if (interfaceNetwork_->size(i)>0) {
                std::set<int> toLevels;
                for (size_t j=i+1; j<=maxLevel_; j++) {
                    toLevels.insert(j);
                }

                interfaceNetwork_->prolongLevelInterfaces(i, toLevels);
            }
        }
        interfaceNetwork_->build();
    }*/

    /*void prolongToAll() {
        // prolong all faults to all subsequent levels
        for (int i=interfaceNetwork_->size()-1; i>=0; i--) {
            interfaceNetwork_->prolongLevelInterfaces(i, maxLevel_);
        }
        interfaceNetwork_->build();
    }*/

    const GridType& grid() const {
        return *grid_;
	}
	
    /*const InterfaceNetwork<GridType>& interfaceNetwork() {
        return *interfaceNetwork_;
    }*/

    InterfaceNetwork<GridType>& interfaceNetwork() {
        return *interfaceNetwork_;
    }
};
#endif