source: ThirdParty/vmg/src/discretization/boundary_value_setter.cpp@ e9f307

/*
 *    vmg - a versatile multigrid solver
 *    Copyright (C) 2012-2013 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   boundary_value_setter.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 15 16:21:40 2013
 *
 * @brief  Class for setting boundary values with open
 *         boundary conditions.
 *
 */

#ifdef HAVE_CONFIG_H
#include "libvmg_config.h"
#endif

#include "comm/comm.hpp"
#include "discretization/boundary_value_setter.hpp"
#include "grid/grid.hpp"
#include "grid/grid_index_translations.hpp"
#include "grid/multigrid.hpp"
#include "mg.hpp"

using namespace VMG;

void BoundaryValueSetter::ComputeBoundaryValues()
{
  const Grid& grid_ml = (*MG::GetRhs())(MG::GetRhs()->MaxLevel());
  Grid& rhs_gml = (*MG::GetRhs())(MG::GetRhs()->GlobalMaxLevel());
  Grid& sol_gml = (*MG::GetSol())(MG::GetSol()->GlobalMaxLevel());

  for (std::vector<BoundaryValue>::iterator iter_b = boundary_values.begin(); iter_b != boundary_values.end(); ++iter_b)
    iter_b->Val() = 0.0;

  for (GridIterator iter_g = grid_ml.Iterators().Local().Begin(); iter_g != grid_ml.Iterators().Local().End(); ++iter_g)
    for (std::vector<BoundaryValue>::iterator iter_b = boundary_values.begin(); iter_b != boundary_values.end(); ++iter_b)
      iter_b->Val() += grid_ml.GetVal(*iter_g) * greens_function((iter_b->GetPos() - grid_ml.GetSpatialPos(*iter_g)).Length());

  std::vector<BoundaryValue> list_bv = MG::GetComm()->CommBoundaryValues();

  for (std::vector<BoundaryValue>::iterator iter = list_bv.begin(); iter != list_bv.end(); ++iter) {
    const Index i = GridIndexTranslations::GlobalToLocal(rhs_gml, iter->GetIndex());
    const vmg_float val = grid_ml.Extent().MeshWidth().Product() * iter->Val();
    rhs_gml(i) = val;
    sol_gml(i) = val;
  }

  // MG::GetComm()->PrintGrid(rhs_gml, "DEBUG GREEN BOUNDARY");
}

void BoundaryValueSetter::InitBoundaryValueSetter()
{
  Comm& comm = *MG::GetComm();

  const Grid& grid = (*MG::GetRhs())(MG::GetRhs()->GlobalMaxLevel());

  Index begin_local, end_local, i;
  for (int i=0; i<3; ++i) {
    begin_local[i] = grid.Global().GlobalBegin()[i] + (comm.BoundaryConditions()[i] == Open ? 1 : 0);
    end_local[i] = grid.Global().GlobalEnd()[i] - (comm.BoundaryConditions()[i] == Open ? 1 : 0);
  }

  if (comm.BoundaryConditions().X() == Open) {

    i.X() = grid.Global().GlobalBegin().X();
    for (i.Y() = begin_local.Y(); i.Y() < end_local.Y(); ++i.Y())
      for (i.Z() = begin_local.Z(); i.Z() < end_local.Z(); ++i.Z())
        boundary_values.push_back(BoundaryValue(grid, i));

    i.X() = grid.Global().GlobalEnd().X()-1;
    for (i.Y() = begin_local.Y(); i.Y() < end_local.Y(); ++i.Y())
      for (i.Z() = begin_local.Z(); i.Z() < end_local.Z(); ++i.Z())
        boundary_values.push_back(BoundaryValue(grid, i));

  }

  if (comm.BoundaryConditions().Y() == Open) {

    i.Y() = grid.Global().GlobalBegin().Y();
    for (i.X() = begin_local.X(); i.X() < end_local.X(); ++i.X())
      for (i.Z() = begin_local.Z(); i.Z() < end_local.Z(); ++i.Z())
        boundary_values.push_back(BoundaryValue(grid, i));

    i.Y() = grid.Global().GlobalEnd().Y() - 1;
    for (i.X() = begin_local.X(); i.X() < end_local.X(); ++i.X())
      for (i.Z() = begin_local.Z(); i.Z() < end_local.Z(); ++i.Z())
        boundary_values.push_back(BoundaryValue(grid, i));

  }

  if (comm.BoundaryConditions().Z() == Open) {

    i.Z() = grid.Global().GlobalBegin().Z();
    for (i.X() = begin_local.X(); i.X() < end_local.X(); ++i.X())
      for (i.Y() = begin_local.Y(); i.Y() < end_local.Y(); ++i.Y())
        boundary_values.push_back(BoundaryValue(grid, i));

    i.Z() = grid.Global().GlobalEnd().Z() - 1;
    for (i.X() = begin_local.X(); i.X() < end_local.X(); ++i.X())
      for (i.Y() = begin_local.Y(); i.Y() < end_local.Y(); ++i.Y())
        boundary_values.push_back(BoundaryValue(grid, i));

  }

}

BoundaryValueSetter::~BoundaryValueSetter()
{

}