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

#include "grid/grid_index_translations.hpp"
#include "samples/discretization_poisson_fv.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
{
/*
  //FIXME
  Index i_c, i_f;

  Index b1_c, b2_c, b1_f, b2_f;
  GridIndexTranslations::GetInnerBoundary(sol_c, b1_c, b2_c, sol_f, b1_f, b2_f);
*/
/*
  const Index b1_c = sol_c.Local().FinerBegin();
  const Index b2_c = sol_c.Local().FinerEnd() - 1;
  const Index b1_f = 0;
  const Index b2_f = rhs_f.Local().SizeTotal() - 1;
*/

/*
  const Index begin_f = sol_f.Local().Begin();
  const Index end_f = sol_f.Local().End();
  const Index begin_c = sol_c.Local().FinerBegin();

  const Vector h2_inv = 0.5 / sol_f.Extent().MeshWidth();
  const vmg_float c_1_3 = 1.0 / 3.0;
  const Vector c_2_3_sp = 2.0 / 3.0 * sol_f.Extent().MeshWidth();
  const vmg_float c_4_3 = 4.0 / 3.0;

  Helper::AssertVectorsEqual(sol_c.GetSpatialPos(sol_c.Local().FinerBegin()), sol_f.GetSpatialPos(sol_f.Local().Begin()));
  Helper::AssertVectorsEqual(sol_c.GetSpatialPos(sol_c.Local().FinerEnd()-1), sol_f.GetSpatialPos(sol_f.Local().End()-1));

  //
  // X-direction
  //
  i_f = begin_f; i_c = begin_c;
  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()) {
      Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
      rhs_f(b1_f.X(),i_f.Y(),i_f.Z()) = (sol_c(b1_c.X()-1,i_c.Y(),i_c.Z()) - sol_f(b1_f.X()+1,i_f.Y(),i_f.Z())) * h2_inv.X();
      rhs_f(b2_f.X(),i_f.Y(),i_f.Z()) = (sol_c(b2_c.X()+1,i_c.Y(),i_c.Z()) - sol_f(b2_f.X()-1,i_f.Y(),i_f.Z())) * h2_inv.X();
    }

  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()));
      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(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));
      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(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));

      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(b1_f.X(),i_f.Y(),i_f.Z()) = sol_f(b1_f.X(),i_f.Y(),i_f.Z()) = c_2_3_sp.X() * 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());

      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_sp.X() * 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
  //
  i_f = begin_f; i_c = begin_c;
  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()) {
      Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
      rhs_f(i_f.X(),b1_f.Y(),i_f.Z()) = (sol_c(i_c.X(),b1_c.Y()-1,i_c.Z()) - sol_f(i_f.X(),b1_f.Y()+1,i_f.Z())) * h2_inv.Y();
      rhs_f(i_f.X(),b2_f.Y(),i_f.Z()) = (sol_c(i_c.X(),b2_c.Y()+1,i_c.Z()) - sol_f(i_f.X(),b2_f.Y()-1,i_f.Z())) * h2_inv.Y();
    }

  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()));
      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(),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));
      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(),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));

      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(),b1_f.Y(),i_f.Z()) = sol_f(i_f.X(),b1_f.Y(),i_f.Z()) = c_2_3_sp.Y() * 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());

      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_sp.Y() * 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
  //
  i_f = begin_f; i_c = begin_c;
  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()) {
      Helper::AssertVectorsEqual(sol_c.GetSpatialPos(i_c), sol_f.GetSpatialPos(i_f));
      rhs_f(i_f.X(),i_f.Y(),b1_f.Z()) = (sol_c(i_c.X(),i_c.Y(),b1_c.Z()-1) - sol_f(i_f.X(),i_f.Y(),b1_f.Z()+1)) * h2_inv.Z();
      rhs_f(i_f.X(),i_f.Y(),b2_f.Z()) = (sol_c(i_c.X(),i_c.Y(),b2_c.Z()+1) - sol_f(i_f.X(),i_f.Y(),b2_f.Z()-1)) * h2_inv.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(); 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()));
      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(),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()));
      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(),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()));

      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(),b1_f.Z()) = sol_f(i_f.X(),i_f.Y(),b1_f.Z()) = c_2_3_sp.Z() * 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);

      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_sp.Z() * 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
}