source: src/base/interface.cpp@ f57182

/*
 *    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   interface.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 12:55:48 2011
 *
 * @brief  VMG::Interface
 *
 */

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

#include <algorithm>
#include <cmath>

#include "base/helper.hpp"
#include "base/interface.hpp"

using namespace VMG;

static Index GetGlobalIndex(const Vector& pos, const SpatialExtent& extent, const BT& bt)
{
  const Index index = (pos - extent.Begin()) / extent.MeshWidth() + 0.5;
  return index + (bt == LocallyRefined ? 1 : 0);
}

void Interface::InitInterface(const Vector& box_offset, const vmg_float& box_size,
                              const int& coarseningSteps, const vmg_float& alpha)
{
  int i;
  Index size_factor;

  const Vector box_center = box_offset + 0.5 * box_size;

  /*
   * Get Extents
   */
  for (i=0; i<coarseningSteps; ++i) {
    extent.push_back(SpatialExtent());
    for (int j=0; j<3; ++j)
      size_factor[j] = (bc[j] == Periodic ? 1 : Helper::intpow(2, static_cast<int>(log(pow(alpha, i+1)) / log(2.0) + 1.0)));

    extent.back().Size() = box_size * static_cast<Vector>(size_factor);
    extent.back().Begin() = box_center - 0.5 * extent.back().Size();
    extent.back().End() = extent.back().Begin() + extent.back().Size();
    extent.back().MeshWidth() = pow(2.0, i-levelMax) * extent.back().Size() / size_factor;
  }

  if (extent.size() == 0) {
    extent.push_back(SpatialExtent());
    extent.back().Size() = box_size;
    extent.back().Begin() = box_offset;
    extent.back().End() = box_offset + box_size;
    extent.back().MeshWidth() = box_size / Helper::intpow(2,levelMax);
  }

  while ((extent.back().Size() / extent.back().MeshWidth()).Min() > Helper::intpow(2,levelMin) + 1) {
    extent.push_back(SpatialExtent(extent.back()));
    extent.back().MeshWidth() *= 2;
  }

  /*
   * Find GlobalMax
   */
  for (std::vector<SpatialExtent>::const_iterator iter=extent.begin(); iter!=extent.end(); ++iter) {
    global.push_back(GlobalIndices());
    global.back().BoundaryType() = GlobalCoarsened;
  }

  for (i=global.size()-1; i>=0; --i) {
    if (i == 0 || extent[i-1].Size() != extent[i].Size()) {
      global[i].BoundaryType() = GlobalMax;
      break;
    }
  }
  for (--i; i>=0; --i)
    global[i].BoundaryType() = LocallyRefined;

  /*
   * Compute global grid values
   */
  for (i=0; i<global.size(); ++i) {

    const Index num_cells = extent[i].Size() / extent[i].MeshWidth() + 0.5;

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

      if (bc[j] == Dirichlet || bc[j] == Open) {
        /*
         * Dirichlet
         */

        if (global[i].BoundaryType() == LocallyRefined) {
          /*
           * Locally Refined
           */

          global[i].GlobalSize()[j] = num_cells[j] + 3;

          if (i == 0) {
            global[i].GlobalFinerBegin()[j] = GetGlobalIndex(box_offset, extent[i], global[i].BoundaryType())[j];
            global[i].GlobalFinerEnd()[j] = GetGlobalIndex(box_offset + box_size, extent[i], global[i].BoundaryType())[j] + 1;
            global[i].GlobalFinerSize()[j] = global[i].GlobalFinerEnd()[j] - global[i].GlobalFinerBegin()[j];
          } else {
            global[i].GlobalFinerBegin()[j] = GetGlobalIndex(extent[i-1].Begin(), extent[i], global[i].BoundaryType())[j];
            global[i].GlobalFinerEnd()[j] = GetGlobalIndex(extent[i-1].End(), extent[i], global[i].BoundaryType())[j] + 1;
            global[i].GlobalFinerSize()[j] = global[i].GlobalFinerEnd()[j] - global[i].GlobalFinerBegin()[j];
          }

        } else if (global[i].BoundaryType() == GlobalMax) {
          /*
           * Global Max
           */

          global[i].GlobalSize()[j] = num_cells[j] + 1;

          if (i == 0) {
            global[i].GlobalFinerBegin()[j] = 0;
            global[i].GlobalFinerEnd()[j] = num_cells[j] + 1;
            global[i].GlobalFinerSize()[j] = num_cells[j] + 1;
          } else {
            global[i].GlobalFinerBegin()[j] = GetGlobalIndex(extent[i-1].Begin(), extent[i], global[i].BoundaryType())[j];
            global[i].GlobalFinerEnd()[j] = GetGlobalIndex(extent[i-1].End(), extent[i], global[i].BoundaryType())[j] + 1;
            global[i].GlobalFinerSize()[j] = global[i].GlobalFinerEnd()[j] - global[i].GlobalFinerBegin()[j];
          }

        } else {
          /*
           * Global Coarsened
           */

          global[i].GlobalSize()[j] = num_cells[j] + 1;

          global[i].GlobalFinerBegin()[j] = 0;
          global[i].GlobalFinerEnd()[j] = num_cells[j] + 1;
          global[i].GlobalFinerSize()[j] = num_cells[j] + 1;

        }

      } else if (bc[j] == Periodic) {
        /*
         * Periodic
         */

        global[i].GlobalSize()[j] = num_cells[j];

        global[i].GlobalFinerBegin()[j] = 0;
        global[i].GlobalFinerEnd()[j] = num_cells[j];
        global[i].GlobalFinerSize()[j] = num_cells[j];

      }
    }

    global[i].LocalBegin() = 0;
    global[i].LocalEnd() = global[i].GlobalSize();
    global[i].LocalSize() = global[i].GlobalSize();

    global[i].LocalFinerBegin() = global[i].GlobalFinerBegin();
    global[i].LocalFinerEnd() = global[i].GlobalFinerEnd();
    global[i].LocalFinerSize() = global[i].GlobalFinerSize();

  }

  levelMin = levelMax - global.size() + 1;

}