source: src/grid/multigrid.cpp@ 66f24d

/**
 * @file   multigrid.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 12:54:37 2011
 *
 * @brief  VMG::Multigrid
 *
 */

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

#include "base/helper.hpp"
#include "base/index.hpp"
#include "comm/comm.hpp"
#include "grid/grid.hpp"
#include "grid/grid_indexing.hpp"
#include "grid/multigrid.hpp"
#include "interface/interface.hpp"

using namespace VMG;

Multigrid::Multigrid(Comm* comm, const Interface* interface)
{
  bool active;
  Index points, remainder;
  Index procs_l;
  LocalIndices local_l;
  GlobalIndices global_l;

  const Index& pos = comm->Pos();
  const Index& procs = comm->Procs();
  const Index points_min = 5 * procs;
  const std::vector<GlobalIndices>& global = interface->Global();
  const std::vector<SpatialExtent>& extent = interface->Extent();

  for (unsigned int i=0; i<global.size(); ++i) {

    procs_l = global[i].Size() / points_min;

    active = true;
    for (int j=0; j<3; ++j) {
      procs_l[j] = std::max(procs_l[j], 1);
      procs_l[j] = std::min(procs[j], procs_l[j]);
      active &= pos[j] < procs_l[j];
    }

    if (global[i].BoundaryType() == GlobalMax)
      levelGlobalMax = interface->MaxLevel() - i;

    if (active) {

      // Set boundary type
      global_l.BoundaryType() = global[i].BoundaryType();

      // Compute global size
      global_l.Size() = global[i].Size() / procs_l;
      remainder = global[i].Size() % procs_l;

      for (int j=0; j<3; ++j)
        if (pos[j] < remainder[j])
          ++(global_l.Size()[j]);

      // Compute global begin
      global_l.Begin() = pos * global_l.Size();

      for (int j=0; j<3; ++j)
        if (pos[j] >= remainder[j])
          global_l.Begin()[j] += remainder[j];

      // Compute global end
      global_l.End() = global_l.Begin() + global_l.Size();

      // Set alignment
      for (int j=0; j<3; ++j) {
        global_l.AlignmentBegin()[j] = std::max(0, global[i].AlignmentBegin()[j] - global_l.Begin()[j]);
        global_l.AlignmentEnd()[j] = std::min(global_l.Size()[j], global[i].AlignmentEnd()[j] - global_l.Begin()[j]);
      }

      switch (comm->BoundaryConditions())
        {
        case Dirichlet:
          local_l = InitDirichlet(global_l, extent[i], pos, procs_l);
          break;
        case Periodic:
          local_l = InitPeriodic(global_l, extent[i], pos, procs_l);
          break;
        case Quasiperiodic:
          local_l = InitDirichlet(global_l, extent[i], pos, procs_l);
          break;
        }

      grids.push_back(new Grid(global_l, local_l, extent[i]));

    }else {

      global_l.Size() = 0;
      global_l.Begin() = 0;
      global_l.End() = 0;
      global_l.AlignmentBegin() = 0;
      global_l.AlignmentEnd() = 0;
      global_l.BoundaryType() = EmptyGrid;

      local_l.Size() = 0;
      local_l.SizeTotal() = 0;
      local_l.Begin() = 0;
      local_l.End() = 0;
      local_l.HaloBegin1() = 0;
      local_l.HaloBegin2() = 0;
      local_l.HaloEnd1() = 0;
      local_l.HaloEnd2() = 0;
      local_l.BoundaryBegin1() = 0;
      local_l.BoundaryBegin2() = 0;
      local_l.BoundaryEnd1() = 0;
      local_l.BoundaryEnd2() = 0;
      local_l.AlignmentBegin() = 0;
      local_l.AlignmentEnd() = 0;

      grids.push_back(new Grid(global_l, local_l, extent[i]));

    }
  }

  numLevels = grids.size();

  levelMax = interface->MaxLevel();
  levelMin = levelMax - numLevels + 1;

  levelIndex = 0;
  levelCurrent = levelMax;
}

LocalIndices Multigrid::InitDirichlet(const GlobalIndices& global, const SpatialExtent& extent, const Index& pos, const Index& procs)
{
  LocalIndices local;

  local.HaloBegin1() = local.HaloEnd1() = 0;
  local.HaloBegin2() = local.HaloEnd2() = 0;
  local.BoundaryBegin1() = local.BoundaryEnd1() = 0;
  local.BoundaryBegin2() = local.BoundaryEnd2() = 0;

  local.SizeTotal() = global.End() - global.Begin() + (global.BoundaryType() == LocallyRefined ? 2 : 0);

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

    if (pos[i] == 0) {
      local.BoundaryBegin1()[i] = 0;
      local.BoundaryEnd1()[i] = 1;
    }else {
      ++local.SizeTotal()[i];
      local.HaloBegin1()[i] = 0;
      local.HaloEnd1()[i] = 1;
    }

    if (pos[i] == procs[i]-1) {
      local.BoundaryBegin2()[i] = local.SizeTotal()[i] - 1;
      local.BoundaryEnd2()[i] = local.SizeTotal()[i];
    }else {
      ++local.SizeTotal()[i];
      local.HaloBegin2()[i] = local.SizeTotal()[i] - 1;
      local.HaloEnd2()[i] = local.SizeTotal()[i];
    }

  }

  local.Begin() = 1;
  local.End() = local.SizeTotal() - 1;

  local.Size() = local.End() - local.Begin();

  local.AlignmentBegin() = global.AlignmentBegin() + local.Begin();
  local.AlignmentEnd() = global.AlignmentEnd() - local.Begin();

  return local;
}

LocalIndices Multigrid::InitPeriodic(const GlobalIndices& global, const SpatialExtent& extent, const Index& pos, const Index& procs)
{
  LocalIndices local;

  local.SizeTotal() = global.End() - global.Begin() + 2;

  local.BoundaryBegin1() = local.BoundaryEnd1() = 0;
  local.BoundaryBegin2() = local.BoundaryEnd2() = 0;

  local.HaloBegin1() = 0;
  local.HaloEnd1() = 1;
  local.HaloBegin2() = local.SizeTotal() - 1;
  local.HaloEnd2() = local.SizeTotal();

  local.Begin() = 1;
  local.End() = local.SizeTotal() - 1;

  local.Size() = local.End() - local.Begin();

  local.AlignmentBegin() = global.AlignmentBegin() + local.Begin();
  local.AlignmentEnd() = global.AlignmentEnd() + local.Begin();

  return local;
}

void Multigrid::SetLevel(int level)
{
  assert(level >= levelMin && level <= levelMax);
  levelCurrent = level;
  levelIndex = levelMax - level;
}

void Multigrid::ToCoarserLevel()
{
  assert(levelCurrent-1 >= levelMin);
  --levelCurrent;
  ++levelIndex;
}

void Multigrid::ToFinerLevel()
{
  assert(levelCurrent+1 <= levelMax);
  ++levelCurrent;
  --levelIndex;
}

void Multigrid::ClearAll()
{
  for (int i=MinLevel(); i<=MaxLevel(); ++i)
    (*this)(i).Clear();
}

// TODO: diff that in case that breaks QP
void Multigrid::ClearAllCoarseLevels()
{
  for (int i=MinLevel(); i<MaxLevel(); ++i)
    (*this)(i).Clear();
}


void Multigrid::SetCoarserDirichletValues()
{
  Index i_c, i_f;

  for (int i=GlobalMaxLevel(); i>MinLevel(); --i) {

    Grid& grid_f = (*this)(i);
    Grid& grid_c = (*this)(i-1);

    i_c.X() = grid_c.Local().BoundaryBegin1().X();
    i_f.X() = grid_f.Local().BoundaryBegin1().X();
    for (i_c.Y()=grid_c.Local().Begin().Y(); i_c.Y()<grid_c.Local().End().Y(); ++i_c.Y())
      for (i_c.Z()=grid_c.Local().Begin().Z(); i_c.Z()<grid_c.Local().End().Z(); ++i_c.Z())
        grid_c(i_c) = grid_f.GetVal(i_f.X(), 2*i_c.Y(), 2*i_c.Z());

    i_c.X() = grid_c.Local().BoundaryBegin2().X();
    i_f.X() = grid_f.Local().BoundaryBegin2().X();
    for (i_c.Y()=grid_c.Local().Begin().Y(); i_c.Y()<grid_c.Local().End().Y(); ++i_c.Y())
      for (i_c.Z()=grid_c.Local().Begin().Z(); i_c.Z()<grid_c.Local().End().Z(); ++i_c.Z())
        grid_c(i_c) = grid_f.GetVal(i_f.X(), 2*i_c.Y(), 2*i_c.Z());

    i_c.Y() = grid_c.Local().BoundaryBegin1().Y();
    i_f.Y() = grid_f.Local().BoundaryBegin1().Y();
    for (i_c.X()=0; i_c.X()<grid_c.Local().SizeTotal().X(); ++i_c.X())
      for (i_c.Z()=grid_c.Local().Begin().Z(); i_c.Z()<grid_c.Local().End().Z(); ++i_c.Z())
        grid_c(i_c) = grid_f.GetVal(2*i_c.X(), i_f.Y(), 2*i_c.Z());

    i_c.Y() = grid_c.Local().BoundaryBegin2().Y();
    i_f.Y() = grid_f.Local().BoundaryBegin2().Y();
    for (i_c.X()=0; i_c.X()<grid_c.Local().SizeTotal().X(); ++i_c.X())
      for (i_c.Z()=grid_c.Local().Begin().Z(); i_c.Z()<grid_c.Local().End().Z(); ++i_c.Z())
        grid_c(i_c) = grid_f.GetVal(2*i_c.X(), i_f.Y(), 2*i_c.Z());

    i_c.Z() = grid_c.Local().BoundaryBegin1().Z();
    i_f.Z() = grid_f.Local().BoundaryBegin1().Z();
    for (i_c.X()=0; i_c.X()<grid_c.Local().SizeTotal().X(); ++i_c.X())
      for (i_c.Y()=0; i_c.Y()<grid_c.Local().SizeTotal().Y(); ++i_c.Y())
        grid_c(i_c) = grid_f.GetVal(2*i_c.X(), 2*i_c.Y(), i_f.Z());

    i_c.Z() = grid_c.Local().BoundaryBegin2().Z();
    i_f.Z() = grid_f.Local().BoundaryBegin2().Z();
    for (i_c.X()=0; i_c.X()<grid_c.Local().SizeTotal().X(); ++i_c.X())
      for (i_c.Y()=0; i_c.Y()<grid_c.Local().SizeTotal().Y(); ++i_c.Y())
        grid_c(i_c) = grid_f.GetVal(2*i_c.X(), 2*i_c.Y(), i_f.Z());

  }
}

Multigrid::~Multigrid()
{
  for (unsigned int i=0; i<grids.size(); ++i)
    delete grids[i];
}