/*
 * interface_sinus.cpp
 *
 *  Created on: 31.03.2011
 *      Author: Julian Iseringhausen
 */

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

#ifndef _USE_MATH_DEFINES
#define _USE_MATH_DEFINES
#endif

#include <cmath>
#include <iostream>

#include "grid/grid.hpp"
#include "grid/multigrid.hpp"

#include "interface_sinus.hpp"

using namespace VMG;
using VMGInterfaces::InterfaceSinus;

void InterfaceSinus::ImportRightHandSide(Multigrid& grid)
{
  Index i;
  Vector pos;

  Grid& grid_l = grid(grid.GlobalMaxLevel());

  const double& sp = grid_l.MeshWidth();
  const double factor = 2.0 * M_PI * grid_l.Extent().Size().Max();

  Index begin, end, begin_local;

  for (int j=0; j<3; ++j) {

    begin[j] = (grid_l.Local().HaloEnd1()[j] - grid_l.Local().HaloBegin1()[j] > 0 ? grid_l.Local().Begin()[j] : 0);
    end[j] = (grid_l.Local().HaloEnd2()[j] - grid_l.Local().HaloBegin2()[j] > 0 ? grid_l.Local().End()[j] : grid_l.Local().SizeTotal()[j]);
    begin_local[j] = grid_l.Global().BeginLocal()[j] - (grid_l.Local().HaloEnd1()[j] - grid_l.Local().HaloBegin1()[j] > 0 ? 1 : 0);
  }

  grid_l.ClearBoundary();

  for (i.X()=begin.X(); i.X()<end.X(); ++i.X())
    for (i.Y()=begin.Y(); i.Y()<end.Y(); ++i.Y())
      for (i.Z()=begin.Z(); i.Z()<end.Z(); ++i.Z()) {

	pos = grid_l.Extent().Begin() + sp * static_cast<Vector>(begin_local + i);

	grid_l(i) = 3.0 * factor * factor * sin(factor * pos.X()) * sin(factor * pos.Y()) * sin(factor * pos.Z());

      }
}

void InterfaceSinus::ExportSolution(Grid& grid)
{
}