source: src/grid/grid_index_translations.cpp@ 13b272

/*
 *    vmg - a versatile multigrid solver
 *    Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
 *
 *  vmg is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  vmg is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/**
 * @file   grid_index_translations.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Tue May 17 11:46:37 2011
 *
 * @brief  Class to convert different representations of grid
 *         indices.
 *
 */

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

#include "base/helper.hpp"
#include "comm/comm.hpp"
#include "grid/grid_double_iterator.hpp"
#include "grid/grid_index_translations.hpp"
#include "grid/grid.hpp"
#include "grid/multigrid.hpp"
#include "mg.hpp"

using namespace VMG;

bool GridIndexTranslations::IsGridPointOf(const Grid& grid, const Index& index_finest)
{
  const int max_level = MG::GetSol()->MaxLevel();
  return index_finest[0] % Helper::intpow(2, max_level - grid.Level()) == 0 &&
         index_finest[1] % Helper::intpow(2, max_level - grid.Level()) == 0 &&
         index_finest[2] % Helper::intpow(2, max_level - grid.Level()) == 0;
}

Index GridIndexTranslations::LocalToGlobal(const Grid& grid, const Index& index_local)
{
  return index_local - grid.Local().HaloSize1() + grid.Global().LocalBegin();
}

Index GridIndexTranslations::LocalToGlobalFinest(const Grid& grid, const Index& index_local)
{
  return GlobalToGlobalFinest(grid, LocalToGlobal(grid, index_local));
}

Index GridIndexTranslations::GlobalToLocal(const Grid& grid, const Index& index_global)
{
  return index_global - grid.Global().LocalBegin() + grid.Local().HaloSize1();
}

Index GridIndexTranslations::GlobalToGlobalFinest(const Grid& grid, const Index& index_global)
{
  return Helper::intpow(2, MG::GetSol()->MaxLevel() - grid.Level()) * index_global;
}

Index GridIndexTranslations::GlobalFinestToLocal(const Grid& grid, const Index& index_finest)
{
  return GlobalToLocal(grid, GlobalFinestToGlobal(grid, index_finest));
}

Index GridIndexTranslations::GlobalFinestToGlobal(const Grid& grid, const Index& index_finest)
{
  assert(IsGridPointOf(grid, index_finest));
  return index_finest / Helper::intpow(2, MG::GetSol()->MaxLevel() - grid.Level());
}

void GridIndexTranslations::GlobalCoarseToFine(Index& begin, Index& end)
{
  for (int j=0; j<3; ++j) {
    begin[j] = 2 * begin[j];
    end[j] = 2 * (end[j]-1) + 1;
  }
}

void GridIndexTranslations::GlobalFineToCoarse(Index& begin, Index& end)
{
  for (int j=0; j<3; ++j) {
    begin[j] = Helper::RoundUpToNextMultiple(begin[j], 2) / 2;
    end[j] = Helper::RoundDownToNextMultiple(end[j]-1, 2) / 2 + 1;
  }
}

void GridIndexTranslations::GetGridAlignment(const Grid& grid_1, GridIteratorSet& bounds_1,
                                             const Grid& grid_2, GridIteratorSet& bounds_2)
{
  const Boundary& boundary = MG::GetComm()->BoundaryConditions();

  if (grid_1.Level() == grid_2.Level()) {
    Index begin_global = grid_1.Global()
                               .LocalBegin()
                               .Clamp(grid_2.Global().LocalBegin(), grid_2.Global().LocalEnd());

    Index end_global = grid_1.Global()
                             .LocalEnd()
                             .Clamp(grid_2.Global().LocalBegin(), grid_2.Global().LocalEnd());

    for (int j=0; j<3; ++j) {
      if (boundary[j] == Dirichlet || boundary[j] == Open) {
        if (begin_global[j] == grid_1.Global().GlobalBegin()[j] || begin_global[j] == grid_2.Global().GlobalBegin()[j])
          begin_global[j] += 1;
        if (end_global[j] == grid_1.Global().GlobalEnd()[j] || end_global[j] == grid_2.Global().GlobalEnd()[j])
          end_global[j] -= 1;
      }
    }

    bounds_1 = GridIteratorSet(GlobalToLocal(grid_1, begin_global), GlobalToLocal(grid_2, end_global));
    bounds_2 = GridIteratorSet(GlobalToLocal(grid_2, begin_global), GlobalToLocal(grid_2, end_global));

  } else {

    const Grid& grid_c = (grid_1.Level() < grid_2.Level() ? grid_1 : grid_2);
    const int global_mult = Helper::intpow(2, MG::GetSol()->MaxLevel() - grid_c.Level());

    Index begin_finest = GlobalToGlobalFinest(grid_1, grid_1.Global().LocalBegin())
                           .Clamp(GlobalToGlobalFinest(grid_2, grid_2.Global().LocalBegin()),
                                  GlobalToGlobalFinest(grid_2, grid_2.Global().LocalEnd()-1));

    Index end_finest = GlobalToGlobalFinest(grid_1, grid_1.Global().LocalEnd()-1)
                         .Clamp(GlobalToGlobalFinest(grid_2, grid_2.Global().LocalBegin()),
                                GlobalToGlobalFinest(grid_2, grid_2.Global().LocalEnd()-1));

    for (int j=0; j<3; ++j) {
      begin_finest[j] = Helper::RoundUpToNextMultiple(begin_finest[j], global_mult);
      end_finest[j] = Helper::RoundDownToNextMultiple(end_finest[j], global_mult);
    }

    for (int j=0; j<3; ++j) {
      if (boundary[j] == Dirichlet || boundary[j] == Open) {
        if (grid_1.Global().LocalBegin()[j] == grid_1.Global().GlobalBegin()[j] || grid_2.Global().LocalBegin()[j] == grid_2.Global().GlobalBegin()[j])
          begin_finest[j] += global_mult;
        if (grid_1.Global().LocalEnd()[j] == grid_1.Global().GlobalEnd()[j] || grid_2.Global().LocalEnd()[j] == grid_2.Global().GlobalEnd()[j])
          end_finest[j] -= global_mult;
      }
    }

    bounds_1 = GridIteratorSet(GlobalFinestToLocal(grid_1, begin_finest), GlobalFinestToLocal(grid_1, end_finest)+1);
    bounds_2 = GridIteratorSet(GlobalFinestToLocal(grid_2, begin_finest), GlobalFinestToLocal(grid_2, end_finest)+1);

  }
}

void GridIndexTranslations::GetInnerBoundary(const Grid& grid_c, Index& b1_c, Index& b2_c,
                                             const Grid& grid_f, Index& b1_f, Index& b2_f)
{
  if (grid_f.Global().BoundaryType() == LocallyRefined) {

    const Boundary& bc = MG::GetComm()->BoundaryConditions();

    const int global_mult = Helper::intpow(2, MG::GetSol()->MaxLevel() - grid_c.Level());

    for (int i=0; i<3; ++i) {
      if (bc[i] == Open) {
        int begin_f_global = grid_f.Global().LocalBegin()[i] + (grid_f.Global().LocalBegin()[i] % 2 == 0 ? 0 : 1);
        int end_f_global = grid_f.Global().LocalEnd()[i] + (grid_f.Global().LocalEnd()[i] % 2 == 0 ? 1 : 0);
      } else {

      }
    }

    Index begin_finest = GlobalToGlobalFinest(grid_c, grid_c.Global().LocalBegin())
                           .Clamp(GlobalToGlobalFinest(grid_f, grid_f.Global().LocalBegin()),
                                  GlobalToGlobalFinest(grid_f, grid_f.Global().LocalEnd()-1));

    Index end_finest = GlobalToGlobalFinest(grid_c, grid_c.Global().LocalEnd()-1)
                         .Clamp(GlobalToGlobalFinest(grid_f, grid_f.Global().LocalBegin()),
                                GlobalToGlobalFinest(grid_f, grid_f.Global().LocalEnd()-1));

  } else {
    b1_c = -1;
    b2_c = -1;
    b1_f = -1;
    b2_f = -1;
  }
}