/*
 *    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   discretization_poisson_fv.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Apr 18 13:03:47 2011
 *
 * @brief  Finite volume discretization for the Poisson
 *         equation. Absolutely equivalent to the finite
 *         difference discretization unless you use
 *         hierarchically coarsened grids.
 *
 */

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

#include "comm/comm.hpp"
#include "discretization/discretization_poisson_fv.hpp"
#include "grid/grid_index_translations.hpp"

using namespace VMG;

void DiscretizationPoissonFV::InitDiscretizationPoissonFV()
{
  switch (order)
  {
  case 2:
    stencil.SetDiag(6.0);
    stencil.push_back(-1,  0,  0, -1.0);
    stencil.push_back( 1,  0,  0, -1.0);
    stencil.push_back( 0, -1,  0, -1.0);
    stencil.push_back( 0,  1,  0, -1.0);
    stencil.push_back( 0,  0, -1, -1.0);
    stencil.push_back( 0,  0,  1, -1.0);
    break;
  case 4:
    stencil.SetDiag(24.0/6.0);
    stencil.push_back(-1,  0,  0, -2.0/6.0);
    stencil.push_back( 1,  0,  0, -2.0/6.0);
    stencil.push_back( 0, -1,  0, -2.0/6.0);
    stencil.push_back( 0,  1,  0, -2.0/6.0);
    stencil.push_back( 0,  0, -1, -2.0/6.0);
    stencil.push_back( 0,  0,  1, -2.0/6.0);
    stencil.push_back(-1, -1,  0, -1.0/6.0);
    stencil.push_back(-1,  1,  0, -1.0/6.0);
    stencil.push_back( 1, -1,  0, -1.0/6.0);
    stencil.push_back( 1,  1,  0, -1.0/6.0);
    stencil.push_back(-1,  0, -1, -1.0/6.0);
    stencil.push_back(-1,  0,  1, -1.0/6.0);
    stencil.push_back( 1,  0, -1, -1.0/6.0);
    stencil.push_back( 1,  0,  1, -1.0/6.0);
    stencil.push_back( 0, -1, -1, -1.0/6.0);
    stencil.push_back( 0, -1,  1, -1.0/6.0);
    stencil.push_back( 0,  1, -1, -1.0/6.0);
    stencil.push_back( 0,  1,  1, -1.0/6.0);
    break;
  default:
    assert(0 != "vmg choose discretization order 2 or 4");
    break;
  }
}

void DiscretizationPoissonFV::ModifyRightHandSide()
{
  if (order == 4) {

    Grid& rhs = MG::GetRhsMaxLevel();

    Stencil stencil(6.0/12.0);
    stencil.push_back(-1,  0,  0, 1.0/12.0);
    stencil.push_back( 1,  0,  0, 1.0/12.0);
    stencil.push_back( 0, -1,  0, 1.0/12.0);
    stencil.push_back( 0,  1,  0, 1.0/12.0);
    stencil.push_back( 0,  0, -1, 1.0/12.0);
    stencil.push_back( 0,  0,  1, 1.0/12.0);

    stencil.Apply(rhs);

  }
}

void DiscretizationPoissonFV::SetInnerBoundaryCompute(Grid& sol_f, Grid& rhs_f, Grid& sol_c) const
{
  Index i_c, i_f;

  Comm& comm = *MG::GetComm();

  const Boundary& bc = comm.BoundaryConditions();
  const Index off((GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[0] % 2 == 0 ? 0 : 1),
                  (GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[1] % 2 == 0 ? 0 : 1),
                  (GridIndexTranslations::LocalToGlobal(sol_f, sol_f.Local().Begin())[2] % 2 == 0 ? 0 : 1));

  const Index begin_f = sol_f.Local().Begin() - off;
  const Index end_f = sol_f.Local().End();

  const Index begin_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, GridIndexTranslations::LocalToGlobalFinest(sol_f, begin_f));

  const Index b1_f = sol_f.Local().BoundaryBegin1();
  const Index b2_f = sol_f.Local().BoundaryBegin2();

  const Index b1_finest = GridIndexTranslations::LocalToGlobalFinest(sol_f, sol_f.Local().Begin() + off);
  const Index b2_finest = GridIndexTranslations::LocalToGlobalFinest(sol_f, sol_f.Local().End()-1);

  const Index b1_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, b1_finest) - 1;
  const Index b2_c = GridIndexTranslations::GlobalFinestToLocal(sol_c, b2_finest) + 1;

  const vmg_float c_1_3 = 1.0 / 3.0;
  const vmg_float c_2_3 = 2.0 / 3.0;
  const vmg_float c_4_3 = 4.0 / 3.0;

  comm.CommToGhosts(sol_f);
  comm.CommToGhosts(sol_c);

  //
  // X-direction
  //
  if (bc.X() == Open) {

    if (sol_f.Local().BoundarySize1().X() > 0) {

      for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
        for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (sol_c(b1_c.X(),i_c.Y(),i_c.Z()) - sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()));

      for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()) + rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()));

      for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b1_f.X(),i_f.Y(),i_f.Z()-1) + rhs_f(b1_f.X(),i_f.Y(),i_f.Z()+1));

      for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = 0.25 * (rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()-1) +
                                                    rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()-1) +
                                                    rhs_f(b1_f.X(),i_f.Y()-1,i_f.Z()+1) +
                                                    rhs_f(b1_f.X(),i_f.Y()+1,i_f.Z()+1));

      for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
          rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = sol_f(b1_f.X(),i_f.Y(),i_f.Z()) = c_2_3 * rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) +
                                                                              c_4_3 * sol_f(b1_f.X()+1,i_f.Y(),i_f.Z()) -
                                                                              c_1_3 * sol_f(b1_f.X()+2,i_f.Y(),i_f.Z());

    }

    if (sol_f.Local().BoundarySize2().X() > 0) {

      for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
        for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (sol_c(b2_c.X(),i_c.Y(),i_c.Z()) - sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()));

      for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()) + rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()));

      for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.5 * (rhs_f(b2_f.X(),i_f.Y(),i_f.Z()-1) + rhs_f(b2_f.X(),i_f.Y(),i_f.Z()+1));

      for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = 0.25 * (rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()-1) +
                                                    rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()-1) +
                                                    rhs_f(b2_f.X(),i_f.Y()-1,i_f.Z()+1) +
                                                    rhs_f(b2_f.X(),i_f.Y()+1,i_f.Z()+1));

      for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
          rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = sol_f(b2_f.X(),i_f.Y(),i_f.Z()) = c_2_3 * rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) +
                                                                              c_4_3 * sol_f(b2_f.X()-1,i_f.Y(),i_f.Z()) -
                                                                              c_1_3 * sol_f(b2_f.X()-2,i_f.Y(),i_f.Z());

    }
  }

  //
  // Y-direction
  //
  if (bc.Y() == Open) {

    if (sol_f.Local().BoundarySize1().Y() > 0) {

      for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
        for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (sol_c(i_c.X(),b1_c.Y(),i_c.Z()) - sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()));

      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()));

      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X(),b1_f.Y(),i_f.Z()-1) + rhs_f(i_f.X(),b1_f.Y(),i_f.Z()+1));

      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()-1) +
                                                    rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()-1) +
                                                    rhs_f(i_f.X()-1,b1_f.Y(),i_f.Z()+1) +
                                                    rhs_f(i_f.X()+1,b1_f.Y(),i_f.Z()+1));

      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
          rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = sol_f(i_f.X(),b1_f.Y(),i_f.Z()) = c_2_3 * rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) +
                                                                              c_4_3 * sol_f(i_f.X(),b1_f.Y()+1,i_f.Z()) -
                                                                              c_1_3 * sol_f(i_f.X(),b1_f.Y()+2,i_f.Z());

    }

    if (sol_f.Local().BoundarySize2().Y() > 0) {

      for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
        for (i_f.Z()=begin_f.Z(), i_c.Z()=begin_c.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2, ++i_c.Z())
          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (sol_c(i_c.X(),b2_c.Y(),i_c.Z()) - sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()));

      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()) + rhs_f(i_f.X()+1,b2_f.Z(),i_f.Z()));

      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.5 * (rhs_f(i_f.X(),b2_f.Y(),i_f.Z()-1) + rhs_f(i_f.X(),b2_f.Y(),i_f.Z()+1));

      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Z()=begin_f.Z()+1; i_f.Z()<end_f.Z(); i_f.Z()+=2)
          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()-1) +
                                                    rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()-1) +
                                                    rhs_f(i_f.X()-1,b2_f.Y(),i_f.Z()+1) +
                                                    rhs_f(i_f.X()+1,b2_f.Y(),i_f.Z()+1));

      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
        for (i_f.Z()=begin_f.Z(); i_f.Z()<end_f.Z(); ++i_f.Z())
          rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = sol_f(i_f.X(),b2_f.Y(),i_f.Z()) = c_2_3 * rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) +
                                                                              c_4_3 * sol_f(i_f.X(),b2_f.Y()-1,i_f.Z()) -
                                                                              c_1_3 * sol_f(i_f.X(),b2_f.Y()-2,i_f.Z());

    }

  }

  //
  // Z-direction
  //
  if (bc.Z() == Open) {

    if (sol_f.Local().BoundarySize1().Z() > 0) {

      for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
        for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (sol_c(i_c.X(),i_c.Y(),b1_c.Z()) - sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1));

      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,i_f.Y(),b1_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b1_f.Z()));

      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.5 * (rhs_f(i_f.X(),i_f.Y()-1,b1_f.Z()) + rhs_f(i_f.X(),i_f.Y()+1,b1_f.Z()));

      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,i_f.Y()-1,b1_f.Z()) +
                                                    rhs_f(i_f.X()+1,i_f.Y()-1,b1_f.Z()) +
                                                    rhs_f(i_f.X()-1,i_f.Y()+1,b1_f.Z()) +
                                                    rhs_f(i_f.X()+1,i_f.Y()+1,b1_f.Z()));

      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
        for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
          rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = sol_f(i_f.X(),i_f.Y(),b1_f.Z()) = c_2_3 * rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) +
                                                                              c_4_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1) -
                                                                              c_1_3 * sol_f(i_f.X(),i_f.Y(),b1_f.Z()+2);

    }

    if (sol_f.Local().BoundarySize2().Z() > 0) {

      for (i_f.X()=begin_f.X(), i_c.X()=begin_c.X(); i_f.X()<end_f.X(); i_f.X()+=2, ++i_c.X())
        for (i_f.Y()=begin_f.Y(), i_c.Y()=begin_c.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2, ++i_c.Y())
          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (sol_c(i_c.X(),i_c.Y(),b2_c.Z()) - sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1));

      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); i_f.Y()+=2)
          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (rhs_f(i_f.X()-1,i_f.Y(),b2_f.Z()) + rhs_f(i_f.X()+1,i_f.Y(),b2_f.Z()));

      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.5 * (rhs_f(i_f.X(),i_f.Y()-1,b2_f.Z()) + rhs_f(i_f.X(),i_f.Y()+1,b2_f.Z()));

      for (i_f.X()=begin_f.X()+1; i_f.X()<end_f.X(); i_f.X()+=2)
        for (i_f.Y()=begin_f.Y()+1; i_f.Y()<end_f.Y(); i_f.Y()+=2)
          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = 0.25 * (rhs_f(i_f.X()-1,i_f.Y()-1,b2_f.Z()) +
                                                    rhs_f(i_f.X()+1,i_f.Y()-1,b2_f.Z()) +
                                                    rhs_f(i_f.X()-1,i_f.Y()+1,b2_f.Z()) +
                                                    rhs_f(i_f.X()+1,i_f.Y()+1,b2_f.Z()));

      for (i_f.X()=begin_f.X(); i_f.X()<end_f.X(); ++i_f.X())
        for (i_f.Y()=begin_f.Y(); i_f.Y()<end_f.Y(); ++i_f.Y())
          rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = sol_f(i_f.X(),i_f.Y(),b2_f.Z()) = c_2_3 * rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) +
                                                                              c_4_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1) -
                                                                              c_1_3 * sol_f(i_f.X(),i_f.Y(),b2_f.Z()-2);

    }

  }

#ifdef DEBUG_MATRIX_CHECKS
  rhs_f.IsConsistent();
  sol_f.IsConsistent();
#endif
}