#ifndef SUPPORT_PATCH_FACTORY_HH
#define SUPPORT_PATCH_FACTORY_HH

#include <queue>
#include <set>

#include <dune/common/bitsetvector.hh>
#include <dune/common/fvector.hh>

#include <dune/fufem/referenceelementhelper.hh>
#include <dune/grid/common/mcmgmapper.hh>

#include <dune/faultnetworks/hierarchicleveliterator.hh>

template <class BasisType>
class SupportPatchFactory
{
    protected:
        typedef typename BasisType::GridView GV;
        typedef typename BasisType::LocalFiniteElement LFE;
	
	typedef typename GV::Grid GridType;
	typedef typename GridType::ctype ctype;
	static const int dim = GridType::dimension;
	typedef typename GridType::template Codim<0>::Entity Element;
	typedef typename GridType::template Codim<dim>::Entity Vertex;

	typedef HierarchicLevelIterator<GridType> HierarchicLevelIteratorType;
	typedef typename Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,  Dune::MCMGElementLayout > ElementMapper;
	
    private:
	const BasisType& coarseBasis_;
	const BasisType& fineBasis_;
	
	const int coarseLevel_;
	const int fineLevel_;
	
 	const std::vector< std::vector <Element>>& vertexInElements_;
	const int centerVertexIdx_;
	const int patchDepth_;

  
	const GridType& grid;
	const GV& coarseGridView;
	const GV& fineGridView;
	
	std::vector<int> localToGlobal;
	Dune::BitSetVector<1> boundaryDofs;
	std::vector<Element> fineRegionElements;
	
	ElementMapper coarseMapper;

    public:
        //setup 
    SupportPatchFactory(const BasisType& coarseBasis, const BasisType& fineBasis, const int coarseLevel, const int fineLevel, const std::vector< std::vector <Element>>& vertexInElements, const int centerVertexIdx, const int patchDepth) :
            coarseBasis_(coarseBasis),
            fineBasis_(fineBasis),
            coarseLevel_(coarseLevel),
            fineLevel_(fineLevel),
            vertexInElements_(vertexInElements),
            centerVertexIdx_(centerVertexIdx),
            patchDepth_(patchDepth),
            grid(coarseBasis_.getGridView().grid()),
	    coarseGridView(coarseBasis_.getGridView()),
	    fineGridView(fineBasis_.getGridView()),
	    coarseMapper(grid, coarseLevel_)
        {   
	    fineRegionElements.resize(0);
	  
            std::vector<Element> coarseElements;
	    Dune::BitSetVector<1> visited(coarseMapper.size());
	    Dune::BitSetVector<1> vertexVisited(vertexInElements_.size());
            visited.unsetAll();
	    vertexVisited.unsetAll();
	    
	    std::set<int> localDofs;
	    std::set<int> localBoundaryDofs;

	    addVertexSupport(centerVertexIdx_, coarseElements, visited, vertexVisited);
	    
	    // construct coarse patch
	    for (int depth=1; depth<patchDepth_; depth++) {
		std::vector<Element> coarseElementNeighbors;
		
		for (size_t i=0; i<coarseElements.size(); ++i) {
		    const Element& elem = coarseElements[i];     
		    const LFE& coarseLFE = coarseBasis_.getLocalFiniteElement(elem);
	      
		    for (size_t j=0; j<coarseLFE.localBasis().size(); ++j) {
            int dofIndex = coarseBasis_.indexInGridView(elem, j);
			
			if (!vertexVisited[dofIndex][0]) {
			    addVertexSupport(dofIndex, coarseElementNeighbors, visited, vertexVisited);
			}
		    }
		}	        
		    
		coarseElements.insert(coarseElements.begin(), coarseElementNeighbors.begin(), coarseElementNeighbors.end());
	    }
	    

	    // construct fine patch
	    for (size_t i=0; i<coarseElements.size(); ++i) {
		const Element& elem = coarseElements[i];
		addFinePatchElements(elem, visited, localDofs, localBoundaryDofs);
	    }

	 
	 
        localToGlobal.resize(localDofs.size());
	    boundaryDofs.resize(localDofs.size());
	    boundaryDofs.unsetAll();
	    
	    std::set<int>::iterator dofIt = localDofs.begin();
	    std::set<int>::iterator dofEndIt = localDofs.end();
	    size_t i=0;
	    for (; dofIt != dofEndIt; ++dofIt) {
		int localDof = *dofIt;
		localToGlobal[i] = localDof;
		
		if (localBoundaryDofs.count(localDof)) {
		    boundaryDofs[i][0] = true;
		}
		i++;
	    }
	}

	std::vector<int>& getLocalToGlobal() {
	    return localToGlobal;
	}
	
	std::vector<Element>& getRegionElements() {
	    return fineRegionElements;
	}
	
	Dune::BitSetVector<1>& getBoundaryDofs() {
	    return boundaryDofs;
	}
	
	
private:


    void print(const Dune::BitSetVector<1>& x, std::string message){
       std::cout << message << std::endl;
       for (size_t i=0; i<x.size(); i++) {
           std::cout << x[i][0] << " ";
       }
       std::cout << std::endl << std::endl;
   }

	void addVertexSupport(const int vertexIdx, std::vector<Element>& coarseElements, Dune::BitSetVector<1>& visited, Dune::BitSetVector<1>& vertexVisited) {
	    const std::vector<Element>& vertexInElement = vertexInElements_[vertexIdx];
	    
	    for (size_t i=0; i<vertexInElement.size(); i++) {
		const Element& elem = vertexInElement[i];
		const int elemIdx = coarseMapper.index(elem);
		
		if (!visited[elemIdx][0]) {
		    coarseElements.push_back(elem);
		    visited[elemIdx][0] = true;
		}
	    }
	    vertexVisited[vertexIdx][0] = true;
	}
      
	bool containsInsideSubentity(const Element& elem, const typename GridType::LevelIntersection& intersection, int subEntity, int codim) {
	    return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
	}
      
      
	void addFinePatchElements(const Element& coarseElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
	    HierarchicLevelIteratorType endHierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::end, fineLevel_);
	    for (HierarchicLevelIteratorType hierIt(grid, coarseElem, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel_); hierIt!=endHierIt; ++hierIt) {

		const Element& fineElem = *hierIt;
		fineRegionElements.push_back(fineElem);

		addLocalDofs(coarseElem, fineElem, inCoarsePatch, localDofs, localBoundaryDofs);
	    }
	}	  
	
	void addLocalDofs(const Element& coarseElem, const Element& fineElem, const Dune::BitSetVector<1>& inCoarsePatch, std::set<int>& localDofs, std::set<int>& localBoundaryDofs) {
	    const LFE& fineLFE = fineBasis_.getLocalFiniteElement(fineElem);
	    
        /*
	    const auto& fineGeometry = fineElem.geometry();
	    const auto& coarseGeometry = coarseElem.geometry();
		
        const auto& ref = Dune::ReferenceElements<ctype, dim>::general(fineElem.type()); */
	    
	    // insert local dofs
	    for (size_t i=0; i<fineLFE.localBasis().size(); ++i) {
		int dofIndex = fineBasis_.index(fineElem, i);
		localDofs.insert(dofIndex);
	    }

	    
	    // search for parts of the global and inner boundary 
	    typename GridType::LevelIntersectionIterator neighborIt = fineGridView.ibegin(fineElem);
	    typename GridType::LevelIntersectionIterator neighborItEnd = fineGridView.iend(fineElem);
	    for (; neighborIt!=neighborItEnd; ++neighborIt) {
		const typename GridType::LevelIntersection& intersection = *neighborIt;

		bool isInnerBoundary = false;
		if (intersection.neighbor()) {
		    Element outsideFather = intersection.outside();
		    
		    while (outsideFather.level()>coarseLevel_) {
			outsideFather = outsideFather.father();
		    }
		    
		    if (!inCoarsePatch[coarseMapper.index(outsideFather)][0]) {
			isInnerBoundary = true;
		    }
		}
		
		if (intersection.boundary() or isInnerBoundary) {
		    typedef typename LFE::Traits::LocalCoefficientsType LocalCoefficients;
		    const LocalCoefficients* localCoefficients = &fineBasis_.getLocalFiniteElement(intersection.inside()).localCoefficients();

		    for (size_t i=0; i<localCoefficients->size(); i++) {
			unsigned int entity = localCoefficients->localKey(i).subEntity();
			unsigned int codim  = localCoefficients->localKey(i).codim();

			if (containsInsideSubentity(fineElem, intersection, entity, codim))
			    localBoundaryDofs.insert(fineBasis_.index(intersection.inside(), i));
		    } 
		}
	    }	
	}
};
#endif