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


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

#include <cmath>

#include "base/helper.hpp"
#include "base/index.hpp"
#include "base/vector.hpp"
#include "comm/comm.hpp"
#include "grid/grid.hpp"
#include "grid/multigrid.hpp"
#include "grid/tempgrid.hpp"
#include "interface/interface_particles.hpp"
#include "mg.hpp"

using namespace VMG;

void InterfaceParticles::ImportRightHandSide(Multigrid& multigrid)
{
  Factory& factory = MG::GetFactory();

  const int& deg = factory.GetObjectStorageVal<int>("SPLINE_DEGREE_FCS");
  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("NUM_PARTICLES_LOCAL_FCS");
  const int deg_2 = deg / 2;

  vmg_float* x = factory.GetObjectStorageVal<vmg_float*>("X_FCS");
  vmg_float* q = factory.GetObjectStorageVal<vmg_float*>("Q_FCS");

  Index index_global, index_local, index;
  Vector pos_rel, pos_abs, index_float, grid_val;
  vmg_float foo;
  vmg_float MMtbl[3][deg];

  TempGrid* temp_grid = GetTempGrid(0);
  Grid& grid = multigrid(multigrid.GlobalMaxLevel());
  LocalIndices local = grid.Local();

  // We don't need a boundary on this grid
  local.End() -= local.Begin();
  local.Begin() = 0;
  local.BoundaryBegin1() = local.BoundaryEnd1() = 0;
  local.BoundaryBegin2() = local.BoundaryEnd2() = 0;

  //
  // Check if we need to communicate and resize halo sizes
  //
  if ((local.HaloEnd1() - local.HaloBegin1()).Product() > 0) {
    local.HaloBegin1() = 0;
    local.HaloEnd1() = local.Begin() = deg_2;
    local.End() = local.Begin() + local.Size();
  }

  if ((local.HaloEnd2() - local.HaloBegin2()).Product() > 0) {
    local.HaloBegin2() = local.End();
    local.HaloEnd2() = local.HaloBegin2() + deg_2;
  }

  local.SizeTotal() = local.Size() +
                      local.HaloEnd1() - local.HaloBegin1() +
                      local.HaloEnd2() - local.HaloBegin2();

  temp_grid->SetProperties(grid.Global(), local, grid.Extent());

  temp_grid->Clear();

  for (vmg_int i=0; i<num_particles_local; ++i) {

    //TODO: Check for correct conversions here
    index_float = (static_cast<Vector>(&x[3*i]) - grid.Extent().Begin()) / grid.Extent().MeshWidth();
    index_global = static_cast<Index>(index_float);
    index_local = temp_grid->GlobalToLocalIndex(index_global);

    Helper::MM(deg, std::modf(index_float[0], &foo), MMtbl[0]);
    Helper::MM(deg, std::modf(index_float[1], &foo), MMtbl[1]);
    Helper::MM(deg, std::modf(index_float[2], &foo), MMtbl[2]);

    for (index.X() = 0; index.X() < deg; ++index.X()) {

      grid_val.X() = q[i] * MMtbl[0][index.X()];

      for (index.Y() = 0; index.Y() < deg; ++index.Y()) {

        grid_val.Y() = grid_val.X() * MMtbl[1][index.Y()];

        for (index.Z() = 0; index.Z() < deg; ++index.Z()) {

          (*temp_grid)(index_local + index - deg_2) -= grid_val.Y() * MMtbl[2][index.Z()];

        }
      }
    }
  }

  MG::GetComm()->CommFromGhosts(*temp_grid);

  for (index.X()=0; index.X()<grid.Local().Size().X(); ++index.X())
    for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
      for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
        grid(index + grid.Local().Begin()) = (*temp_grid)(index + temp_grid->Local().Begin());
}

void InterfaceParticles::ExportSolution(Grid& grid)
{
  Factory& factory = MG::GetFactory();

  const int& deg = factory.GetObjectStorageVal<int>("SPLINE_DEGREE_FCS");
  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("NUM_PARTICLES_LOCAL_FCS");
  const Vector h_inv = 1.0 / grid.Extent().MeshWidth();

  vmg_float* x = factory.GetObjectStorageVal<vmg_float*>("X_FCS");
  vmg_float* q = factory.GetObjectStorageVal<vmg_float*>("Q_FCS");
  vmg_float* f = factory.GetObjectStorageVal<vmg_float*>("F_FCS");

  Index index, index_global, index_local;
  Vector index_float, spl, der;
  vmg_float foo;
  vmg_float MMtbl[3][deg], MdMtbl[3][deg];

  TempGrid* temp_grid = GetTempGrid(0);

  for (index.X()=0; index.X()<grid.Local().Size().X(); ++index.X())
    for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
      for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
        (*temp_grid)(index + temp_grid->Local().Begin()) = grid(index + grid.Local().Begin());

  MG::GetComm()->CommToGhosts(*temp_grid);

  for (vmg_int i=0; i<num_particles_local; ++i) {

    index_float = (static_cast<Vector>(&x[3*i]) - grid.Extent().Begin()) / grid.Extent().MeshWidth();
    index_global = static_cast<Index>(index_float);
    index_local = temp_grid->GlobalToLocalIndex(index_global);

    Helper::MdM(deg, std::modf(index_float[0], &foo), MMtbl[0], MdMtbl[0]);
    Helper::MdM(deg, std::modf(index_float[1], &foo), MMtbl[1], MdMtbl[1]);
    Helper::MdM(deg, std::modf(index_float[2], &foo), MMtbl[2], MdMtbl[2]);

    for (int j=0; j<3; ++j)
      f[3*i+j] = 0.0;

    for (index.X() = 0; index.X() < deg; ++index.X()) {

      spl.X() = MMtbl[0][index.X()];
      der.X() = MdMtbl[0][index.X()];

      for (index.Y() = 0; index.Y() < deg; ++index.Y()) {

        spl.Y() = MMtbl[1][index.Y()];
        der.Y() = MdMtbl[1][index.Y()];

        for (index.Z() = 0; index.Z() < deg; ++index.Z()) {

          spl.Z() = MMtbl[2][index.Z()];
          der.Z() = MdMtbl[2][index.Z()];

          f[3*i  ] += temp_grid->GetVal(index + temp_grid->Local().Begin()) *
                      der[0]*spl[1]*spl[2] * h_inv[0];

          f[3*i+1] += temp_grid->GetVal(index + temp_grid->Local().Begin()) *
                      spl[0]*der[1]*spl[2] * h_inv[1];

          f[3*i+2] += temp_grid->GetVal(index + temp_grid->Local().Begin()) *
                      spl[0]*spl[1]*der[2] * h_inv[2];

        }
      }
    }

    for (int j=0; j<3; ++j)
      f[3*i+j] *= q[i];

  }

}