/*
 * smoother_test.cpp
 *
 *  Created on: 20.09.2010
 *      Author: Julian Iseringhausen
 */

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

#include "comm/comm_serial.hpp"
#include "interface/interface.hpp"
#include "samples/discretization_poisson_fd.hpp"
#include "smoother/gs.hpp"
#include "smoother/gsrb.hpp"
#include "mg.hpp"

#include "interface_sinus.hpp"
#include "smoother_test.hpp"

using namespace VMG;
using VMGTests::SmootherTestSuite;

CPPUNIT_TEST_SUITE_REGISTRATION(SmootherTestSuite);

void SmootherTestSuite::setUp()
{
  Boundary boundary(Dirichlet, Dirichlet, Dirichlet);

  Comm* comm = new CommSerial(boundary);
  comm->Register("COMM");

  Discretization* discretization = new DiscretizationPoissonFD();
  discretization->Register("DISCRETIZATION");

  Interface* interface = new VMGInterfaces::InterfaceSinus(boundary, 4, 4, 0.0, 1.0);
  MG::SetInterface(interface, comm);

  interface->ImportRightHandSide(*MG::GetRhs());

  gs = new GaussSeidel();
  gsrb = new GaussSeidelRB();
}

void SmootherTestSuite::tearDown()
{
  MG::Destroy();

  delete gs;
  delete gsrb;
}

void SmootherTestSuite::GaussSeidelLEXTest()
{
  double norm;
  Multigrid& sol = *MG::GetSol();
  Multigrid& rhs = *MG::GetRhs();
  Comm& comm = *MG::GetComm();

  sol.ClearAll();

  for (int i=0; i<20; ++i) {

    gs->Run(sol, rhs, 50);
    norm = comm.ComputeResidualNorm(sol, rhs);

    if (norm < 1e-10)
      break;
  }

  CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE("Gauss-Seidel LEX smoother does not converge.", 0.0, norm, 1e-10);
}


void SmootherTestSuite::GaussSeidelRBTest()
{
  double norm;

  MG::GetSol()->ClearAll();

  for (int i=0; i<20; ++i) {

    gsrb->Run(*MG::GetSol(), *MG::GetRhs(), 50);
    norm = MG::GetComm()->ComputeResidualNorm(*MG::GetSol(), *MG::GetRhs());

    if (norm < 1e-10)
      break;
  }

  CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE("Gauss Seidel RB smoother does not converge.", 0.0, norm, 1e-10);
}