// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#pragma once

#include <vector>

#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/quadraturerules.hh>
#include <boost/numeric/mtl/mtl.hpp>

using namespace Dune;

template <class LocalView, class FineBasis>
void getLocalSystem(const LocalView& localView, const FineBasis& fine_basis,
                    mtl::dense_vector<double> const& phase,
                    mtl::dense2D<double>& elementMatrix,
                    mtl::dense_vector<double>& elementVector,
                    double eps)
{
  using namespace Dune::Indices;

  using MultiElement = MultiEntity<typename LocalView::GridView::Grid>;
  const int dim = MultiElement::dimension;
  elementMatrix.change_dim(localView.size(), localView.size());
  elementVector.change_dim(localView.size());
  set_to_zero(elementMatrix); // fills the entire matrix with zeroes
  elementVector = 0;

  // geometry of the basis element
  auto geometry = localView.element().geometry();
  const auto& localFE = localView.tree().finiteElement();

  auto phaseLocalView = fine_basis.localView();
  auto phaseIndexSet = fine_basis.localIndexSet();
  for (auto const& e : childs(localView.element(), fine_basis.gridView()))
  {
    phaseLocalView.bind(e);
    phaseIndexSet.bind(phaseLocalView);

    const auto& phaseLocalFE = phaseLocalView.tree().finiteElement();

    auto quadGeometry = e.geometry();
    const auto& quad = QuadratureRules<double, dim>::rule(quadGeometry.type(), 2*localFE.localBasis().order());
    for (const auto& quadPoint : quad)
    {
      auto local = MultiElement::local(e, localView.element(), quadPoint.position());

      const auto jacobian = geometry.jacobianInverseTransposed(local);
      const double dx = quadGeometry.integrationElement(quadPoint.position());

      static std::vector<FieldMatrix<double,1,dim> > referenceGradients;
      localFE.localBasis().evaluateJacobian(local, referenceGradients);

      static std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
      for (std::size_t i = 0; i < gradients.size(); ++i)
        jacobian.mv(referenceGradients[i][0], gradients[i]);

      static std::vector<FieldVector<double,1> > shapeValues;
      localFE.localBasis().evaluateFunction(local, shapeValues);

      static std::vector<FieldVector<double,1> > phaseShapeValues;
      phaseLocalFE.localBasis().evaluateFunction(quadPoint.position(), phaseShapeValues);

      double phase_atqp = 0.0;
      for (std::size_t i = 0; i < phaseLocalFE.size(); ++i) {
        auto row = phaseIndexSet.index(phaseLocalView.tree().localIndex(i));
        phase_atqp += phase[row[0]] * phaseShapeValues[i];
      }

      for (size_t i=0; i<localFE.size(); i++) {
        size_t row = localView.tree().localIndex(i);
        for (size_t j=0; j<localFE.size(); j++ ) {
          size_t col = localView.tree().localIndex(j);
          elementMatrix[row][col] += (
            phase_atqp * (gradients[i]*gradients[j])
            + (1.0/(eps*eps*eps)) * (1.0-phase_atqp) * (shapeValues[i] * shapeValues[j]) ) * quadPoint.weight() * dx;
        }

        elementVector[row] += phase_atqp * shapeValues[i] * quadPoint.weight() * dx;
      }
    }
  }
}


template <class Basis, class FineBasis, class Matrix, class Vector>
void assembleSystem(const Basis& basis, const FineBasis& fine_basis, Vector const& phase, Matrix& matrix, Vector& rhs, double eps)
{
  using namespace Dune::Indices;

  matrix.change_dim(basis.dimension(), basis.dimension());
  matrix = 0;

  rhs.change_dim(basis.dimension());
  rhs = 0;

  auto localView = basis.localView();
  auto localIndexSet = basis.localIndexSet();

  mtl::mat::inserter<Matrix, mtl::update_plus<double> > ins(matrix, 20);
  for (const auto& element : elements(basis.gridView())) {
    localView.bind(element);
    localIndexSet.bind(localView);

    mtl::dense2D<double> elementMatrix;
    mtl::dense_vector<double> elementVector;
    getLocalSystem(localView, fine_basis, phase, elementMatrix, elementVector, eps);

    for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
      auto row = localIndexSet.index(i);
      if (elementVector[i] != 0)
        rhs[row[0]] += elementVector[i];

      for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
        auto col = localIndexSet.index(j);
        if (elementMatrix[i][j] != 0)
          ins[row[0]][col[0]] += elementMatrix[i][j];
      }
    }
  }
}