/**
 * @file   comm_serial.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 12:28:12 2011
 *
 * @brief  VMG::CommSerial
 *
 */

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

#ifdef HAVE_BOOST_FILESYSTEM
#include <boost/filesystem.hpp>
namespace fs = boost::filesystem;
#endif

#include <cstdarg>
#include <cstdio>
#include <sstream>

#include "base/discretization.hpp"
#include "base/helper.hpp"
#include "base/index.hpp"
#include "base/linked_cell_list.hpp"
#include "base/stencil.hpp"
#include "base/vector.hpp"
#include "comm/comm_serial.hpp"
#include "comm/comm.hpp"
#include "grid/multigrid.hpp"
#include "grid/tempgrid.hpp"
#include "mg.hpp"

using namespace VMG;

static char print_buffer[512];

void CommSerial::InitCommSerial()
{
  error_norm_count = 0;
  residual_count = 0;
}

void CommSerial::CommSubgrid(Grid& grid_old, Grid& grid_new)
{
  Grid::iterator iter;

  if (&grid_old == &grid_new)
    return;

  for (iter = grid_old.Iterators().CompleteGrid().Begin(); iter != grid_old.Iterators().CompleteGrid().End(); ++iter)
    grid_new(*iter) = grid_old.GetVal(*iter);
}

void CommSerial::CommAddSubgrid(Grid& grid_old, Grid& grid_new)
{
  Grid::iterator iter;

  if (&grid_old == &grid_new)
    return;

  for (iter = grid_old.Iterators().CompleteGrid().Begin(); iter != grid_old.Iterators().CompleteGrid().End(); ++iter)
    grid_new(*iter) += grid_old.GetVal(*iter);
}

Grid& CommSerial::GetFinerGrid(Multigrid& multigrid)
{
  return multigrid(multigrid.Level()+1);
}

Grid& CommSerial::GetCoarserGrid(Multigrid& multigrid)
{
  return multigrid(multigrid.Level()-1);
}

void CommSerial::CommFromGhosts(Grid& mesh)
{
  Grid::iterator iter;

  const GridIteratorSuite& iterators = mesh.Iterators();
  const Index& size = mesh.Local().Size();

  if (BoundaryConditions().X() == Periodic) {

    for (iter = iterators.Halo1().X().Begin(); iter != iterators.Halo1().X().End(); ++iter)
      mesh((*iter).X()+size.X(), (*iter).Y(), (*iter).Z()) += mesh.GetVal(*iter);

    for (iter = iterators.Halo2().X().Begin(); iter != iterators.Halo2().X().End(); ++iter)
      mesh((*iter).X()-size.X(), (*iter).Y(), (*iter).Z()) += mesh.GetVal(*iter);
  }

  if (BoundaryConditions().Y() == Periodic) {

    for (iter = iterators.Halo1().Y().Begin(); iter != iterators.Halo1().Y().End(); ++iter)
      mesh((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z()) += mesh.GetVal(*iter);

    for (iter = iterators.Halo2().Y().Begin(); iter != iterators.Halo2().Y().End(); ++iter)
      mesh((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z()) += mesh.GetVal(*iter);

  }

  if (BoundaryConditions().Z() == Periodic) {

    for (iter = iterators.Halo1().Z().Begin(); iter != iterators.Halo1().Z().End(); ++iter)
      mesh((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z()) += mesh.GetVal(*iter);

    for (iter = iterators.Halo2().Z().Begin(); iter != iterators.Halo2().Z().End(); ++iter)
      mesh((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z()) += mesh.GetVal(*iter);

  }

#ifdef DEBUG_MATRIX_CHECKS
  mesh.IsConsistent();
#endif
}

void CommSerial::CommToGhosts(Grid& mesh)
{
  Grid::iterator iter;

  const GridIteratorSuite& iterators = mesh.Iterators();
  const Index& size = mesh.Local().Size();

  if (BoundaryConditions().Z() == Periodic) {

    for (iter = iterators.Halo1().Z().Begin(); iter != iterators.Halo1().Z().End(); ++iter)
      mesh(*iter) = mesh.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()+size.Z());

    for (iter = iterators.Halo2().Z().Begin(); iter != iterators.Halo2().Z().End(); ++iter)
      mesh(*iter) = mesh.GetVal((*iter).X(), (*iter).Y(), (*iter).Z()-size.Z());

  }

  if (BoundaryConditions().Y() == Periodic) {

    for (iter = iterators.Halo1().Y().Begin(); iter != iterators.Halo1().Y().End(); ++iter)
      mesh(*iter) = mesh.GetVal((*iter).X(), (*iter).Y()+size.Y(), (*iter).Z());

    for (iter = iterators.Halo2().Y().Begin(); iter != iterators.Halo2().Y().End(); ++iter)
      mesh(*iter) = mesh.GetVal((*iter).X(), (*iter).Y()-size.Y(), (*iter).Z());

  }

  if (BoundaryConditions().X() == Periodic) {

    for (iter = iterators.Halo1().X().Begin(); iter != iterators.Halo1().X().End(); ++iter)
      mesh(*iter) = mesh.GetVal((*iter).X()+size.X(), (*iter).Y(), (*iter).Z());

    for (iter = iterators.Halo2().X().Begin(); iter != iterators.Halo2().X().End(); ++iter)
      mesh(*iter) = mesh.GetVal((*iter).X()-size.X(), (*iter).Y(), (*iter).Z());

  }

#ifdef DEBUG_MATRIX_CHECKS
  mesh.IsConsistent();
#endif
}

void CommSerial::CommToGhostsRB(Grid& grid, const int& offset)
{
  CommToGhosts(grid);
}

Grid& CommSerial::GetGlobalCoarseGrid(Multigrid& multigrid)
{
  return multigrid(multigrid.MinLevel());
}

void CommSerial::CommParticles(const Grid& grid, std::list<Particle::Particle>& particles)
{
  Factory& factory = MG::GetFactory();
  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
  vmg_float* x = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POS_ARRAY");
  vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");

  particles.clear();

  for (vmg_int i=0; i<num_particles_local; ++i)
    particles.push_back(Particle::Particle(&x[3*i], q[i], 0.0, 0.0, 0, i));
}

void CommSerial::CommLCListGhosts(const Grid& grid, Particle::LinkedCellList& lc)
{
  std::list<Particle::Particle>::iterator p_iter;
  Grid::iterator g_iter;

  for (int i=2; i>=0; --i) {

    if (BoundaryConditions()[i] == Periodic) {

      for (g_iter=lc.Iterators().Boundary1()[i].Begin(); g_iter!=lc.Iterators().Boundary1()[i].End(); ++g_iter)
	lc(*g_iter).clear();

      for (g_iter=lc.Iterators().NearBoundary1()[i].Begin(); g_iter!=lc.Iterators().NearBoundary1()[i].End(); ++g_iter)
	for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter)
	  lc.AddParticleToComplete(p_iter->Pos() + grid.Extent().Size(), p_iter->Charge(),
				   p_iter->Pot(), p_iter->Force(),
				   p_iter->Rank(), p_iter->Index());

      for (g_iter=lc.Iterators().Boundary2()[i].Begin(); g_iter!=lc.Iterators().Boundary2()[i].End(); ++g_iter)
	lc(*g_iter).clear();

      for (g_iter=lc.Iterators().NearBoundary2()[i].Begin(); g_iter!=lc.Iterators().NearBoundary2()[i].End(); ++g_iter)
	for (p_iter=lc(*g_iter).begin(); p_iter!=lc(*g_iter).end(); ++p_iter)
	  lc.AddParticleToComplete(p_iter->Pos() - grid.Extent().Size(), p_iter->Charge(),
				   p_iter->Pot(), p_iter->Force(),
				   p_iter->Rank(), p_iter->Index());

    }
  }
}

void CommSerial::CommAddPotential(std::list<Particle::Particle>& particles)
{
  std::list<Particle::Particle>::iterator iter;
  vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");

  for (iter=particles.begin(); iter!=particles.end(); ++iter)
    p[iter->Index()] += iter->Pot();
}

void CommSerial::CommAddPotential(Particle::LinkedCellList& lc)
{
  Grid::iterator lc_iter;
  std::list<Particle::Particle>::iterator p_iter;
  vmg_float* p = MG::GetFactory().GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");

  for (lc_iter=lc.Iterators().Local().Begin(); lc_iter!=lc.Iterators().Local().End(); ++lc_iter)
    for (p_iter=lc(*lc_iter).begin(); p_iter!=lc(*lc_iter).end(); ++p_iter)
      p[p_iter->Index()] += p_iter->Pot();
}

void CommSerial::PrintString(const char* format, ...)
{
  va_list args;
  va_start(args, format);
  vsprintf(print_buffer, format, args);
  printf("VMG: %s\n", print_buffer);
  va_end(args);
}

void CommSerial::PrintStringOnce(const char* format, ...)
{
  va_list args;
  va_start(args, format);
  vsprintf(print_buffer, format, args);
  printf("VMG: %s\n", print_buffer);
  va_end(args);
}

void CommSerial::PrintXML(const std::string& filename, const std::string& xml_data)
{
  pugi::xml_document doc;
  doc.load(xml_data.c_str());

  pugi::xml_node node_global = doc.prepend_child("Global");
  pugi::xml_node node_num_processes = node_global.append_child("NumProcesses");
  pugi::xml_node node_num_processes_data = node_num_processes.append_child(pugi::node_pcdata);
  node_num_processes_data.set_value(Helper::ToString(GlobalProcs().Product()).c_str());

  doc.save_file(filename.c_str());
}

void CommSerial::PrintXMLAll(const std::string& filename, const std::string& xml_data)
{
  PrintXML(filename, xml_data);
}

void CommSerial::OpenFileAndPrintHeader(std::ofstream& out, const Grid& mesh, const char* information, bool inner)
{
  char path_str[129];
  int count = OutputCount();

  sprintf(path_str, "%s%04d.dat", OutputPath().c_str(), count);

  out.open(path_str, std::ios::trunc);

  if (!out.good()) {
#ifdef DEBUG_OUTPUT
    printf("Multigrid: File %s not accessible.\n", path_str);
#endif /* DEBUG_OUTPUT */
    return;
  }

  out << "# vtk DataFile Version 2.0" << std::endl
      << count << ": " << information << std::endl
      << "ASCII" << std::endl
      << "DATASET STRUCTURED_POINTS" << std::endl;

  if (inner) {
    out << "DIMENSIONS "
	<< mesh.Local().Size().X() << " "
	<< mesh.Local().Size().Y() << " "
	<< mesh.Local().Size().Z() << std::endl;
  }else {
    out << "DIMENSIONS "
	<< mesh.Local().SizeTotal().X() << " "
	<< mesh.Local().SizeTotal().Y() << " "
	<< mesh.Local().SizeTotal().Z() << std::endl;
  }

  out	<< "ORIGIN 0 0 0" << std::endl
	<< "SPACING " << mesh.MeshWidth() << " "
	<< mesh.MeshWidth() << " "
	<< mesh.MeshWidth() << std::endl;

  if (inner)
    out << "POINT_DATA " << mesh.Local().Size().Product() << std::endl;
  else
    out << "POINT_DATA " << mesh.Local().SizeTotal().Product() << std::endl;

  out << "SCALARS residual double 1" << std::endl
      << "LOOKUP_TABLE default" << std::endl;
}

void CommSerial::PrintDefect(Grid& sol, Grid& rhs, const char* information)
{
  std::ofstream out_file;
  std::stringstream out;

  TempGrid *temp = MG::GetTempGrid();
  temp->SetProperties(sol);
  temp->ImportFromResidual(sol, rhs);

#ifdef DEBUG_MATRIX_CHECKS
  temp->IsConsistent();
#endif

  OpenFileAndPrintHeader(out_file, *temp, information, true);

  assert(out_file.good());

  if (!out_file.good())
    return;

  for (int k=temp->Local().Begin().Z(); k<temp->Local().End().Z(); k++)
    for (int j=temp->Local().Begin().Y(); j<temp->Local().End().Y(); j++)
      for (int i=temp->Local().Begin().X(); i<temp->Local().End().X(); i++)
	out << std::scientific << temp->GetVal(i,j,k) << std::endl;

  out_file << out.str();

  out_file.close();
}

void CommSerial::PrintGrid(Grid& mesh, const char* information)
{
  std::stringstream out;
  std::ofstream out_file;

  CommToGhosts(mesh);

  OpenFileAndPrintHeader(out_file, mesh, information, false);

  for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
    for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
      for (int i=0; i<mesh.Local().SizeTotal().X(); i++)
	out << std::scientific << mesh.GetVal(i,j,k) << std::endl;

  out_file << out.str();
  out_file.close();
}

void CommSerial::PrintCorrectedGrid(Grid& mesh, const char* information)
{
  std::ofstream out_file;
  std::stringstream out;

  CommToGhosts(mesh);

  OpenFileAndPrintHeader(out_file, mesh, information, false);

  for (int k=0; k<mesh.Local().SizeTotal().Z(); k++)
    for (int j=0; j<mesh.Local().SizeTotal().Y(); j++)
      for (int i=0; i<mesh.Local().SizeTotal().X(); i++)
	out << std::scientific << mesh.GetCorrectedVal(i,j,k) << std::endl;

  out_file << out.str();

  out_file.close();
}

void CommSerial::PrintInnerGrid(Grid& mesh, const char* information)
{
  std::ofstream out_file;
  std::stringstream out;

  CommToGhosts(mesh);

  OpenFileAndPrintHeader(out_file, mesh, information, true);

  for (int k=mesh.Local().Begin().Z(); k<mesh.Local().End().Z(); k++)
    for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
      for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
	out << std::scientific << mesh.GetVal(i,j,k) << std::endl;

  out_file << out.str();

  out_file.close();
}

void CommSerial::PrintInnerCorrectedGrid(Grid& mesh, const char* information)
{
  std::ofstream out_file;
  std::stringstream out;

  CommToGhosts(mesh);

  OpenFileAndPrintHeader(out_file, mesh, information, true);

  for (int k=mesh.Local().Begin().Z(); k<mesh.Local().End().Z(); k++)
    for (int j=mesh.Local().Begin().Y(); j<mesh.Local().End().Y(); j++)
      for (int i=mesh.Local().Begin().X(); i<mesh.Local().End().X(); i++)
	out << std::scientific << mesh.GetCorrectedVal(i,j,k) << std::endl;

  out_file << out.str();

  out_file.close();
}

void CommSerial::PrintAllSettings()
{
  Multigrid* mg = MG::GetFactory().Get("SOL")->Cast<Multigrid>();
  std::stringstream buf;
  std::ofstream out;

  buf << OutputPath() << "settings.txt";

  out.open(buf.str().c_str(), std::ios::trunc);

  for (int i=mg->MinLevel(); i<=mg->MaxLevel(); ++i) {

    out << "GRID LEVEL " << i << std::endl

	<< "GLOBAL" << std::endl
	<< "BEGINFINEST: " << (*mg)(i).Global().BeginFinest() << std::endl
	<< "ENDFINEST: " << (*mg)(i).Global().EndFinest() << std::endl
	<< "SIZEFINEST: " << (*mg)(i).Global().SizeFinest() << std::endl
	<< "BEGINLOCAL: " << (*mg)(i).Global().BeginLocal() << std::endl
	<< "ENDLOCAL: " << (*mg)(i).Global().EndLocal() << std::endl
	<< "SIZELOCAL: " << (*mg)(i).Global().SizeLocal() << std::endl
	<< "SIZEGLOBAL: " << (*mg)(i).Global().SizeGlobal() << std::endl
	<< "BOUNDARYTYPE: " << (*mg)(i).Global().BoundaryType() << std::endl

	<< "LOCAL" << std::endl
	<< "SIZE: " << (*mg)(i).Local().Size() << std::endl
	<< "SIZETOTAL: " << (*mg)(i).Local().SizeTotal() << std::endl
	<< "BEGIN: " << (*mg)(i).Local().Begin() << std::endl
	<< "END: " << (*mg)(i).Local().End() << std::endl
	<< "HALOBEGIN1: " << (*mg)(i).Local().HaloBegin1() << std::endl
	<< "HALOEND1: " << (*mg)(i).Local().HaloEnd1() << std::endl
	<< "HALOBEGIN2: " << (*mg)(i).Local().HaloBegin2() << std::endl
	<< "HALOEND2: " << (*mg)(i).Local().HaloEnd2() << std::endl
	<< "BOUNDARYBEGIN1: " << (*mg)(i).Local().BoundaryBegin1() << std::endl
	<< "BOUNDARYEND1: " << (*mg)(i).Local().BoundaryEnd1() << std::endl
	<< "BOUNDARYBEGIN2: " << (*mg)(i).Local().BoundaryBegin2() << std::endl
	<< "BOUNDARYEND2: " << (*mg)(i).Local().BoundaryEnd2() << std::endl
	<< "ALIGNMENTBEGIN: " << (*mg)(i).Local().AlignmentBegin() << std::endl
	<< "ALIGNMENTEND: " << (*mg)(i).Local().AlignmentEnd() << std::endl

	<< "EXTENT" << std::endl
	<< "SIZE: " << (*mg)(i).Extent().Size() << std::endl
	<< "SIZEFACTOR: " << (*mg)(i).Extent().SizeFactor() << std::endl
	<< "BEGIN: " << (*mg)(i).Extent().Begin() << std::endl
	<< "END: " << (*mg)(i).Extent().End() << std::endl
	<< "MESHWIDTH: " << (*mg)(i).Extent().MeshWidth() << std::endl;

  }

  assert(out.good());
  out.close();
}

void CommSerial::DebugPrintError(const Grid& sol, const char* information)
{
  Index i;
  std::ofstream out_file;
  std::stringstream out;

  OpenFileAndPrintHeader(out_file, sol, information, false);

  for (i.Z()=0; i.Z()<sol.Local().SizeTotal().Z(); ++i.Z())
    for (i.Y()=0; i.Y()<sol.Local().SizeTotal().Y(); ++i.Y())
      for (i.X()=0; i.X()<sol.Local().SizeTotal().X(); ++i.X())
	out << std::scientific << sol.DebugKnownSolution(i) - sol.GetVal(i) << std::endl;

  out_file << out.str();

  out_file.close();
}

void CommSerial::DebugPrintErrorNorm(Grid& sol)
{
#ifdef HAVE_LIBGSL
  char path_str[129];
  std::ofstream out;
  const vmg_float sp = sol.MeshWidth();
  const vmg_float sp_3_inv = 1.0/(sp*sp*sp);
  const int numsamples = 1e6;
  const vmg_float numsamples_inv = 1.0 / numsamples;
  vmg_float val1, val2, diff;

  vmg_float norm_L1 = 0.0;
  vmg_float norm_L2 = 0.0;
  vmg_float norm_Linfty = 0.0;
  vmg_float norm_u_L1 = 0.0;
  vmg_float norm_u_L2 = 0.0;
  vmg_float norm_u_Linfty = 0.0;

  Vector x, x_le, x_ri;
  Index ind;

  gsl_qrng *q = gsl_qrng_alloc(gsl_qrng_niederreiter_2, 3);

  sprintf(path_str, "%serror.txt", defects_path_str.c_str());

  out.open(path_str, std::ios::app);

  assert(out.good());

  CommToGhosts(sol);

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

    gsl_qrng_get(q, x.vec());

    ind = x / sp;

    x_le = x - sp * static_cast<Vector>(ind);
    x_ri = sp - x_le;

    for (int j=0; j<3; ++j)
      if (this->BoundaryConditions()[j] == Periodic)
	ind[j] += sol.Local().Begin()[j];

    val1 = ( sol.GetVal(ind) * x_ri[0] * x_ri[1] * x_ri[2] +
	     sol.GetVal(ind[0]+1, ind[1]  , ind[2]  ) * x_le[0] * x_ri[1] * x_ri[2] +
	     sol.GetVal(ind[0]  , ind[1]+1, ind[2]  ) * x_ri[0] * x_le[1] * x_ri[2] +
	     sol.GetVal(ind[0]  , ind[1]  , ind[2]+1) * x_ri[0] * x_ri[1] * x_le[2] +
	     sol.GetVal(ind[0]+1, ind[1]+1, ind[2]  ) * x_le[0] * x_le[1] * x_ri[2] +
	     sol.GetVal(ind[0]+1, ind[1]  , ind[2]+1) * x_le[0] * x_ri[1] * x_le[2] +
	     sol.GetVal(ind[0]  , ind[1]+1, ind[2]+1) * x_ri[0] * x_le[1] * x_le[2] +
	     sol.GetVal(ind[0]+1, ind[1]+1, ind[2]+1) * x_le[0] * x_le[1] * x_le[2] ) * sp_3_inv;

    val2 = sol.DebugKnownSolution(x);

    diff = fabs(val1 - val2);

    norm_L1 += diff * numsamples_inv;
    norm_L2 += diff*diff * numsamples_inv;
    norm_Linfty = std::max(norm_Linfty, diff);

    norm_u_L1 += fabs(val2) * numsamples_inv;
    norm_u_L2 += val2*val2 * numsamples_inv;
    norm_u_Linfty = std::max(norm_u_Linfty, fabs(val2));
  }

  norm_L2 = sqrt(norm_L2);
  norm_u_L2 = sqrt(norm_u_L2);

  std::cout << std::scientific << std::endl
	    << "Error L1-norm:              " << norm_L1 << std::endl
	    << "Error L2-norm:              " << norm_L2 << std::endl
	    << "Error Linfty-norm:          " << norm_Linfty << std::endl
	    << "Relative error L1-norm:     " << norm_L1 / norm_u_L1 << std::endl
	    << "Relative error L2-norm:     " << norm_L2 / norm_u_L2 << std::endl
	    << "Relative error Linfty-norm: " << norm_Linfty / norm_u_Linfty << std::endl;

  gsl_qrng_free(q);

  out << std::scientific << error_norm_count++ << " "
      << norm_L1 << " "
      << norm_L2 << " "
      << norm_Linfty << " "
      << norm_L1/norm_u_L1 << " "
      << norm_L2/norm_u_L2 << " "
      << norm_Linfty/norm_u_Linfty << std::endl;

  out.close();

#endif
}

inline int GetIndex(const Grid& grid, int i, int j, int k)
{
  if (grid.Global().BoundaryType() == LocallyRefined)
    return k + grid.Local().Size().Z() * (j + grid.Local().Size().Y() * i);
  else
    return k + grid.Local().SizeTotal().Z() * (j + grid.Local().SizeTotal().Y() * i);
}

void CommSerial::PrintGridStructureLevel(Grid& grid, std::ofstream& out)
{
  const vmg_float sp = grid.MeshWidth();
  int numLines;

  if (grid.Global().BoundaryType() == LocallyRefined)
    numLines = grid.Local().Size().X() * grid.Local().Size().Y() +
      grid.Local().Size().Y() * grid.Local().Size().Z() +
      grid.Local().Size().X() * grid.Local().Size().Z();
  else
    numLines = grid.Local().SizeTotal().X() * grid.Local().SizeTotal().Y() +
      grid.Local().SizeTotal().Y() * grid.Local().SizeTotal().Z() +
      grid.Local().SizeTotal().X() * grid.Local().SizeTotal().Z();

  out << "    <Piece";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    out << " NumberOfPoints=\"" << grid.Local().Size().Product()  << "\""
	<< " NumberOfVerts=\"" << grid.Local().Size().Product()  << "\"";
  }else {
    out << " NumberOfPoints=\"" << grid.Local().SizeTotal().Product()  << "\""
	<< " NumberOfVerts=\"" << grid.Local().SizeTotal().Product()  << "\"";
  }

  out << " NumberOfLines=\"" << numLines << "\""
      << " NumberOfStrips=\"0\""
      << " NumberOfPolys=\"0\"" << ">" << std::endl
      << "      <Points>" << std::endl
      << "        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl
      << "          ";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    for (int i=0; i<grid.Local().Size().X(); i++)
      for (int j=0; j<grid.Local().Size().Y(); j++)
	for (int k=0; k<grid.Local().Size().Z(); k++)
	  out << grid.Extent().Begin().X() + i * sp << " "
	      << grid.Extent().Begin().Y() + j * sp << " "
	      << grid.Extent().Begin().Z() + k * sp << " ";
  }else {
    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
      for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
	for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
	  out << grid.Extent().Begin().X() + i * sp << " "
	      << grid.Extent().Begin().Y() + j * sp << " "
	      << grid.Extent().Begin().Z() + k * sp << " ";
  }

  out << std::endl
      << "        </DataArray>" << std::endl
      << "      </Points>" << std::endl
      << "      <Verts>" << std::endl
      << "        <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl
      << "          ";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    for (int i=0; i<grid.Local().Size().Product(); i++)
      out << i << " ";
  }else {
    for (int i=0; i<grid.Local().SizeTotal().Product(); i++)
      out << i << " ";
  }

  out << std::endl
      << "        </DataArray>" << std::endl
      << "        <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl
      << "          ";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    for (int i=1; i<=grid.Local().Size().Product(); i++)
      out  << i << " ";
  }else {
    for (int i=1; i<=grid.Local().SizeTotal().Product(); i++)
      out  << i << " ";
  }

  out << std::endl
      << "        </DataArray>" << std::endl
      << "      </Verts>" << std::endl
      << "      <Lines>" << std::endl
      << "        <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl
      << "          ";

  if (grid.Global().BoundaryType() == LocallyRefined) {
    for (int i=0; i<grid.Local().Size().X(); i++)
      for (int j=0; j<grid.Local().Size().Y(); j++)
	out << GetIndex(grid,i,j,0) << " "
	    << GetIndex(grid,i,j,grid.Local().Size().Z()-1) << " ";

    for (int j=0; j<grid.Local().Size().Y(); j++)
      for (int k=0; k<grid.Local().Size().Z(); k++)
	out << GetIndex(grid,0,j,k) << " "
	    << GetIndex(grid,grid.Local().Size().X()-1,j,k) << " ";

    for (int i=0; i<grid.Local().Size().X(); i++)
      for (int k=0; k<grid.Local().Size().Z(); k++)
	out << GetIndex(grid,i,0,k) << " "
	    << GetIndex(grid,i,grid.Local().Size().Y()-1,k) << " ";
  }else {
    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
      for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
	out << GetIndex(grid,i,j,0) << " "
	    << GetIndex(grid,i,j,grid.Local().SizeTotal().Z()-1) << " ";

    for (int j=0; j<grid.Local().SizeTotal().Y(); j++)
      for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
	out << GetIndex(grid,0,j,k) << " "
	    << GetIndex(grid,grid.Local().SizeTotal().X()-1,j,k) << " ";

    for (int i=0; i<grid.Local().SizeTotal().X(); i++)
      for (int k=0; k<grid.Local().SizeTotal().Z(); k++)
	out << GetIndex(grid,i,0,k) << " "
	    << GetIndex(grid,i,grid.Local().SizeTotal().Y()-1,k) << " ";
  }

  out << std::endl
      << "        </DataArray>" << std::endl
      << "        <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl
      << "          ";

  for (int i=1; i<=numLines; i++)
    out << 2*i << " ";

  out << std::endl
      << "        </DataArray>" << std::endl
      << "      </Lines>" << std::endl
      << "    </Piece>" << std::endl;
}

void CommSerial::DebugPrintGridStructure(Multigrid& grid)
{
  std::ofstream out;
  char path_str[129];

  sprintf(path_str, "%sgrid.vtp", OutputPath().c_str());

  out.open(path_str, std::ios::trunc);

  if (!out.good()) {
#ifdef DEBUG_OUTPUT
    printf("Multigrid: File %s not accessible.\n", path_str);
#endif /* DEBUG_OUTPUT */
    return;
  }

  out << "<?xml version=\"1.0\"?>" << std::endl
      << "<VTKFile type=\"PolyData\">" << std::endl
      << "  <PolyData>" << std::endl;

  for (int i=grid.MinLevel(); i<=grid.MaxLevel(); i++)
    PrintGridStructureLevel(grid(i), out);

  out << "  </PolyData>" << std::endl
      << "</VTKFile>" << std::endl;

  out.close();
}

std::string CommSerial::CreateOutputDirectory()
{
#ifdef HAVE_BOOST_FILESYSTEM
  std::string path, unique_path;
  std::stringstream unique_suffix;
  int suffix_counter = 0;
  char buffer[129];
  time_t rawtime;
  struct tm *timeinfo;

  time(&rawtime);
  timeinfo = localtime(&rawtime);
  strftime(buffer, 128, "./output/%Y_%m_%d_%H_%M_%S/", timeinfo);
  path = buffer;

  if (!fs::exists("output"))
    fs::create_directory("output");

  unique_path = path;

  while (fs::exists(unique_path.c_str())) {

    unique_suffix.str("");
    unique_suffix << "_" << suffix_counter++ << "/";

    unique_path = path;
    unique_path.replace(unique_path.size()-1, 1, unique_suffix.str());

  }

  fs::create_directory(unique_path.c_str());

  return unique_path;

#else

  return "./";

#endif
}