#ifndef PARABOLICPDESOLVER_HPP
#define PARABOLICPDESOLVER_HPP

#include "DGReactionAssemble.hpp"
#include "PGReactionAssemble.hpp"

#include "PDESolver.hpp"


class ParabolicPDESolver : public PDESolver {
private:
  std::string quantity = "L2";

  std::string costMeasure = "size";

  std::unique_ptr<NonLinearTimeIntegrator> timeInt;

  IStochasticReactionAssemble *assemble;

protected:
  void run(SampleSolution &solution) override;

  void computeQ(SampleSolution &solution) override;

  void computeCost(SampleSolution &solution) override;

  void plotSolution(SampleSolution &solution) override;

public:
  ParabolicPDESolver(IStochasticReactionAssemble *assemble,
                     const string &quantity, const string &costMeasure) :
    PDESolver(), assemble(assemble), quantity(quantity), costMeasure(costMeasure) {

    if (typeid(*assemble).name() == typeid(PGReactionAssemble).name())
      timeInt = TimeIntegratorCreator(IMPLICIT_EULER).
        WithNonLinearSolver(new Newton(std::make_unique<GMRES>(GetPC("SuperLU")))).
        CreateUniqueNonLinearTimeIntegrator();

    if (typeid(*assemble).name() == typeid(DGReactionAssemble).name())
      timeInt = TimeIntegratorCreator(IMPLICIT_EULER).
        WithNonLinearSolver(new Newton(std::make_unique<GMRES>(GetPC(
        "PointBlockJacobi_dG")))).
        CreateUniqueNonLinearTimeIntegrator();

    if (verbose)
      mout.PrintInfo("TransportPDESolver", verbose,
                     PrintInfoEntry("Quantity", quantity),
                     PrintInfoEntry("Cost measure", costMeasure));
  }

  ~ParabolicPDESolver() override {
    delete assemble;
  }

  IAssemble *GetAssemble() const override {
    return nullptr;
//        return assemble;
  }

  IDiscretization *GetDisc() const override { return nullptr; }

  IStochasticProblem *GetProblem() const override { return nullptr; }

  void DrawSample(const SampleID &id) override {}

  std::string Name() const override { return "ParabolicPDESolver"; }
};

#endif //PARABOLICPDESOLVER_HPP