source: ThirdParty/vmg/src/base/interface.cpp@ 10b1ef

/*
 *    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 <libvmg_config.h>
#endif

#include <algorithm>
#include <cmath>

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

using namespace VMG;

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

  const Index add_node = Index(bc[0]==Periodic ? 0 : 1,
                               bc[1]==Periodic ? 0 : 1,
                               bc[2]==Periodic ? 0 : 1);

  const Index inner_boundary = Index(bc[0]==Open ? 2 : 0,
                                     bc[1]==Open ? 2 : 0,
                                     bc[2]==Open ? 2 : 0);


  const Vector box_center = box_offset + 0.5 * box_size;

  /*
   * Get Extents
   */
   if (bc[0] == Open || bc[1] == Open || bc[2] == Open) {

     //TODO: Change this to max_boundary_nodes at one point
     while (global.size() == 0 ||
         (bc[0] == Open && global.back().GlobalSize()[0] > Helper::intpow(2, levelMin)+1) ||
         (bc[1] == Open && global.back().GlobalSize()[1] > Helper::intpow(2, levelMin)+1) ||
         (bc[2] == Open && global.back().GlobalSize()[2] > Helper::intpow(2, levelMin)+1)) {

       global.push_back(GlobalIndices());
       extent.push_back(SpatialExtent());

       for (int j=0; j<3; ++j)
         size_factor[j] = (bc[j] == Open ? Helper::intpow(2, static_cast<int>(log(pow(alpha, global.size())) / log(2.0) + 1.0)) : 1);

       num_cells = static_cast<Vector>(std::pow(2.0, levelMax-static_cast<int>(global.size())+1)) * size_factor + 0.5;

       extent.back().MeshWidth() = box_size * static_cast<Vector>(size_factor) / num_cells;
       extent.back().Size() = (num_cells + inner_boundary) * extent.back().MeshWidth();
       extent.back().Begin() = box_center - 0.5 * extent.back().Size();
       extent.back().End() = extent.back().Begin() + extent.back().Size();

       global.back().LocalSize() = num_cells + add_node + inner_boundary;
       global.back().LocalBegin() = -1 * (num_cells + inner_boundary) / 2;
       global.back().LocalEnd() = (num_cells + inner_boundary) / 2 + add_node;

       global.back().GlobalSize() = global.back().LocalSize();
       global.back().GlobalBegin() = global.back().LocalBegin();
       global.back().GlobalEnd() = global.back().LocalEnd();

       global.back().GlobalSizeFinest() = Helper::intpow(2, global.size()-1) * (num_cells+inner_boundary) + add_node;
       global.back().GlobalBeginFinest() = -1 * ((Helper::intpow(2, global.size()-1) * (num_cells + inner_boundary)) / 2);
       global.back().GlobalEndFinest() = (Helper::intpow(2, global.size()-1) * (num_cells + inner_boundary)) / 2 + add_node;

       global.back().BoundaryType() = LocallyRefined;

     }

     global.back().BoundaryType() = GlobalMax;

   } else {

     num_cells = Helper::intpow(2, levelMax);

     global.push_back(GlobalIndices());
     extent.push_back(SpatialExtent());

     extent.back().Size() = box_size;
     extent.back().Begin() = box_center - 0.5 * extent.back().Size();
     extent.back().End() = extent.back().Begin() + extent.back().Size();
     extent.back().MeshWidth() = extent.back().Size() / num_cells;

     global.back().LocalSize() = num_cells + add_node;
     global.back().LocalBegin() = -1 * num_cells/2;
     global.back().LocalEnd() = num_cells/2 + add_node;

     global.back().GlobalSize() = global.back().LocalSize();
     global.back().GlobalBegin() = global.back().LocalBegin();
     global.back().GlobalEnd() = global.back().LocalEnd();

     global.back().GlobalSizeFinest() = global.back().LocalSize();
     global.back().GlobalBeginFinest() = global.back().LocalBegin();
     global.back().GlobalEndFinest() = global.back().LocalEnd();

     global.back().BoundaryType() = GlobalMax;

   }

   while (global.back().GlobalSize().Min() > Helper::intpow(2, levelMin)+1) {

     num_cells /= 2;

     global.push_back(GlobalIndices());
     extent.push_back(SpatialExtent());

     extent.back().Size() = (++extent.rbegin())->Size();
     extent.back().Begin() = (++extent.rbegin())->Begin();
     extent.back().End() = (++extent.rbegin())->End();
     extent.back().MeshWidth() = 2.0 * (++extent.rbegin())->MeshWidth();

     global.back().LocalSize() = num_cells + add_node;
     global.back().LocalBegin() = -1 * num_cells/2;
     global.back().LocalEnd() = num_cells/2 + add_node;

     global.back().GlobalSize() = global.back().LocalSize();
     global.back().GlobalBegin() = global.back().LocalBegin();
     global.back().GlobalEnd() = global.back().LocalEnd();

     global.back().GlobalSizeFinest() = (++global.rbegin())->GlobalSizeFinest();
     global.back().GlobalBeginFinest() = (++global.rbegin())->GlobalBeginFinest();
     global.back().GlobalEndFinest() = (++global.rbegin())->GlobalEndFinest();

     global.back().BoundaryType() = GlobalCoarsened;

   }

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