#ifndef CELL_PATCH_FACTORY_HH
#define CELL_PATCH_FACTORY_HH

#include <set>
#include <queue>

#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 CellPatchFactory
{
    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 Dune::LevelMultipleCodimMultipleGeomTypeMapper<GridType,  Dune::MCMGElementLayout > ElementMapper;
	
    private:
    const BasisType& basis_;
    const int level_;

	const GridType& grid;
    const GV& gridView;
	
    std::vector<std::vector<int>> cellLocalToGlobal;
    std::vector<Dune::BitSetVector<1>> cellBoundaryDofs;
    std::vector<std::vector<Element>> cellElements;
	
    ElementMapper mapper;

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;
    }*/

    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);
    }

public:
        //setup 
    CellPatchFactory(const BasisType& basis, const int level) :
            basis_(basis),
            level_(level),
            grid(basis_.getGridView().grid()),
            gridView(basis_.getGridView()),
            mapper(grid, level_) {

        const auto& interfaces = basis_.faultNetwork();

        cellElements.resize(0);
	  
        Dune::BitSetVector<1> visited(mapper.size());
        visited.unsetAll();

        std::queue<Element> cellSeeds;
        cellSeeds.push(*gridView. template begin <0>());
        size_t cellCount = 0;

        while (!cellSeeds.empty()) {
            const auto cellSeed = cellSeeds.front();
            cellSeeds.pop();

            if (visited[mapper.index(cellSeed)][0] == true)
                continue;

            cellCount++;

            cellElements.resize(cellCount);
            auto& elements = cellElements.back();

            std::queue<Element> elementQueue;
            elementQueue.push(cellSeed);


            cellLocalToGlobal.resize(cellCount);
            cellBoundaryDofs.resize(cellCount);
            std::set<int> dofs;
            std::set<int> boundaryDofs;

            while (!elementQueue.empty()) {
                const auto elem = elementQueue.front();
                const auto elemID = mapper.index(elem);
                elementQueue.pop();

                if (visited[elemID][0] == true)
                    continue;

                elements.push_back(elem);
                visited[elemID][0] = true;

                const LFE& lfe = basis_.getLocalFiniteElement(elem);
                // insert dofs
                for (size_t i=0; i<lfe.localBasis().size(); ++i) {
                    int dofIndex = basis_.index(elem, i);
                    dofs.insert(dofIndex);
                }

                for (const auto& is : intersections(gridView, elem)) {
                    if (is.neighbor()) {
                        const auto neighbor = is.outside();
                        const auto neighborId = mapper.index(neighbor);

                        if (!(visited[neighborId][0])) {
                            if (interfaces.isInterfaceIntersection(is)) {
                                cellSeeds.push(neighbor);
                            } else {
                                elementQueue.push(neighbor);
                            }
                        }
                    }

                    // set boundaryDofs
                    if (is.boundary()) {
                        const auto& localCoefficients = lfe.localCoefficients();

                        for (size_t i=0; i<localCoefficients.size(); i++) {
                            auto entity = localCoefficients.localKey(i).subEntity();
                            auto codim  = localCoefficients.localKey(i).codim();

                            if (containsInsideSubentity(elem, is, entity, codim))
                                boundaryDofs.insert(basis_.index(elem, i));
                        }
                    }
                }
            }

            auto& localToGlobal = cellLocalToGlobal.back();
            auto& localBoundaryDofs = cellBoundaryDofs.back();

            localToGlobal.resize(dofs.size());
            localBoundaryDofs.resize(dofs.size());
            localBoundaryDofs.unsetAll();

            size_t i=0;
            for (const auto& dof : dofs) {
                localToGlobal[i] = dof;

                if (boundaryDofs.count(dof)) {
                    localBoundaryDofs[i][0] = true;
                }
                i++;
            }
        }
	}

    auto& getLocalToGlobal() {
        return cellLocalToGlobal;
	}
	
    auto& getRegionElements() {
        return cellElements;
	}
	
    auto& getBoundaryDofs() {
        return cellBoundaryDofs;
	}
};
#endif