source: src/grid/tempgrid.cpp@ 64ba929

/**
 * @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/grid_index_translations.hpp"
#include "grid/tempgrid.hpp"
#include "interface/interface.hpp"
#include "mg.hpp"

using namespace VMG;

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

  iterators.SetSubgrids(rhs.Local());

  level = rhs.Level();

  Allocate();
}

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

  iterators.SetSubgrids(local_);

  Allocate();
}

void TempGrid::SetPropertiesToGlobalCoarseGrid()
{
  Interface* interface = MG::GetInterface();

  global = interface->Global().back();
  extent = interface->Extent().back();

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

    if (interface->BoundaryConditions()[i] == Dirichlet ||
        interface->BoundaryConditions()[i] == Open) {

      local.Size()[i] = global.SizeGlobal()[i] - 2;
      local.SizeTotal()[i] = global.SizeGlobal()[i];
      local.Begin()[i] = 1;
      local.End()[i] = local.SizeTotal()[i] - 1;
      local.HaloBegin1()[i] = 0;
      local.HaloEnd1()[i] = 0;
      local.HaloBegin2()[i] = 0;
      local.HaloEnd2()[i] = 0;
      local.BoundaryBegin1()[i] = 0;
      local.BoundaryEnd1()[i] = 1;
      local.BoundaryBegin2()[i] = local.SizeTotal()[i] - 1;
      local.BoundaryEnd2()[i] = local.SizeTotal()[i];

    }else if (interface->BoundaryConditions()[i] == Periodic) {

      local.Size()[i] = global.SizeGlobal()[i];
      local.SizeTotal()[i] = global.SizeGlobal()[i];
      local.Begin()[i] = 0;
      local.End()[i] = global.SizeGlobal()[i];
      local.HaloBegin1()[i] = 0;
      local.HaloEnd1()[i] = 0;
      local.HaloBegin2()[i] = 0;
      local.HaloEnd2()[i] = 0;
      local.BoundaryBegin1()[i] = 0;
      local.BoundaryEnd1()[i] = 0;
      local.BoundaryBegin2()[i] = 0;
      local.BoundaryEnd2()[i] = 0;

    }

  }

  level = interface->MinLevel();

  iterators.SetSubgrids(local);

  Allocate();
}

void TempGrid::SetPropertiesToFiner(const Grid& grid, const Boundary& boundary)
{
  assert(grid.Father() != NULL);
  assert(grid.Level() < grid.Father()->MaxLevel());

  const Index off = GridIndexTranslations::EndOffset(boundary);

  this->Global().BeginFinest() = grid.Global().BeginFinest();
  this->Global().EndFinest() = grid.Global().EndFinest();
  this->Global().SizeFinest() = grid.Global().SizeFinest();
  this->Global().SizeGlobal() = 2 * (grid.Global().SizeGlobal()-off) + off;
  this->Global().BeginLocal() = grid.Global().BeginLocal();
  this->Global().EndLocal() = grid.Global().EndLocal();

  GridIndexTranslations::CoarseToFine(this->Global().BeginLocal(), this->Global().EndLocal(), this->Global().SizeGlobal());

  this->Global().SizeLocal() = this->Global().EndLocal() - this->Global().BeginLocal();
  this->Global().BoundaryType() = (*(grid.Father()))(grid.Level()+1).Global().BoundaryType();

  this->Local().Size() = this->Global().SizeLocal();
  this->Local().SizeTotal() = this->Global().SizeLocal();
  this->Local().Begin() = 0;
  this->Local().End() = this->Global().SizeLocal();

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

    if (grid.Local().HaloEnd1()[i] - grid.Local().HaloBegin1()[i] > 0) {
      this->Local().SizeTotal()[i] += 1;
      this->Local().Begin()[i] += 1;
      this->Local().End()[i] += 1;
      this->Local().HaloBegin1()[i] = 0;
      this->Local().HaloEnd1()[i] = 1;
    }else {
      this->Local().HaloBegin1()[i] = 0;
      this->Local().HaloEnd1()[i] = 0;
    }

    if (grid.Local().BoundaryEnd1()[i] - grid.Local().BoundaryBegin1()[i] > 0) {
      this->Local().Size()[i] -= 1;
      this->Local().Begin()[i] += 1;
      this->Local().BoundaryBegin1()[i] = 0;
      this->Local().BoundaryEnd1()[i] = 1;
    }else {
      this->Local().BoundaryBegin1()[i] = 0;
      this->Local().BoundaryEnd1()[i] = 0;
    }

    if (grid.Local().HaloEnd2()[i] - grid.Local().HaloBegin2()[i] > 0) {
      this->Local().SizeTotal()[i] += 1;
      this->Local().HaloBegin2()[i] = this->Local().End()[i];
      this->Local().HaloEnd2()[i] = this->Local().End()[i] + 1;
    }else {
      this->Local().HaloBegin2()[i] = 0;
      this->Local().HaloEnd2()[i] = 0;
    }

    if (grid.Local().BoundaryEnd2()[i] - grid.Local().BoundaryBegin2()[i] > 0) {
      this->Local().Size()[i] -= 1;
      this->Local().End()[i] -= 1;
      this->Local().BoundaryBegin2()[i] = this->Local().End()[i];
      this->Local().BoundaryEnd2()[i] = this->Local().End()[i]+1;
    }else {
      this->Local().BoundaryBegin2()[i] = 0;
      this->Local().BoundaryEnd2()[i] = 0;
    }

  }

  this->Local().AlignmentBegin() = this->Local().Begin();
  this->Local().AlignmentEnd() = this->local.End();

  this->Extent().Size() = grid.Extent().Size();
  this->Extent().SizeFactor() = grid.Extent().SizeFactor();
  this->Extent().Begin() = grid.Extent().Begin();
  this->Extent().End() = grid.Extent().End();
  this->Extent().MeshWidth() = 0.5 * grid.Extent().MeshWidth();

  this->level = grid.Level() + 1;

  iterators.SetSubgrids(local);

  Allocate();
}

void TempGrid::SetPropertiesToCoarser(const Grid& grid, const Boundary& boundary)
{
  assert(grid.Father() != NULL);
  assert(grid.Level() > grid.Father()->MinLevel());

  const Index off = GridIndexTranslations::EndOffset(boundary);

  this->level = grid.Level() - 1;

  this->Global().BeginFinest() = grid.Global().BeginFinest();
  this->Global().EndFinest() = grid.Global().EndFinest();
  this->Global().SizeFinest() = grid.Global().SizeFinest();

  this->Global().BeginLocal() = grid.Global().BeginLocal();
  this->Global().EndLocal() = grid.Global().EndLocal();

  GridIndexTranslations::FineToCoarse(this->Global().BeginLocal(), this->Global().EndLocal());

  this->Global().SizeLocal() = this->Global().EndLocal() - this->Global().BeginLocal();
  this->Global().SizeGlobal() = (grid.Global().SizeGlobal()-off)/2 + off;
  this->Global().BoundaryType() = (*(grid.Father()))(grid.Level()-1).Global().BoundaryType();

  this->Local().SizeTotal() = this->Global().SizeLocal();
  this->Local().Size() = this->Global().SizeLocal();
  this->Local().Begin() = 0;
  this->Local().End() = this->Global().SizeLocal();

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

    if (grid.Local().HaloEnd1()[i] - grid.Local().HaloBegin1()[i] > 0) {
      this->Local().SizeTotal()[i] += 1;
      this->Local().Begin()[i] += 1;
      this->Local().End()[i] += 1;
      this->Local().HaloBegin1()[i] = 0;
      this->Local().HaloEnd1()[i] = 1;
    }else {
      this->Local().HaloBegin1()[i] = 0;
      this->Local().HaloEnd1()[i] = 0;
    }

    if (grid.Local().BoundaryEnd1()[i] - grid.Local().BoundaryBegin1()[i] > 0) {
      this->Local().Size()[i] -= 1;
      this->Local().Begin()[i] += 1;
      this->Local().BoundaryBegin1()[i] = 0;
      this->Local().BoundaryEnd1()[i] = 1;
    }else {
      this->Local().BoundaryBegin1()[i] = 0;
      this->Local().BoundaryEnd1()[i] = 0;
    }

    if (grid.Local().HaloEnd2()[i] - grid.Local().HaloBegin2()[i] > 0) {
      this->Local().SizeTotal()[i] += 1;
      this->Local().HaloBegin2()[i] = this->Local().End()[i];
      this->Local().HaloEnd2()[i] = this->Local().End()[i] + 1;
    }else {
      this->Local().HaloBegin2()[i] = 0;
      this->Local().HaloEnd2()[i] = 0;
    }

    if (grid.Local().BoundaryEnd2()[i] - grid.Local().BoundaryBegin2()[i] > 0) {
      this->Local().Size()[i] -= 1;
      this->Local().End()[i] -= 1;
      this->Local().BoundaryBegin2()[i] = this->Local().End()[i];
      this->Local().BoundaryEnd2()[i] = this->Local().End()[i]+1;
    }else {
      this->Local().BoundaryBegin2()[i] = 0;
      this->Local().BoundaryEnd2()[i] = 0;
    }

  }

  this->Local().AlignmentBegin() = this->Local().Begin();
  this->Local().AlignmentEnd() = this->local.End();

  this->Extent().Size() = grid.Extent().Size();
  this->Extent().SizeFactor() = grid.Extent().SizeFactor();
  this->Extent().Begin() = grid.Extent().Begin();
  this->Extent().End() = grid.Extent().End();
  this->Extent().MeshWidth() = 2.0 * grid.Extent().MeshWidth();

  iterators.SetSubgrids(local);

  Allocate();
}

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

  Grid::iterator grid_iter;
  Stencil::iterator stencil_iter;
  vmg_float residual;
  Index i;

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

  this->Clear();

  MG::GetComm()->CommToGhosts(sol);

  for (grid_iter=Iterators().Local().Begin(); grid_iter!=Iterators().Local().End(); ++grid_iter) {

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

    for (stencil_iter=A.begin(); stencil_iter!=A.end(); ++stencil_iter)
      residual -= prefactor * stencil_iter->Val() * sol.GetVal(*grid_iter + stencil_iter->Disp());

    (*this)(*grid_iter) = residual;

  }

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

void TempGrid::Allocate()
{
  const int size = local.SizeTotal().Product();

  if (size > size_max) {
    size_max = size;
    delete [] grid;
    grid = new vmg_float[size];
  }
}

TempGrid::TempGrid()
{
  size_max = 0;
}

TempGrid::~TempGrid()
{
}