/**
 * @file   tempgrid.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 12:55:05 2011
 *
 * @brief  VMG::TempGrid
 *
 */


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

#include "base/discretization.hpp"
#include "base/stencil.hpp"
#include "comm/comm.hpp"
#include "grid/tempgrid.hpp"
#include "mg.hpp"

using namespace VMG;

void TempGrid::SetProperties(const Grid& rhs)
{
  local = rhs.Local();
  global = rhs.Global();
  extent = rhs.Extent();

  delete [] grid;

  grid = new vmg_float[local.SizeTotal().Product()];
}

void TempGrid::SetProperties(const GlobalIndices& global_, const LocalIndices& local_, const SpatialExtent& extent_)
{
  local = local_;
  global = global_;
  extent = extent_;

  delete [] grid;

  grid = new vmg_float[local.SizeTotal().Product()];
}

void TempGrid::ImportFromResidual(const Grid& sol, const Grid& rhs)
{
#ifdef DEBUG
  this->IsCompatible(sol);
  this->IsCompatible(rhs);
#endif

  vmg_float residual;
  Index i;

  const vmg_float prefactor = MG::GetDiscretization()->OperatorPrefactor(sol);
  const Stencil& A = MG::GetDiscretization()->GetStencil();

  this->Clear();

  for (i.X()=Local().Begin().X(); i.X()<Local().End().X(); ++i.X())
    for (i.Y()=Local().Begin().Y(); i.Y()<Local().End().Y(); ++i.Y())
      for (i.Z()=Local().Begin().Z(); i.Z()<Local().End().Z(); ++i.Z()) {

	residual = rhs.GetVal(i) - prefactor * A.GetDiag() * sol.GetVal(i);

	for (Stencil::iterator iter=A.begin(); iter!=A.end(); iter++)
	  residual -= prefactor * iter->val * sol.GetVal(i+iter);

	(*this)(i) = residual;

      }

  MG::GetComm()->CommToGhosts(*this);
}

TempGrid::~TempGrid()
{
}