/**
 * @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 "samples/discretization_poisson_fv.hpp"

using namespace VMG;

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

  const vmg_float h2_inv = 0.5 / sol_f.MeshWidth();

  const Index b1_c = sol_c.Local().AlignmentBegin();
  const Index b2_c = sol_c.Local().AlignmentEnd() - 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().AlignmentBegin();

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

  //
  // X-direction
  //
  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()) {
      assert(b2_c.X()+1 >= 0 && b2_c.X() < sol_c.Local().SizeTotal().X());
      assert(i_c.Y() >= 0 && i_c.Y() < sol_c.Local().SizeTotal().Y());
      assert(i_c.Z() >= 0 && i_c.Z() < sol_c.Local().SizeTotal().Z());
      assert(b2_f.X()-1 >= 0 && b2_f.X()-1 < sol_f.Local().SizeTotal().X());
      assert(i_f.Y() >= 0 && i_f.Y() < sol_f.Local().SizeTotal().Y());
      assert(i_f.Z() >= 0 && i_f.Z() < sol_f.Local().SizeTotal().Z());
      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;
      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;
    }

  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 * 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 * 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
  //

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

  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 * 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 * 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
  //

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

  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 * 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 * 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
}