#ifndef DUNE_MMDG_DG_HH
#define DUNE_MMDG_DG_HH

#include <cmath>

#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>

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

#include <dune/geometry/quadraturerules.hh>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/solver.hh>
#include <dune/istl/umfpack.hh>

#include <dune/mmdg/nonconformingp1vtkfunction.hh>

template<class GridView, class Mapper, class Problem>
class DG
{
public:
  static constexpr int dim = GridView::dimension;
  static constexpr auto simplex = Dune::GeometryTypes::simplex(dim);
  static constexpr auto edge = Dune::GeometryTypes::simplex(dim-1);

  using Scalar = typename GridView::ctype;
  // using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
  // using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;

  using Matrix = Dune::DynamicMatrix<Scalar>;
  using Vector = Dune::DynamicVector<Scalar>;

  using VTKFunction = Dune::VTKFunction<GridView>;
  using P1Function = Dune::NonconformingP1VTKFunction<GridView,
    Dune::DynamicVector<Scalar>>;
  using QRuleElem = Dune::QuadratureRule<Scalar, dim>;
  using QRulesElem = Dune::QuadratureRules<Scalar, dim>;
  using QRuleEdge = Dune::QuadratureRule<Scalar, dim-1>;
  using QRulesEdge = Dune::QuadratureRules<Scalar, dim-1>;

  //constructor
  DG (const GridView& gridView, const Mapper& mapper, const Problem& problem) :
    gridView_(gridView), mapper_(mapper), problem_(problem),
    dof_((1 + dim) * gridView.size(0))
    {

      A = std::make_shared<Matrix>(dof_, dof_, 0.0); //stiffness matrix
      b = Vector(dof_, 0.0); //load vector
      d = Vector(dof_, 0.0); //solution vector

      //initialize stiffness matrix A, load vector b and solution vector d
      // A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
      // b = Vector(dof_);
      // d = Vector(dof_);
    }

  void operator() (const Scalar mu)
  {
    //assemble stiffness matrix A and load vector b
    assembleSLE(mu);
    //A->compress();

    (*A).solve(d, b);

    for (int i = 0; i < dof_; i++)
      for (int j = 0; j < i; j++)
        /*assert*/if(std::abs((*A)[i][j] - (*A)[j][i]) >
          std::numeric_limits<Scalar>::epsilon())
            std::cout << i << ", " << j << std::endl;

    // Dune::InverseOperatorResult result;
    // Dune::UMFPack<Matrix> solver(*A);
    // solver.apply(d, b, result);

    //write solution to a vtk file
    writeVTKoutput();
  }

  //compute L2 error using a quadrature rule
  Scalar computeL2error(const int quadratureOrder = 5) const
  {
    if (!problem_.hasExactSolution())
    {
      DUNE_THROW(Dune::Exception, "ERROR: Problem has no exact solution or "
        "exact solution is not implemented");
    }

    Scalar error(0.0);

    for (const auto& elem : elements(gridView_))
    {
      const int elemIdxSLE = (dim + 1) * mapper_.index(elem);
      const auto& geo = elem.geometry();
      const QRuleElem& qrule = QRulesElem::rule(simplex, quadratureOrder);

      //determine L2-error on element elem using a quadrature rule
      for (const auto& ip : qrule)
      {
        const auto& qpLocal = ip.position();
        const auto& qpGlobal = geo.global(qpLocal);

        const Scalar quadratureFactor =
          ip.weight() * geo.integrationElement(qpLocal);

        Scalar numericalSolutionAtQP = d[elemIdxSLE];

        for (int i = 0; i < dim; i++)
        {
          numericalSolutionAtQP += d[elemIdxSLE + i + 1] * qpGlobal[i];
        }

        const Scalar errorEvaluation =
          numericalSolutionAtQP - problem_.exactSolution(qpGlobal);

        error += quadratureFactor * errorEvaluation * errorEvaluation;
      }
    }

    return std::sqrt(error);
  }

  const GridView& gridView_;
  const Mapper& mapper_; //element mapper for gridView_
  const Problem& problem_; //the DG problem to be solved

  const int dof_; //degrees of freedom

protected:
  DG (const GridView& gridView, const Mapper& mapper, const Problem& problem,
    const int dof) :
    gridView_(gridView), mapper_(mapper), problem_(problem), dof_(dof)
    {
      A = std::make_shared<Matrix>(dof_, dof_, 0.0); //stiffness matrix
      b = Vector(dof_, 0.0); //load vector
      d = Vector(dof_, 0.0); //solution vector

      //initialize stiffness matrix A, load vector b and solution vector d
      // A = std::make_shared<Matrix>(dof_, dof_, 10, 1, Matrix::implicit); //TODO
      // b = Vector(dof_);
      // d = Vector(dof_);
    }

  //assemble stiffness matrix A and load vector b
  void assembleSLE (const Scalar mu)
  {
    //quadrature rules
    const QRuleElem& rule_q =
      QRulesElem::rule(simplex, problem_.quadratureOrder_q());
    const QRuleElem& rule_K_Elem =
      QRulesElem::rule(simplex, problem_.quadratureOrder_K());
    const QRuleEdge& rule_K_Edge =
      QRulesEdge::rule(edge, problem_.quadratureOrder_K());
    const QRuleEdge& rule_Kplus1_Edge =
      QRulesEdge::rule(edge, problem_.quadratureOrder_K() + 1);
    const QRuleEdge& rule_2ndOrder = QRulesEdge::rule(edge, 2);
    const QRuleEdge& rule_boundary =
      QRulesEdge::rule(edge, problem_.quadratureOrderBoundary());

    //we use the basis
    // phi_elem,0 (x) = indicator(elem);
    // phi_elem,i (x) = x[i]*indicator(elem);
    //for elem in elements(gridView) and i = 1,...,dim
    for (const auto& elem : elements(gridView_))
    {
      const int elemIdx = mapper_.index(elem);
      const auto& geo = elem.geometry();

      //in the system of linear equations (SLE) Ad = b,
      //the index elemIdxSLE refers to the basis function phi_elem,0
      //and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
      //basis function phi_elem,i
      const int elemIdxSLE = (dim + 1)*elemIdx;

      //fill load vector b using a quadrature rule
      for (const auto& ip : rule_q)
      {
        //quadrature point in local and global coordinates
        const auto& qpLocal = ip.position();
        const auto& qpGlobal = geo.global(qpLocal);

        const Scalar quadratureFactor =
          ip.weight() * geo.integrationElement(qpLocal);
        const Scalar qEvaluation(problem_.q(qpGlobal));

        //quadrature for int_elem q*phi_elem,0 dV
        b[elemIdxSLE] += quadratureFactor * qEvaluation;

        //quadrature for int_elem q*phi_elem,i dV
        for (int i = 0; i < dim; i++)
        {
          b[elemIdxSLE + i + 1] += quadratureFactor * qpGlobal[i] * qEvaluation;
        }
      }

      //evaluate
      // int_elem K*grad(phi_elem,i)*grad(phi_elem,j) dV
      // = int_elem K[j][i] dV
      //using a quadrature rule
      for (const auto& ip : rule_K_Elem)
      {
        //quadrature point in local and global coordinates
        const auto& qpLocal = ip.position();
        const auto& qpGlobal = geo.global(qpLocal);

        const Scalar quadratureFactor =
          ip.weight() * geo.integrationElement(qpLocal);

        for (int i = 0; i < dim; i++)
        {
          for (int j = 0; j <= i; j++)
          {
            (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
              quadratureFactor * problem_.K(qpGlobal)[j][i];
          }
        }
      }

      //iterate over all intersection with the boundary of elem
      for (const auto& intersct : intersections(gridView_, elem))
      {
        if constexpr (dim != 1)
        {
          if (gridView_.grid().isInterface(intersct))
          {
            continue;
          }
        }

        const auto& normal = intersct.centerUnitOuterNormal();
        const auto& intersctGeo = intersct.geometry();
        const double intersctVol = intersctGeo.volume();
        const auto& intersctCenter = intersctGeo.center();

        //create storage for multiply used integral values
        //linearIntegrals[i] = int_intersct x[i] ds
        Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);

        //quadraticIntregrals[i][j] = int_intersct x[i]*x[j] ds
        //NOTE: is there a better type for a symmetric matrix?
        //note that only the lower diagonal and diagonal entries of
        //quadraticIntregrals are used
        Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);

        //KIntegrals[i] = int_intersct K[:,i]*normal ds
        Dune::FieldVector<Scalar, dim> KIntegrals(0.0);

        //KIntegralsX[i][j] = int_intersct x[j]*(K[:,i]*normal) ds
        Dune::FieldMatrix<Scalar, dim, dim> KIntegralsX(0.0);

        //fill vector linearIntegrals and matrix quadraticIntregrals
        for (const auto& ip : rule_2ndOrder)
        {
          //quadrature point in local and global coordinates
          const auto& qpLocal = ip.position();
          const auto& qpGlobal = intersctGeo.global(qpLocal);

          const Scalar quadratureFactor =
            ip.weight() * intersctGeo.integrationElement(qpLocal);

          for (int i = 0; i < dim; i++)
          {
            //use midpoint rule for exact evaluation of
            // int_intersct x[i] ds //TODO
            linearIntegrals[i] = intersctVol * intersctCenter[i];

            for (int j = 0; j <= i; j++)
            {
              quadraticIntregrals[i][j] +=
                quadratureFactor * qpGlobal[i] * qpGlobal[j];
              }
            }
          }

          //quadraticIntregrals is symmetric
          for (int i = 0; i < dim; i++)
          {
            for (int j = 0; j < i; j++)
            {
              quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
            }
          }

        //fill vector KIntegrals using a quadrature rule
        for (const auto& ip : rule_K_Edge)
        {
          //quadrature point in local and global coordinates
          const auto& qpLocal = ip.position();
          const auto& qpGlobal = intersctGeo.global(qpLocal);

          const Scalar quadratureFactor =
            ip.weight() * intersctGeo.integrationElement(qpLocal);
          const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
            problem_.K(qpGlobal);

          for (int i = 0; i < dim; i++)
          {
            KIntegrals[i] += quadratureFactor * (KEvaluation[i] * normal);
          }
        }

        //fill vector KIntegralsX using a quadrature rule
        for (const auto& ip : rule_Kplus1_Edge)
        {
          //quadrature point in local and global coordinates
          const auto& qpLocal = ip.position();
          const auto& qpGlobal = intersctGeo.global(qpLocal);

          const Scalar quadratureFactor =
            ip.weight() * intersctGeo.integrationElement(qpLocal);
          const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
            problem_.K(qpGlobal);

          for (int i = 0; i < dim; i++)
          {
            for (int j = 0; j < dim; j++)
            {
              KIntegralsX[i][j] += quadratureFactor *
                qpGlobal[j] * (KEvaluation[i] * normal);
            }
          }
        }

        //exact evaluation of
        // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
        (*A)[elemIdxSLE][elemIdxSLE] += mu * intersctVol;

        if (intersct.neighbor()) //intersct has neighboring element
        {
          //index of the neighboring element
          const int neighborIdx = mapper_.index(intersct.outside());
          const int neighborIdxSLE = (dim + 1) * neighborIdx;

          for (int i = 0; i < dim; i++)
          { //we use the relations
            // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
            // = mu * int_intersct x[i] ds
            //and
            // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
            // = 0.5 * int_intersct K[:,i] * normal ds
            (*A)[elemIdxSLE + i + 1][elemIdxSLE] +=
              mu * linearIntegrals[i] - 0.5 * KIntegrals[i];

            for (int j = 0; j <= i; j++)
            {
              //we use the relations
              // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
              // = mu * int_intersct x[i] * x[j] ds
              //and
              // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
              // = 0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
              (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
                += mu * quadraticIntregrals[i][j]
                  - 0.5 * ( KIntegralsX[i][j] + KIntegralsX[j][i] );
            }
          }

          if (neighborIdx > elemIdx)
          { //make sure that each facet is considered only once
            continue;
          }

          //exact evaluation of
          // int_intersct mu * jump(phi_elem,0) * jump(phi_neighbor,0) ds
          (*A)[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;

          //stiffness matrix A is symmetric
          (*A)[neighborIdxSLE][elemIdxSLE] +=
            (*A)[elemIdxSLE][neighborIdxSLE];

          for (int i = 0; i < dim; i++)
          {
            //we use the relations
            // int_intersct mu * jump(phi_elem,i) * jump(phi_neighbor,0) ds
            // = -mu * int_intersct x[i] ds
            //and
            // int_intersct avg(K*grad(phi_neighbor,i)) * jump(phi_elem,0) ds
            // = -0.5 * int_intersct K[:,i] * normal ds
            (*A)[elemIdxSLE + i + 1][neighborIdxSLE] +=
              -mu * linearIntegrals[i] + 0.5 * KIntegrals[i];

            //we use the relations
            // int_intersct mu * jump(phi_elem,0) * jump(phi_neighbor,i) ds
            // = -mu * int_intersct x[i] ds
            //and
            // int_intersct avg(K*grad(phi_neighbor,i)) * jump(phi_elem,0) ds
            // = 0.5 * int_intersct K[:,i] * normal ds
            (*A)[elemIdxSLE][neighborIdxSLE + i + 1] +=
              -mu * linearIntegrals[i] - 0.5 * KIntegrals[i];

            //stiffness matrix A is symmetric
            (*A)[neighborIdxSLE][elemIdxSLE + i + 1] +=
              (*A)[elemIdxSLE + i + 1][neighborIdxSLE];
            (*A)[neighborIdxSLE + i + 1][elemIdxSLE] +=
              (*A)[elemIdxSLE][neighborIdxSLE + i + 1];

            for (int j = 0; j <= i; j++)
            {
              //we use the relations
              // int_intersct mu * jump(phi_elem,i) * jump(phi_neighbor,j) ds
              // = -mu * int_intersct x[i] * x[j] ds
              //and
              // int_intersct avg(K*grad(phi_neighbor,j)) * jump(phi_elem,i) ds
              // = 0.5 * int_intersct x[i] * ( K[:,j] * normal ) ds
              //as well as
              // int_intersct avg(K*grad(phi_elem,i)) * jump(phi_neighbor,j) ds
              // = -0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
              (*A)[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
                -mu * quadraticIntregrals[i][j]
                - 0.5 * ( KIntegralsX[j][i] - KIntegralsX[i][j] );

              //stiffness matrix A is symmetric
              (*A)[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
                (*A)[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];

              if (i != j)
              {
                //we use the relations
                // int_intersct mu * jump(phi_elem,j) * jump(phi_neighbor,i) ds
                // = -mu * int_intersct x[i] * x[j] ds
                //and
                // int_intersct avg(K*grad(phi_neighbor,i))*jump(phi_elem,j) ds
                // = 0.5 * int_intersct x[j] * ( K[:,i] * normal ) ds
                //as well as
                // int_intersct avg(K*grad(phi_elem,j))*jump(phi_neighbor,i) ds
                // = -0.5 * int_intersct x[i] * ( K[:,j] * normal ) ds
                (*A)[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
                  -mu * quadraticIntregrals[i][j]
                  - 0.5 * ( KIntegralsX[i][j] - KIntegralsX[j][i] );

                //stiffness matrix A is symmetric
                (*A)[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
                  (*A)[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
              }
            }
          }
        }
        else //boundary facet
        {
          //fill load vector b using a quadrature rule (boundary terms)
          for (const auto& ip : rule_boundary)
          {
            //quadrature point in local and global coordinates
            const auto& qpLocal = ip.position();
            const auto& qpGlobal = intersctGeo.global(qpLocal);

            const Scalar quadratureFactor =
              ip.weight() * intersctGeo.integrationElement(qpLocal);

            const Scalar boundaryEvaluation(problem_.boundary(qpGlobal));
            const Dune::FieldMatrix<Scalar, dim, dim> KEvaluation =
              problem_.K(qpGlobal);

            //quadrature for int_intersct mu * g * phi_elem,0 ds
            b[elemIdxSLE] += quadratureFactor * mu * boundaryEvaluation;

            //quadrature for
            // int_intersct (mu * phi_elem,i - K[:,i] * normal) * g ds
            //with K*grad(phi_elem,i) = K[:,i]
            for (int i = 0; i < dim; i++)
            {
              b[elemIdxSLE + i + 1] +=
                quadratureFactor * (mu * qpGlobal[i] -
                  (KEvaluation[i] * normal)) * boundaryEvaluation;
            }
          }

          for (int i = 0; i < dim; i++)
          {
            //we evaluate the integrals
            // int_intersct mu * phi_elem,0 * phi_elem,i ds
            //and
            // int_intersct phi_elem,0 * K*grad(phi_elem,i) * normal ds
            //with K*grad(phi_elem,i) = K[:,i]
            (*A)[elemIdxSLE + i + 1][elemIdxSLE] +=
              mu * linearIntegrals[i] - KIntegrals[i];

            for (int j = 0; j <= i; j++)
            {
              //we evaluate the integrals
                // int_intersct mu * phi_elem,i * phi_elem,j ds
                //and
                // int_intersct (phi_elem,i * K[:,j]
                //    + phi_elem,j * K[:,i]) * normal ds
                //with K*grad(phi_elem,i) = K[:,i]
              (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
                mu * quadraticIntregrals[i][j]
                - ( KIntegralsX[i][j] + KIntegralsX[j][i] );
            }
          }
        }
      }

      //stiffness matrix A is symmetric
      for (int i = 0; i < dim; i++)
      {
        (*A)[elemIdxSLE][elemIdxSLE + i + 1] =
          (*A)[elemIdxSLE + i + 1][elemIdxSLE];

        for(int j = 0; j < i; j++)
        {
          (*A)[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
            (*A)[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
        }
      }
    }
  }

  //writes the solution to a vtk file
  void writeVTKoutput (const int bulkDOF) const
  {
    //storage for pressure data, for each element we store the
    //pressure at the corners of the element, the VTKFunction will
    //create the output using a linear interpolation for each element
    Dune::DynamicVector<Scalar> pressure(bulkDOF, 0.0);
    Dune::DynamicVector<Scalar> exactPressure(bulkDOF, 0.0);

    for (const auto& elem : elements(gridView_))
    {
      const int elemIdxSLE = (dim + 1)*mapper_.index(elem);
      const auto& geo = elem.geometry();

      for (int k = 0; k < geo.corners(); k++)
      {
        if (problem_.hasExactSolution())
        { //evaluation at the corners is slightly shifted towards the center of
          //the element to correctly picture discontinuities of the solution
          exactPressure[elemIdxSLE + k] = problem_.exactSolution(
            0.9999 * geo.corner(k) + 0.0001 * geo.center());
        }

        //contribution of the basis function
        // phi_elem,0 (x) = indicator(elem);
        //at the kth corner of elem
        pressure[elemIdxSLE + k] = d[elemIdxSLE];

        for (int i = 0; i < dim; i++)
        {
          //contribution of the basis function
          // phi_elem,i (x) = x[i]*indicator(elem);
          //at the kth corner of elem
          pressure[elemIdxSLE + k] += d[elemIdxSLE + i + 1] * geo.corner(k)[i];
        }
      }
    }

    //create output directory if necessary
    mkdir("data", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

    Dune::VTKWriter<GridView> vtkWriter(gridView_, Dune::VTK::nonconforming);
    vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
      new P1Function(gridView_, pressure, "pressure")));

    if (problem_.hasExactSolution())
    {
      vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
        new P1Function(gridView_, exactPressure, "exactPressure")));
    }

    vtkWriter.pwrite("pressureData", "", "data");
  }

  std::shared_ptr<Matrix> A; //stiffness matrix
  Vector b; //load vector
  Vector d; //solution vector

private:
  //writes the solution to a vtk file
  void writeVTKoutput () const
  {
    writeVTKoutput(dof_);
  }

};

#endif