#ifndef DEBUG_UTILS_
#define DEBUG_UTILS_

#include <string>

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

#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>

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

//using namespace std;


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

    template <int dim, typename ctype=double>
    void print(const Dune::FieldVector<ctype, dim>& x, std::string message){
       if (message != "")
        std::cout << message << std::endl;

       for (size_t i=0; i<x.size(); i++) {
           std::cout << x[i] << ";" << std::endl;
       }

   }

    template <int dim, typename ctype=double>
    void print(const Dune::BlockVector<Dune::FieldVector<ctype, dim>>& x, std::string message){
       std::cout << message  << " " << "size " << x.size() << std::endl;
       for (size_t i=0; i<x.size(); i++) {
           print(x[i], "");
       }
       std::cout << std::endl << std::endl;
   }

   template <int dim, typename ctype=double>
   void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, bool endOfLine = true){
        if (message != "")
            std::cout << message << std::endl;

        for (size_t i=0; i<mat.N(); i++) {
            for (size_t j=0; j<mat.M(); j++) {
                if (mat.exists(i,j))
                    std::cout << mat[i][j] << " ";
                else
                    std::cout << "0 ";
            }

            if (endOfLine)
                std::cout << ";"<< std::endl;
        }
   }

    template <int dim, typename ctype=double>
    void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, const size_t line, bool endOfLine = true){
        if (message != "")
            std::cout << message << std::endl;

        assert(line<mat.N());
        for (size_t j=0; j<mat.M(); j++) {
            if (mat.exists(line,j))
                    std::cout << mat[line][j] << " ";
                else
                    std::cout << "0 ";
            }

            if (endOfLine)
                std::cout << ";"<< std::endl;
   }

    template <int dim, typename ctype=double>
    void print(const Dune::BCRSMatrix<Dune::FieldMatrix<ctype, dim, dim>>& mat, std::string message){
        std::cout << message << " " << "size (" << mat.N() << ", " << mat.M() << ")" <<std::endl;

        Dune::FieldMatrix<ctype, dim, dim> zeroEntry = 0;

        for (size_t i=0; i<mat.N(); i++) {
            for (size_t d=0; d<dim; d++) {
                for (size_t j=0; j<mat.M(); j++) {
                    if (mat.exists(i,j))
                        print(mat[i][j], "", d, j==mat.M()-1);
                    else
                        print(zeroEntry, "", d, j==mat.M()-1);
                }
            }
       }
       std::cout << std::endl;
   }

   template <class T=Dune::FieldMatrix<double,1,1>>
   void print(const Dune::Matrix<T>& mat, std::string message){
       std::cout << message << std::endl;
       for (size_t i=0; i<mat.N(); i++) {
           for (size_t j=0; j<mat.M(); j++) {
                   std::cout << mat[i][j] << " ";
           }
           std::cout << std::endl;
       }
       std::cout << std::endl;
   }

   template <class T>
   void print(const std::vector<T>& x, std::string message){
      std::cout << message << std::endl;
      for (size_t i=0; i<x.size(); i++) {
          std::cout << x[i] << std::endl;
      }
      std::cout << std::endl << std::endl;
   }

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

   void print(const std::set<size_t>& x, std::string message){
      std::cout << message << std::endl;
      std::set<size_t>::iterator dofIt = x.begin();
      std::set<size_t>::iterator dofEndIt = x.end();
      for (; dofIt != dofEndIt; ++dofIt) {
          std::cout << *dofIt << " ";
      }
      std::cout << std::endl << std::endl;
   }

   void print(const std::set<int>& x, std::string message){
      std::cout << message << std::endl;
      std::set<int>::iterator dofIt = x.begin();
      std::set<int>::iterator dofEndIt = x.end();
      for (; dofIt != dofEndIt; ++dofIt) {
          std::cout << *dofIt << " ";
      }
      std::cout << std::endl << std::endl;
   }

   void step(const double stepNumber, std::string message=""){
      std::cout << message << " Step " << stepNumber << "!" << std::endl;
   }*/

   template <class GridView, class VectorType>
   void writeToVTK(const GridView& gridView, const VectorType& x, const std::string path, const std::string name) {
       Dune::VTKWriter<GridView> vtkwriter(gridView);

       std::ofstream lStream( "garbage.txt" );
       std::streambuf* lBufferOld = std::cout.rdbuf();
       std::cout.rdbuf(lStream.rdbuf());

       vtkwriter.addVertexData(x,"data");
       vtkwriter.pwrite(name, path, path);

       std::cout.rdbuf( lBufferOld );
   }

   template <class GridView>
   void writeToVTK(const GridView& gridView, const std::string path, const std::string name) {
       Dune::VTKWriter<GridView> vtkwriter(gridView);

       std::ofstream lStream( "garbage.txt" );
       std::streambuf* lBufferOld = std::cout.rdbuf();
       std::cout.rdbuf(lStream.rdbuf());

       vtkwriter.pwrite(name, path, path);

       std::cout.rdbuf( lBufferOld );
   }

   /*
   template <class VectorType, class DGBasisType>
   void writeToVTK(const DGBasisType& basis, const VectorType& x, const std::string path, const std::string name) {
       Dune::VTKWriter<typename DGBasisType::GridView> vtkwriter(basis.getGridView());
       VectorType toBePlotted(x);

       toBePlotted.resize(basis.getGridView().indexSet().size(DGBasisType::GridView::Grid::dimension));

       std::ofstream lStream( "garbage.txt" );
       std::streambuf* lBufferOld = std::cout.rdbuf();
       std::cout.rdbuf( lStream.rdbuf() );

       vtkwriter.addVertexData(toBePlotted,"data");
       vtkwriter.pwrite(name, path, path);

       std::cout.rdbuf( lBufferOld );
   }*/

   template <class Intersection, class Basis>
   void printIntersection(const Intersection& it, const Basis& basis, std::string message = "") {
       if (message != "") {
           std::cout << message << std::endl;
       }

       const auto& inside = it.inside();
       const auto& localCoefficients = basis.getLocalFiniteElement(inside).localCoefficients();

       std::cout << "dofs: ";
       for (size_t i=0; i<localCoefficients.size(); i++) {
           unsigned int entity = localCoefficients.localKey(i).subEntity();
           unsigned int codim  = localCoefficients.localKey(i).codim();

           if (it.containsInsideSubentity(entity, codim)) {
               int dofIndex = basis.index(it.inside(), i);
               std::cout << dofIndex << " ";
               break;
           }
       }
       std::cout << std::endl;
   }

   template <class BoundaryPatch, class Basis>
   void printBoundary(const BoundaryPatch& patch, const Basis& basis, std::string message = "") {
       if (message != "") {
           std::cout << message << std::endl;
       }

       for (auto bIt = patch.begin(); bIt != patch.end(); ++bIt) {
            printIntersection(*bIt, basis, "");
       }

       std::cout << std::endl;
   }

   template <class BasisType, typename ctype=double>
   void printBasisDofLocation(const BasisType& basis) {
     /*  typedef typename BasisType::GridView GridView;

       const int dim = GridView::dimension;

       std::map<int, int> indexTransformation;
       std::map<int, Dune::FieldVector<ctype, dim>> indexCoords;


       const GridView& gridView = basis.getGridView();
       const int gridN = std::pow(gridView.indexSet().size(dim), 1.0/dim)-1;

       for (auto& it : elements(gridView)) {
           const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
           const auto& geometry = it.geometry();

           for(int i=0; i<geometry.corners(); ++i) {
               const auto& vertex = geometry.corner(i);
               const auto& local = geometry.local(vertex);

               int vertexIndex = vertex[0]*gridN;
               for (int j=1; j<dim; ++j){
                   vertexIndex += vertex[j]*gridN*std::pow(gridN+1, j);
               }

               const int localBasisSize = tFE.localBasis().size();
               std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
               tFE.localBasis().evaluateFunction(local, localBasisRep);

               for(int k=0; k<localBasisSize; ++k) {
                   if (localBasisRep[k]==1) {
                       int dofIndex = basis.index(it, k);

                       if (indexTransformation.count(dofIndex)==0) {
                           indexTransformation[dofIndex] = vertexIndex;
                           indexCoords[dofIndex] = vertex;
                       }

                       break;
                   }
               }
           }
       }

       std::cout << std::endl;
       std::cout << "Coarse level: Geometric vertex indices vs. associated basis dof indices " << indexTransformation.size() << std::endl;
       std::map<int,int>::iterator mapIt = indexTransformation.begin();
       std::map<int,int>::iterator mapEndIt = indexTransformation.end();
       for (; mapIt!=mapEndIt; ++mapIt) {
           std::cout << mapIt->second << " => " << mapIt->first << " Coords: ";

           const Dune::FieldVector<ctype, dim>& coords = indexCoords[mapIt->first];
           for (size_t i=0; i<coords.size(); i++) {
               std::cout << coords[i] << " ";
           }
           std::cout << std::endl;
       }
       std::cout << std::endl;*/

       const int dim = BasisType::GridView::dimension;

       std::map<int, Dune::FieldVector<ctype, dim>> indexCoords;


       const auto& gridView = basis.getGridView();

       for (auto& it : elements(gridView)) {
            const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
            const auto& geometry = it.geometry();

            for(int i=0; i<geometry.corners(); ++i) {
                const auto& vertex = geometry.corner(i);
                const auto& local = geometry.local(vertex);

                const int localBasisSize = tFE.localBasis().size();
                std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
                tFE.localBasis().evaluateFunction(local, localBasisRep);

                for(int k=0; k<localBasisSize; ++k) {
                    if (1.0 - localBasisRep[k]< 1e-14) {
                        int dofIndex = basis.index(it, k);
                        indexCoords[dofIndex] = vertex;
                        break;
                    }
                }
              }
       }

       std::cout << "Coords of basis dofs: " << indexCoords.size() << std::endl;
       auto&& mapIt = indexCoords.begin();
       const auto& mapEndIt = indexCoords.end();
       for (; mapIt!=mapEndIt; ++mapIt) {
            std::cout << mapIt->first << " Coords: ";

            const auto& coords = mapIt->second;
            for (size_t i=0; i<coords.size(); i++) {
                std::cout << coords[i] << " ";
            }
            std::cout << std::endl;
       }

       std::cout << std::endl;
   }

   template <class GridView>
   void printDofLocation(const GridView& gridView) {

       std::cout << "-- GridView vertices --" << std::endl;
       std::cout << "Index   Coords         " << std::endl;
       std::cout << "-----------------------" << std::endl;

       Dune::MultipleCodimMultipleGeomTypeMapper<
           GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());

       for (auto&& it : vertices(gridView)) {
            const auto i = vertexMapper.index(it);
            const auto& geometry = it.geometry();
            const auto& corner = geometry.corner(0);

            std::cout << std::setw(5) << i << "   ";
            for(size_t i=0; i<corner.size(); ++i) {
                std::cout << std::setw(5) << std::setprecision(3) << corner[i] << " ";
            }
            std::cout << std::endl;
       }
       std::cout << std::endl;
   }

   template <class Vector, class NBodyAssembler>
   void printMortarBasis(const NBodyAssembler& nBodyAssembler) {
       std::cout << "-- Mortar basis in canonical coords --" << std::endl;
       std::cout << "--------------------------------------" << std::endl;

       const auto& dim = nBodyAssembler.getTransformationMatrix().N();

       Vector mortarBasisVector(dim);
       std::vector<Vector> canonicalBasisVector(nBodyAssembler.nGrids());
       for (size_t i=0; i<dim; i++) {
            mortarBasisVector = 0;
            mortarBasisVector[i] = 1;

            nBodyAssembler.postprocess(mortarBasisVector, canonicalBasisVector);
            print(canonicalBasisVector, "canonicalBasisVector " + std::to_string(i));
       }
       std::cout << std::endl;
   }
#endif