/**
 * @file   interface_fcs.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 12:56:20 2011
 *
 * @brief  VMG::InterfaceFCS
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#ifndef HAVE_MPI
#error MPI is needed to use the Scafacos interface.
#endif

#include <mpi.h>

#include "base/object.hpp"
#include "comm/comm_serial.hpp"
#include "interface/interface_fcs.h"
#include "interface/interface_particles.hpp"
#include "level/level_operator_cs.hpp"
#include "samples/discretization_poisson_fd.hpp"
#include "samples/stencils.hpp"
#include "samples/techniques.hpp"
#include "smoother/gsrb.cpp"
#ifdef HAVE_LAPACK
#include "solver/dgesv.hpp"
#include "solver/dsysv.hpp"
#else
#include "solver/givens.hpp"
#endif
#include "solver/solver_dirichlet.hpp"
#include "solver/solver_periodic.hpp"

using namespace VMG;

void VMGInterfaceFCSInit(vmg_int level, vmg_int periodic, vmg_int spline_degree,
			 vmg_int max_iter, vmg_int smoothing_steps, vmg_int gamma,
			 vmg_float precision, vmg_float box_begin, vmg_float box_end,
			 MPI_Comm mpi_comm)
{
  const BC bc = (periodic ? Periodic : Dirichlet);

  /*
   * Register all the parameters for later use.
   */
  new ObjectStorage<int>("LEVEL_FCS", level);
  new ObjectStorage<int>("PERIODIC_FCS", periodic);
  new ObjectStorage<int>("SPLINE_DEGREE_FCS", spline_degree);
  new ObjectStorage<int>("MAX_ITER_FCS", max_iter);
  new ObjectStorage<int>("SMOOTHING_STEPS_FCS", smoothing_steps);
  new ObjectStorage<int>("GAMMA_FCS", gamma);
  new ObjectStorage<vmg_float>("PRECISION_FCS", precision);
  new ObjectStorage<vmg_float>("BOX_BEGIN_FCS", box_begin);
  new ObjectStorage<vmg_float>("BOX_END_FCS", box_end);
  new ObjectStorage<MPI_Comm>("MPI_COMM_FCS", mpi_comm);

  /*
   * Register communication class.
   * For now, parallel version is not yet released.
   */
  Comm* comm = new CommSerial(bc);
  comm->Register("COMM");

  /*
   * Register particle interface.
   */
  Interface* interface = new VMG::InterfaceParticles(bc, 2, level, box_begin, box_end);
  MG::SetInterface(interface, comm);

  /*
   * Use a finite difference discretization of the Poisson equation
   */
  Discretization* discretization = new DiscretizationPoissonFD();
  discretization->Register("DISCRETIZATION");

  /*
   * Use a correction scheme multigrid algorithm with full weight stencils.
   * Since we use the Gauss-Seidel RB smoother here, this may be replaced
   * with a half-weight version once debugging is finished.
   */
  LevelOperator* lop = new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
  lop->Register("LEVEL_OPERATOR");

  /*
   * Use Gauss-Seidel Red-Black ordering
   */
  Smoother* smoother = new GaussSeidelRB();
  smoother->Register("SMOOTHER");

  /*
   * Use LAPACK solver when available, otherwise use Givens rotations.
   */
#ifdef HAVE_LAPACK
  Solver* solver = (periodic ?
		    static_cast<Solver*>(new DSYSV<SolverPeriodic>()) :
		    static_cast<Solver*>(new DGESV<SolverDirichlet>()));
#else
  Solver* solver = (periodic ?
		    static_cast<Solver*>(new Givens<SolverPeriodic>()) :
		    static_cast<Solver*>(new Givens<SolverDirichlet>()));
#endif
  solver->Register("SOLVER");

  /*
   * Set commands for the actual multigrid cycle
   */
  if (periodic)
    Techniques::SetCorrectionSchemePeriodic(2, level, gamma);
  else
    Techniques::SetCorrectionSchemeDirichlet(2, level, gamma);

  /*
   * Register required parameters
   */
  new ObjectStorage<int>("PRESMOOTHSTEPS", smoothing_steps);
  new ObjectStorage<int>("POSTSMOOTHSTEPS", smoothing_steps);
  new ObjectStorage<vmg_float>("PRECISION", precision);
  new ObjectStorage<int>("MAX_ITERATION", max_iter);

  /*
   * Check whether the library is correctly initialized now.
   */
  MG::IsInitialized();
}

void VMGInterfaceFCSRun(vmg_float* x, vmg_float* q, vmg_float* p, vmg_float* f, vmg_int num_particles_local)
{
  /*
   * Register parameters for later use.
   */
  new ObjectStorage<vmg_float*>("X_FCS", x);
  new ObjectStorage<vmg_float*>("Q_FCS", q);
  new ObjectStorage<vmg_float*>("P_FCS", p);
  new ObjectStorage<vmg_float*>("F_FCS", f);
  new ObjectStorage<vmg_int>("NUM_PARTICLES_LOCAL_FCS", num_particles_local);

  /*
   * Start the multigrid solver
   */
  MG::Solve();
}

void VMGInterfaceFCSDestroy(void)
{
  /*
   * Delete all data.
   */
  MG::Destroy();
}