/** \file PhaseFieldCrystal.h */

#ifndef PHASE_FIELD_CRYSTAL_H
#define PHASE_FIELD_CRYSTAL_H

#include "AMDiS.h"
#include "ExtendedProblemStat.h"
#include "Helpers.h"
#include "POperators.h"

using namespace AMDiS;

const Flag PFC_MESH_ADOPTED = 1<<2;
const Flag PFC_ADOPTED = 1<<3;

/** Phase-field Crystal problem
 */
class PhaseFieldCrystal : public ProblemIterationInterface,
                          public ProblemInstatBase
{
public:

  PhaseFieldCrystal(const std::string &name_);
  ~PhaseFieldCrystal();

  /// Initialisation of the problem.
  virtual void initialize(Flag initFlag,
                  ProblemStat *adoptProblem = NULL,
                  Flag adoptFlag = INIT_NOTHING);
  virtual void initTimeInterface();

  virtual Flag initDataFromFile(AdaptInfo *adaptInfo);
  virtual void solveInitialProblem(AdaptInfo *adaptInfo);
  virtual void transferInitialSolution(AdaptInfo *adaptInfo);

  /// Implementation of ProblemTimeInterface::initTimestep().
  virtual void initTimestep(AdaptInfo *adaptInfo);

  /// Implementation of ProblemTimeInterface::closeTimestep().
  virtual void closeTimestep(AdaptInfo *adaptInfo);

  virtual void beginIteration(AdaptInfo *adaptInfo);
  virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo= FULL_ITERATION);
  virtual void endIteration(AdaptInfo *adaptInfo);

  /// getting methods
  SystemVector *getSolution()
  {
    return prob->getSolution();
  }

  double *getTempParameter() { return &tempParameter; }

  inline Mesh* getMesh(int comp = 0) 
  {
    return prob->getMesh(comp); 
  }

  inline FiniteElemSpace* getFeSpace(int comp = 0) 
  {
    return prob->getFeSpace(comp); 
  }

  std::string getName()
  {
    return name;
  };

  int getNumProblems() 
  {
    return 1;
  };

  ProblemStat *getProblem(int number= 0)
  {
    if (number<0 || number >= getNumProblems())
      throw(std::runtime_error("problem with given number does not exist"));
    if (number == 0)
      return prob;
  };

  ProblemStat *getProblem(std::string name)
  {
    if (name == "pfc")
      return prob;
    else
      throw(std::runtime_error("problem with given name '" + name + "' does not exist"));
  };

  /// setting methods

  void setAssembleMatrixOnlyOnce_butTimestepChange(int i, int j)
  {
    fixedMatrixTimestep.push_back(std::make_pair(i,j));
  }

  void setNumberOfTimesteps(int nTimesteps_)
  {
    nTimesteps= nTimesteps_;
  }

  void writeFiles(AdaptInfo *adaptInfo, bool force);
  void serialize(std::ostream&) {};
  void deserialize(std::istream&) {};

protected:

  virtual void fillOperators();
  virtual void fillBoundaryConditions();

protected:

  ProblemStat *prob;

  bool useMobility;

  unsigned dim;
  unsigned dow;

  int nTimesteps;
  int oldMeshChangeIdx;

  double oldTimestep;
  double tempParameter;
  double r;
  double rho0;
  double density;
  double two;
  double minus2;

private:

  std::vector<std::pair<int,int> > fixedMatrixTimestep;

};


/** \ingroup MainInstat
 * \brief
 * Abstract function to calculate the pure PFC-Energy
 */
class Energy : public BinaryAbstractFunction<double,double,double>
{
  public:
    Energy() : BinaryAbstractFunction<double,double,double>(4) { }

    double operator()(const double &rho, const double &mu) const {
      return -0.25*sqr(sqr(rho)) + 0.5*rho*mu; }
};


class MobilityPfc : public AbstractFunction<double,double>
{
  public:
    MobilityPfc(double density_ = -0.3, double factor_ = 1.0) : 
      AbstractFunction<double,double>(1),
      density(density_),
      factor(factor_),
      delta(1.e-6) { }

    double operator()(const double &rho) const 
    {
      double mobility= abs(rho - 3.0*density)*factor;
      return std::max(mobility, delta);
    }
  
  protected:
    double density;
    double factor;
    double delta;
};

#endif // PHASE_FIELD_CRYSTAL_H