/**
 * @file   domain_decomposition_mpi.cpp
 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
 * @date   Mon Jun 27 12:53:50 2011
 *
 * @brief  Computes a domain decomposition which separates
 *         the finest grid equally for all processes.
 *
 */

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

#include "comm/comm.hpp"
#include "comm/domain_decomposition_mpi.hpp"
#include "grid/grid.hpp"
#include "grid/multigrid.hpp"
#include "interface/interface.hpp"

using namespace VMG;

void DomainDecompositionMPI::Compute(Comm* comm, const Interface* interface, std::vector<GlobalIndices>& global)
{
  GlobalIndices global_l;
  Index remainder, procs;
  Index begin_finest, end_finest;
  Index last_procs = comm->GlobalProcs();

  global.clear();

  for (unsigned int i=0; i<interface->Global().size(); ++i) {

    global_l.BeginFinest() = interface->Global()[i].BeginFinest();
    global_l.EndFinest() = interface->Global()[i].EndFinest();
    global_l.SizeFinest() = interface->Global()[i].SizeFinest();
    global_l.SizeGlobal() = interface->Global()[i].SizeGlobal();

    global_l.BoundaryType() = interface->Global()[i].BoundaryType();

    if (IsActive(comm, global_l.SizeGlobal(), procs)) {

      if (i == 0) {

	global_l.SizeLocal() = global_l.SizeGlobal() / procs;
	remainder = global_l.SizeGlobal() % procs;

	for (int j=0; j<3; ++j)
	  if (comm->GlobalPos()[j] < remainder[j])
	    ++(global_l.SizeLocal()[j]);

	global_l.BeginLocal() = comm->GlobalPos() * global_l.SizeLocal();

	for (int j=0; j<3; ++j)
	  if (comm->GlobalPos()[j] >= remainder[j])
	  global_l.BeginLocal()[j] += remainder[j];

	global_l.EndLocal() = global_l.BeginLocal() + global_l.SizeLocal();

      }else {

	for (int j=0; j<3; ++j) {

	  if (procs[j] == last_procs[j]) {

	    begin_finest[j] = (global.back().BeginFinest() +
			       GridIndexTranslations::GlobalToFiner(global.back().BeginLocal(), i-1))[j];
	    end_finest[j] = (global.back().BeginFinest() +
			     GridIndexTranslations::GlobalToFiner(global.back().EndLocal(), i-1))[j];

	    if (global.back().BeginLocal()[j] == 0)
	      begin_finest[j] = global_l.BeginFinest()[j];

	    if (global.back().EndLocal()[j] == global.back().SizeGlobal()[j])
	      end_finest[j] = global_l.EndFinest()[j];

	    begin_finest[j] -= global_l.BeginFinest()[j];
	    end_finest[j] -= global_l.BeginFinest()[j];

	    FineToCoarse(comm, begin_finest[j], end_finest[j], i);

	    global_l.BeginLocal()[j] = begin_finest[j];
	    global_l.EndLocal()[j] = end_finest[j];
	    global_l.SizeLocal()[j] = end_finest[j] - begin_finest[j];

	  }else {

	    global_l.SizeLocal()[j] = global_l.SizeGlobal()[j] / procs[j];
	    remainder[j] = global_l.SizeGlobal()[j] % procs[j];


	    if (comm->GlobalPos()[j] < remainder[j])
	      ++(global_l.SizeLocal()[j]);

	    global_l.BeginLocal()[j] = comm->GlobalPos()[j] * global_l.SizeLocal()[j];


	    if (comm->GlobalPos()[j] >= remainder[j])
	      global_l.BeginLocal()[j] += remainder[j];

	    global_l.EndLocal()[j] = global_l.BeginLocal()[j] + global_l.SizeLocal()[j];

	  }
	}
      }
    }else {

      global_l.BeginLocal() = 0;
      global_l.EndLocal() = 0;
      global_l.SizeLocal() = 0;

    }

    last_procs = procs;

    global.push_back(global_l);

  }
}

bool DomainDecompositionMPI::IsActive(Comm* comm, const Index& size_global, Index& procs)
{
  bool is_active = true;
  const int points_min = 5;

  procs = size_global / points_min + 1;

  for (int i=0; i<3; ++i) {
    procs[i] = std::min(procs[i], comm->GlobalProcs()[i]);
    is_active &= comm->GlobalPos()[i] < procs[i];
  }

  return is_active;
}

void DomainDecompositionMPI::FineToCoarse(Comm* comm, int& begin, int& end, int levels)
{
  int last_point = end - 1;

  for (int i=0; i<levels; ++i) {

    if (begin % 2 == 0)
      begin /= 2;
    else
      begin = (begin+1) / 2;

    if (last_point % 2 == 0)
      last_point /= 2;
    else
      last_point = (last_point-1) / 2;

  }

end = last_point + 1;
}