#include "VecelliptProject.h"
#include "AdaptStationary.h"

/// Dirichlet boundary function
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
    return exp(-10.0 * (x * x));
  }
};

/// RHS function
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}

  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
    int dim = x.getSize();
    double r2 = (x * x);
    double ux = exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * dim) * ux;
  };
};

Vecelliptdemo::Vecelliptdemo():
  vecellipt("vecellipt"),
  adaptInfo(NULL),
  adapt(NULL),
  matrixOperator00(NULL),
  matrixOperator10(NULL),
  matrixOperator11(NULL),
  rhsOperator0(NULL)
{}

Vecelliptdemo::~Vecelliptdemo() {
  if(matrixOperator00 != NULL)
    delete matrixOperator00;
  if(matrixOperator10 != NULL)
    delete matrixOperator10;
  if(matrixOperator11 != NULL)
    delete matrixOperator11;
  if(rhsOperator0 != NULL)
    delete rhsOperator0;
  if(adapt != NULL)
    delete adapt;
  if(adaptInfo != NULL)
    delete adaptInfo;
}

void Vecelliptdemo::create(string& filename) {
    // ===== init parameters =====
  Parameters::init(true, filename);

  // ===== create and init the scalar problem ===== 
  vecellipt.initialize(INIT_ALL);


  // === create adapt info ===
  adaptInfo = new AdaptInfo("vecellipt->adapt", vecellipt.getNumComponents());


  // === create adapt ===
  adapt = new AdaptStationary("vecellipt->adapt", &vecellipt, adaptInfo);

  
  // ===== create matrix operators =====
  matrixOperator00 = new Operator(vecellipt.getFeSpace(0), vecellipt.getFeSpace(0));
  matrixOperator00->addSecondOrderTerm(new Laplace_SOT);
  vecellipt.addMatrixOperator(matrixOperator00, 0, 0);

  matrixOperator10 = new Operator(vecellipt.getFeSpace(1), vecellipt.getFeSpace(0));
  matrixOperator10->addZeroOrderTerm(new Simple_ZOT);
  vecellipt.addMatrixOperator(matrixOperator10, 1, 0);

  matrixOperator11 = new Operator(vecellipt.getFeSpace(1), vecellipt.getFeSpace(1));
  matrixOperator11->addZeroOrderTerm(new Simple_ZOT(-1.0));
  vecellipt.addMatrixOperator(matrixOperator11, 1, 1);


  // ===== create rhs operator =====
  int degree = vecellipt.getFeSpace(0)->getBasisFcts()->getDegree();
  rhsOperator0 = new Operator(vecellipt.getFeSpace(0));
  rhsOperator0->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
  vecellipt.addVectorOperator(rhsOperator0, 0);

  // ===== add boundary conditions =====
  vecellipt.addDirichletBC(1, 0, 0, new G);

}

int Vecelliptdemo::solve(SolutionInformation& info) {
  int retCode = adapt->adapt();
  info.sysVec = vecellipt.getSolution();
  return retCode;
}