source: src/grid/multigrid.cpp@ facba0

/**
 * @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 <iostream>

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

using namespace VMG;

Multigrid::Multigrid(Comm* comm, const Interface* interface)
{
  Index points, remainder;
  LocalIndices local_l;
  std::vector<GlobalIndices> global_separated;

  const std::vector<GlobalIndices>& global = interface->Global();
  const std::vector<SpatialExtent>& extent = interface->Extent();

  comm->GetDomainDecomposition().Compute(comm, interface, global_separated);

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

    /*
     * Check if this level is the finest global level
     */
    if (global[i].BoundaryType() == GlobalMax)
      levelGlobalMax = interface->MaxLevel() - i;


    if (global_separated[i].SizeLocal().Product() > 0) {

      /*
       * Initialize some properties with zero
       */
      local_l.HaloBegin1() = local_l.HaloEnd1() = 0;
      local_l.HaloBegin2() = local_l.HaloEnd2() = 0;
      local_l.BoundaryBegin1() = local_l.BoundaryEnd1() = 0;
      local_l.BoundaryBegin2() = local_l.BoundaryEnd2() = 0;

      /*
       * Set boundary dependant properties
       */
      for (int j=0; j<3; ++j) {

        if (comm->BoundaryConditions()[j] == Dirichlet ||
            comm->BoundaryConditions()[j] == Open) {

          local_l.SizeTotal()[j] = global_separated[i].SizeLocal()[j] +
            (global_separated[i].BoundaryType() == LocallyRefined ? 2 : 0);

          /*
           * We have a boundary at the ends of the process grids
           * and halo grid points otherwise
           */
          if (global_separated[i].BeginLocal()[j] == 0) {
            local_l.BoundaryBegin1()[j] = 0;
            local_l.BoundaryEnd1()[j] = 1;
          }else {
            ++local_l.SizeTotal()[j];
            local_l.HaloBegin1()[j] = 0;
            local_l.HaloEnd1()[j] = 1;
          }

          if (global_separated[i].EndLocal()[j] == global_separated[i].SizeGlobal()[j]) {
            local_l.BoundaryBegin2()[j] = local_l.SizeTotal()[j] - 1;
            local_l.BoundaryEnd2()[j] = local_l.SizeTotal()[j];
          }else {
            ++local_l.SizeTotal()[j];
            local_l.HaloBegin2()[j] = local_l.SizeTotal()[j] - 1;
            local_l.HaloEnd2()[j] = local_l.SizeTotal()[j];
          }

        }else if (comm->BoundaryConditions()[j] == Periodic) {

          local_l.SizeTotal()[j] = global_separated[i].SizeLocal()[j] + 2;

          /*
           * No boundary
           */
          local_l.BoundaryBegin1()[j] = local_l.BoundaryEnd1()[j] = 0;
          local_l.BoundaryBegin2()[j] = local_l.BoundaryEnd2()[j] = 0;

          /*
           * Halo grid points on all processes
           */
          local_l.HaloBegin1()[j] = 0;
          local_l.HaloEnd1()[j] = 1;
          local_l.HaloBegin2()[j] = local_l.SizeTotal()[j] - 1;
          local_l.HaloEnd2()[j] = local_l.SizeTotal()[j];

        }

      }

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

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

      if (i==0)
        local_l.AlignmentBegin() = (global_separated[i].SizeGlobal() - static_cast<int>( pow(2.0, interface->MaxLevel() ) + 0.5 ) - 1) / 2;
      else if (global_separated[i].BoundaryType() == LocallyRefined)
        local_l.AlignmentBegin() = (extent[i].Size() - extent[i-1].Size()) * pow(2.0, interface->MaxLevel()-i-1) + 0.5;
      else if (global_separated[i].BoundaryType() == GlobalMax)
        local_l.AlignmentBegin() = (extent[i].Size() - extent[i-1].Size()) * pow(2.0, interface->MaxLevel()-i-1) - 0.5;
      else
        local_l.AlignmentBegin() = 0;

      local_l.AlignmentBegin() += local_l.Begin();
      local_l.AlignmentEnd() = local_l.AlignmentBegin() + global_separated[i].SizeLocal();

      local_l.AlignmentBegin().Clamp(local_l.Begin(), local_l.End());
      local_l.AlignmentEnd().Clamp(local_l.Begin(), local_l.End());


      grids.push_back(new Grid(global_separated[i], local_l, extent[i], interface->MaxLevel()-i, this));

    }else {

      global_separated[i].SizeGlobal() = global[i].SizeGlobal();
      global_separated[i].BeginFinest() = global[i].BeginFinest();
      global_separated[i].EndFinest() = global[i].EndFinest();
      global_separated[i].SizeFinest() = global[i].SizeFinest();
      global_separated[i].SizeLocal() = 0;
      global_separated[i].BeginLocal() = 0;
      global_separated[i].EndLocal() = 0;
      global_separated[i].BoundaryType() = global[i].BoundaryType();

      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_separated[i], local_l, extent[i], interface->MaxLevel()-i, this));

    }
  }

  numLevels = grids.size();

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

  levelIndex = 0;
  levelCurrent = levelMax;

}

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];
}