#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <dune/common/exceptions.hh>

#include <dune/solvers/common/arithmetic.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>

#include "enums.hh"
#include "enumparser.hh"

#include "fixedpointiterator.hh"

template <class Factory, class StateUpdater, class VelocityUpdater,
          class ErrorNorm>
FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::
    FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
                       std::shared_ptr<Nonlinearity> globalFriction,
                       ErrorNorm const &errorNorm)
    : factory_(factory),
      parset_(parset),
      globalFriction_(globalFriction),
      fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
      fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
      lambda_(parset.get<double>("v.fpi.lambda")),
      velocityMaxIterations_(parset.get<size_t>("v.solver.maximumIterations")),
      velocityTolerance_(parset.get<double>("v.solver.tolerance")),
      verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
      errorNorm_(errorNorm) {}

template <class Factory, class StateUpdater, class VelocityUpdater,
          class ErrorNorm>
int FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run(
    std::shared_ptr<StateUpdater> stateUpdater,
    std::shared_ptr<VelocityUpdater> velocityUpdater,
    Matrix const &velocityMatrix, Vector const &velocityRHS,
    Vector &velocityIterate) {
  auto multigridStep = factory_.getSolver();

  EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
  LoopSolver<Vector> velocityProblemSolver(
      multigridStep, velocityMaxIterations_, velocityTolerance_, &energyNorm,
      verbosity_, false); // absolute error

  size_t fixedPointIteration;
  ScalarVector alpha;
  stateUpdater->extractAlpha(alpha);
  for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
       ++fixedPointIteration) {
    // solve a velocity problem
    globalFriction_->updateAlpha(alpha);
    ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
                                velocityRHS, velocityIterate);
    BlockProblem velocityProblem(parset_, convexProblem);
    multigridStep->setProblem(velocityIterate, velocityProblem);
    velocityProblemSolver.preprocess();
    velocityProblemSolver.solve();

    Vector v_m;
    velocityUpdater->extractOldVelocity(v_m);
    v_m *= 1.0 - lambda_;
    Arithmetic::addProduct(v_m, lambda_, velocityIterate);

    // solve a state problem
    stateUpdater->solve(v_m);
    ScalarVector newAlpha;
    stateUpdater->extractAlpha(newAlpha);

    if (errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
      fixedPointIteration++;
      break;
    }
    alpha = newAlpha;
  }
  if (fixedPointIteration == fixedPointMaxIterations_)
    DUNE_THROW(Dune::Exception, "FPI failed to converge");

  velocityUpdater->postProcess(velocityIterate);

  return fixedPointIteration;
}

#include "fixedpointiterator_tmpl.cc"