#include "AMDiS.h"

using namespace std;
using namespace AMDiS;

// ===========================================================================
// ===== function definitions ================================================
// ===========================================================================

/** \brief
 * RHS function
 */
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
  MEMORY_MANAGED(F);

  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {};

  /** \brief
   * Implementation of AbstractFunction::operator().
   */
  double operator()(const WorldVector<double>& x) const {
    return -2.0 * x[0];
  }
};

/** \brief
 * Implementation of a scalar problem. In \ref buildBeforeCoarsen() parametric
 * coordinates for the vertices are filled in a DOFVector. This DOFVector is
 * used in \ref parametric to parametrize elements while mesh traversal. 
 */
class ParametricSphere : public ProblemScal 
{
public:
  /** \brief
   * constructor
   */
  ParametricSphere(const char *name) 
    : ProblemScal(name),
      parametric(NULL)
  {};

  /** \brief
   * destructor
   */
  ~ParametricSphere() {
    for (int i = 0; i < Global::getGeo(WORLD); i++) {
      DELETE parametricCoords[i];
    }
    DELETE parametric;
  }

  /** \brief
   * initialization of the base class and creation of \ref parametric.
   */
  void initialize(Flag initFlag,
		  ProblemScal *adoptProblem = NULL,
		  Flag adoptFlag = INIT_NOTHING)
  {
    ProblemScal::initialize(initFlag, adoptProblem, adoptFlag);

    // ===== create projection =====
    WorldVector<double> ballCenter;
    ballCenter.set(0.0);
    NEW BallProject(1, 
		    VOLUME_PROJECTION, 
		    ballCenter, 
		    1.0);

    // ===== create coords vector =====
    int i;
    for(i = 0; i < Global::getGeo(WORLD); i++) {
      parametricCoords[i] = NEW DOFVector<double>(this->getFESpace(), 
						  "parametric coords");
    }

    // ===== create parametric object =====
    parametric = NEW ParametricFirstOrder(&parametricCoords);

    // ===== enable parametric traverse =====
    this->getMesh()->setParametric(parametric);
  };

  /** \brief
   * Implementation of ProblemStatBase::buildBeforeCoarsen().
   */
  void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag flag) {
    FUNCNAME("ParametricSphere::buildAfterCoarsen()");
    MSG("calculation of parametric coordinates\n");
    int i, j;
    int preDOFs = this->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX);
    int dim = this->getMesh()->getDim();
    int dow = Global::getGeo(WORLD);
    WorldVector<double> vertexCoords;
    const DegreeOfFreedom **dof;
    DegreeOfFreedom vertexIndex;

    // ===== disable parametric traverse =====
    this->getMesh()->setParametric(NULL);

    TraverseStack stack;
    ElInfo *elInfo = NULL;
    elInfo = stack.traverseFirst(this->getMesh(), -1, 
				 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while(elInfo) {
      dof = elInfo->getElement()->getDOF();
      for(i = 0; i < dim + 1; i++) {
	vertexCoords = elInfo->getCoord(i);
	Projection::getProjection(1)->project(vertexCoords);
	for(j = 0; j < dow; j++) {
	  (*(parametricCoords[j]))[dof[i][preDOFs]] = vertexCoords[j];
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }
    

    // ===== enable parametric traverse =====
    this->getMesh()->setParametric(parametric);
  };

protected:
  /** \brief
   * DOFVector storing parametric coordinates.
   */
  WorldVector<DOFVector<double>*> parametricCoords;

  /** \brief
   * Parametrizes vertex coordinates while mesh traversal.
   */
  ParametricFirstOrder *parametric;
};

// ===========================================================================
// ===== main program ========================================================
// ===========================================================================

int main(int argc, char* argv[])
{
  FUNCNAME("parametric main");

  // ===== check for init file =====
  TEST_EXIT(argc == 2)("usage: parametric initfile\n");

  // ===== init parameters =====
  Parameters::init(false, argv[1]);

  // ===== create and init the scalar problem ===== 
  ParametricSphere parametric("parametric");
  parametric.initialize(INIT_ALL);

  // === create adapt info ===
  AdaptInfo *adaptInfo = NEW AdaptInfo("parametric->adapt", 1);

  // === create adapt ===
  AdaptStationary *adapt = NEW AdaptStationary("parametric->adapt",
					       &parametric,
					       adaptInfo);
  
  // ===== create matrix operator =====
  Operator matrixOperator(Operator::MATRIX_OPERATOR,
			  parametric.getFESpace());
  matrixOperator.addSecondOrderTerm(NEW Laplace_SOT);
  parametric.addMatrixOperator(&matrixOperator);

  // ===== create rhs operator =====
  Operator rhsOperator(Operator::VECTOR_OPERATOR,
		       parametric.getFESpace());

  int degree = parametric.getFESpace()->getBasisFcts()->getDegree();

  rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree)));
  parametric.addVectorOperator(&rhsOperator);

  // ===== start adaption loop =====
  adapt->adapt();

  parametric.writeFiles(adaptInfo, true);

  return 0;
}