#include "ParallelDomainProblem.h"
#include "ProblemScal.h"
#include "ProblemInstat.h"
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
#include "DOFMatrix.h"
#include "DOFVector.h"
#include "VtkWriter.h"

#include "petscksp.h"

namespace AMDiS {

  ParallelDomainProblemBase::ParallelDomainProblemBase(const std::string& name,
						       ProblemIterationInterface *iIF,
						       ProblemTimeInterface *tIF,
						       FiniteElemSpace *fe)
    : iterationIF(iIF),
      timeIF(tIF),
      feSpace(fe),
      mesh(fe->getMesh()),
      initialPartitionMesh(true),
      nRankDOFs(0)
  {
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
  }

  void ParallelDomainProblemBase::initParallelization(AdaptInfo *adaptInfo)
  {
    if (mpiSize <= 1)
      return;

    // create an initial partitioning of the mesh
    partitioner->createPartitionData();
    // set the element weights, which are 1 at the very first begin
    setElemWeights(adaptInfo);
    // and now partition the mesh
    partitionMesh(adaptInfo);   


    // === Create new global and local DOF numbering. ===

    int nRankDOFs = 0;
    int nOverallDOFs = 0;
    createLocalGlobalNumbering(nRankDOFs, nOverallDOFs);


    // === Create interior boundary information ===

    createInteriorBoundaryInfo();


    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
      

    /// === Reset all DOFAdmins of the mesh. ===

    int nAdmins = mesh->getNumberOfDOFAdmin();
    for (int i = 0; i < nAdmins; i++) {
      for (int j = 0; j < mesh->getDOFAdmin(i).getSize(); j++)
 	const_cast<DOFAdmin&>(mesh->getDOFAdmin(i)).setDOFFree(j, true);
      for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
 	const_cast<DOFAdmin&>(mesh->getDOFAdmin(i)).setDOFFree(j, false);
    }


    /// === Create petsc matrix. ===

    int ierr;
    ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
    ierr = MatSetSizes(petscMatrix, nRankDOFs, nRankDOFs, nOverallDOFs, nOverallDOFs);
    ierr = MatSetType(petscMatrix, MATAIJ);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
    ierr = VecSetSizes(petscRhsVec, nRankDOFs, nOverallDOFs);
    ierr = VecSetType(petscRhsVec, VECMPI);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscSolVec);
    ierr = VecSetSizes(petscSolVec, nRankDOFs, nOverallDOFs);
    ierr = VecSetType(petscSolVec, VECMPI);
  }

  void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
  {}


  void ParallelDomainProblemBase::fillPetscMatrix(DOFMatrix *mat, 
						  DOFVector<double> *vec)
  {
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    traits::row<Matrix>::type row(mat->getBaseMatrix());
    traits::col<Matrix>::type col(mat->getBaseMatrix());
    traits::const_value<Matrix>::type value(mat->getBaseMatrix());

    typedef traits::range_generator<major, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

    for (cursor_type cursor = begin<major>(mat->getBaseMatrix()), cend = end<major>(mat->getBaseMatrix()); cursor != cend; ++cursor)
      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
	if (value(*icursor) != 0.0) {
	  int r = mapLocalGlobalDOFs[row(*icursor)];
	  int c = mapLocalGlobalDOFs[col(*icursor)];
	  double v = value(*icursor);

	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
	}


    MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);

    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
    }
  }


  void ParallelDomainProblemBase::solvePetscMatrix(DOFVector<double> *vec)
  {
    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
    KSPGetPC(ksp, &pc);
    PCSetType(pc, PCJACOBI);
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
    KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, 0);
    KSPSolve(ksp, petscRhsVec, petscSolVec);

    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    int counter = 0;
    for (dofIt.reset(); !dofIt.end(); ++dofIt)
      *dofIt = vecPointer[counter++];

    VecRestoreArray(petscSolVec, &vecPointer);

    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
    
    int i = 0;
    for (std::map<int, std::vector<DegreeOfFreedom> >::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end();
	 ++sendIt, i++) {
      sendBuffers[i] = new double[sendIt->second.size()];

      for (int j = 0; j < static_cast<int>(sendIt->second.size()); j++)
	sendBuffers[i][j] = (*vec)[(sendIt->second)[j]];

      mpiComm.Isend(sendBuffers[i], sendIt->second.size(), MPI_DOUBLE, sendIt->first, 0);
    }

    i = 0;
    for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
      recvBuffers[i] = new double[recvIt->second.size()];      

      mpiComm.Irecv(recvBuffers[i], recvIt->second.size(), MPI_DOUBLE, recvIt->first, 0);
    }

    
    mpiComm.Barrier();
    
    i = 0;
    for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
	(*vec)[(recvIt->second)[j]] = recvBuffers[i][j];

      delete [] recvBuffers[i];
    }
    
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];    
  }


  double ParallelDomainProblemBase::setElemWeights(AdaptInfo *adaptInfo) 
  {
    double localWeightSum = 0.0;
    int elNum = -1;

    elemWeights.clear();

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					 Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      Element *element = elInfo->getElement();

      // get partition data
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));

      if (partitionData && partitionData->getPartitionStatus() == IN) {
	if (partitionData->getLevel() == 0) {
	  elNum = element->getIndex();
	}
	TEST_EXIT(elNum != -1)("invalid element number\n");
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }


  void ParallelDomainProblemBase::partitionMesh(AdaptInfo *adaptInfo)
  {
    if (initialPartitionMesh) {
      initialPartitionMesh = false;
      partitioner->fillCoarsePartitionVec(&oldPartitionVec);
      partitioner->partition(&elemWeights, INITIAL);
    } else {
      oldPartitionVec = partitionVec;
      partitioner->partition(&elemWeights, ADAPTIVE_REPART, 100.0 /*0.000001*/);
    }    

    partitioner->fillCoarsePartitionVec(&partitionVec);
  }


  void ParallelDomainProblemBase::createInteriorBoundaryInfo()
  {
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
    while (elInfo) {
      Element *element = elInfo->getElement();

      // Hidde elements which are not part of ranks partition.
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
      if (partitionData->getPartitionStatus() == IN) {
	for (int i = 0; i < 3; i++) {
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->getElementData(PARTITION_ED));
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
 	    AtomicBoundary& bound = interiorBoundary.
	      getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]);
 	    bound.rankObject.el = element;
 	    bound.rankObject.subObjAtBoundary = EDGE;
 	    bound.rankObject.ithObjAtBoundary = i;
 	    bound.neighbourObject.el = elInfo->getNeighbour(i);
 	    bound.neighbourObject.subObjAtBoundary = EDGE;
 	    bound.neighbourObject.ithObjAtBoundary = -1;
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
  }


  void ParallelDomainProblemBase::removeMacroElements()
  {
    std::vector<MacroElement*> macrosToRemove;
    for (std::deque<MacroElement*>::iterator it = mesh->firstMacroElement();
	 it != mesh->endOfMacroElements();
	 ++it) {
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>
	((*it)->getElement()->getElementData(PARTITION_ED));
      if (partitionData->getPartitionStatus() != IN)
	macrosToRemove.push_back(*it);
    }

    mesh->removeMacroElements(macrosToRemove);
  }


  void ParallelDomainProblemBase::createLocalGlobalNumbering(int& nRankDOFs, 
							     int& nOverallDOFs)
  {
    /// === Get rank information about DOFs. ===

    // Stores to each DOF pointer the set of ranks the DOF is part of.
    std::map<const DegreeOfFreedom*, std::set<int> > partitionDOFs;
    // Set of all DOFs of the rank.
    std::vector<const DegreeOfFreedom*> rankDOFs;
    // Set of all interior boundary DOFs in ranks partition which are owned by 
    // another rank.
    std::map<const DegreeOfFreedom*, int> boundaryDOFs;

    createDOFMemberInfo(partitionDOFs, rankDOFs, boundaryDOFs);

    nRankDOFs = rankDOFs.size();
    nOverallDOFs = partitionDOFs.size();

    // === Get starting position for global rank dof ordering ====

    int rstart = 0;
    MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
    rstart -= nRankDOFs;

   
    /// === Create information which dof indices must be send and which must be received. ===

    std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> > sendNewDofs;
    std::map<int, std::vector<DegreeOfFreedom> > recvNewDofs;

    for (std::map<const DegreeOfFreedom*, int>::iterator it = boundaryDOFs.begin();
	 it != boundaryDOFs.end();
	 ++it) {

      if (it->second == mpiRank) {
	// If the boundary dof is a rank dof, it must be send to other ranks.

	// old global index
	int oldDofIndex = (it->first)[0];
	// search for new dof index in ranks partition for this boundary dof
	int newDofIndex = 0;
	for (int i = 0; i < static_cast<int>(rankDOFs.size()); i++) {
	  if (rankDOFs[i] == it->first) {
	    newDofIndex = rstart + i;
	    break;
	  }
	}

	// Search for all ranks that have this dof too.
	for (std::set<int>::iterator itRanks = partitionDOFs[it->first].begin();
	     itRanks != partitionDOFs[it->first].end();
	     ++itRanks) {
	  if (*itRanks != mpiRank)
	    sendNewDofs[*itRanks][oldDofIndex] = newDofIndex;
	}
      } else {
	// If the boundary dof is not a rank dof, its new dof index, and later
	// also the dof values, must be received from another rank.
	recvNewDofs[it->second].push_back((it->first)[0]);
      }
    }


    /// === Send and receive the dof indices at boundary. ===

    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());
    
    int i = 0;
    for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin();
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
      sendBuffers[i] = new int[sendIt->second.size() * 2];
      int c = 0;
      for (std::map<DegreeOfFreedom, DegreeOfFreedom>::iterator dofIt = sendIt->second.begin();
	   dofIt != sendIt->second.end();
	   ++dofIt, c += 2) {
	sendBuffers[i][c] = dofIt->first;
	sendBuffers[i][c + 1] = dofIt->second;

	sendDofs[sendIt->first].push_back(dofIt->second);
      }

      mpiComm.Isend(sendBuffers[i], sendIt->second.size() * 2, MPI_INT, sendIt->first, 0);
    }

    i = 0;
    for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
      recvBuffers[i] = new int[recvIt->second.size() * 2];

      mpiComm.Irecv(recvBuffers[i], recvIt->second.size() * 2, MPI_INT, recvIt->first, 0);
    }


    mpiComm.Barrier();

    
    /// === Delete send buffers. ===

    i = 0;
    for (std::map<int, std::map<DegreeOfFreedom, DegreeOfFreedom> >::iterator sendIt = sendNewDofs.begin();
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) 
      delete [] sendBuffers[i];


    /// === Change dof indices for rank partition. ===

    for (int i = 0; i < static_cast<int>(rankDOFs.size()); i++) {
      const_cast<DegreeOfFreedom*>(rankDOFs[i])[0] = i; 
      mapLocalGlobalDOFs[i] = rstart + i;
      mapGlobalLocalDOFs[rstart + i] = i;
      isRankDOF[i] = true;
    }


    /// === Change dof indices at boundary from other ranks. ===

    i = 0;
    for (std::map<int, std::vector<DegreeOfFreedom> >::iterator recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++) {

	int oldDof = recvBuffers[i][j * 2];
	int newDof = recvBuffers[i][j * 2 + 1];
	int newLocalDof = mapLocalGlobalDOFs.size();

	recvDofs[recvIt->first].push_back(newDof);

	for (std::map<const DegreeOfFreedom*, int>::iterator dofIt = boundaryDOFs.begin();
	     dofIt != boundaryDOFs.end();
	     ++dofIt) {
	  if ((dofIt->first)[0] == oldDof) {
	    const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newLocalDof;
	    mapLocalGlobalDOFs[newLocalDof] = newDof;
	    mapGlobalLocalDOFs[newDof] = newLocalDof;
	    isRankDOF[newLocalDof] = false;
	    break;
	  }
	}
      }

      delete [] recvBuffers[i];
    }


    /// === Create local information from sendDofs and recvDofs

    for (std::map<int, std::vector<DegreeOfFreedom> >::iterator it = sendDofs.begin();
	 it != sendDofs.end();
	 ++it)
      for (std::vector<DegreeOfFreedom>::iterator dofIt = it->second.begin();
	   dofIt != it->second.end();
	   dofIt++)
	*dofIt = mapGlobalLocalDOFs[*dofIt];

    for (std::map<int, std::vector<DegreeOfFreedom> >::iterator it = recvDofs.begin();
	 it != recvDofs.end();
	 ++it)
      for (std::vector<DegreeOfFreedom>::iterator dofIt = it->second.begin();
	   dofIt != it->second.end();
	   dofIt++)
	*dofIt = mapGlobalLocalDOFs[*dofIt];  
  }



  void ParallelDomainProblemBase::createDOFMemberInfo(
		       std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
		       std::vector<const DegreeOfFreedom*>& rankDOFs,
		       std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
  {
    /// === Determine to each dof the set of partitions the dof belongs to. ===

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *element = elInfo->getElement();

      // Determine to each dof the partition(s) it corresponds to.
      for (int i = 0; i < 3; i++) 
	partitionDOFs[element->getDOF(i)].insert(partitionVec[element->getIndex()]);
          
      elInfo = stack.traverseNext(elInfo);
    }

    /// === Determine the set of ranks dofs and the dofs ownership at the boundary. ===

    for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDOFs.begin();
	 it != partitionDOFs.end();
	 ++it) {
      for (std::set<int>::iterator itpart1 = it->second.begin();
	   itpart1 != it->second.end();
	   ++itpart1) {
	if (*itpart1 == mpiRank) {
	  if (it->second.size() == 1) {
	    rankDOFs.push_back(it->first);
	  } else {	    
	    // This dof is at the ranks boundary. It is owned by the rank only if
	    // the rank number is the highest of all ranks containing this dof.

	    bool insert = true;
	    int highestRank = mpiRank;
	    for (std::set<int>::iterator itpart2 = it->second.begin();
		 itpart2 != it->second.end();
		 ++itpart2) {
	      if (*itpart2 > mpiRank)
		insert = false;

	      if (*itpart2 > highestRank)
		highestRank = *itpart2;
	    }

	    if (insert)
	      rankDOFs.push_back(it->first);

	    boundaryDOFs[it->first] = highestRank;
	  }

	  break;
	}
      }
    }
  }


  ParallelDomainProblemScal::ParallelDomainProblemScal(const std::string& name,
						       ProblemScal *problem,
						       ProblemInstatScal *problemInstat)
    : ParallelDomainProblemBase(name, problem, problemInstat, problem->getFESpace()),
      probScal(problem)
  {
  }

  void ParallelDomainProblemScal::initParallelization(AdaptInfo *adaptInfo)
  {
    ParallelDomainProblemBase::initParallelization(adaptInfo);

    probScal->getSystemMatrix()->setIsRankDOF(isRankDOF);
  }

  Flag ParallelDomainProblemScal::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    //    return iterationIF->oneIteration(adaptInfo, toDo);

    Flag flag = dynamic_cast<StandardProblemIteration*>(iterationIF)->buildAndAdapt(adaptInfo, toDo);

    fillPetscMatrix(probScal->getSystemMatrix(), probScal->getRHS());

    solvePetscMatrix(probScal->getSolution());

//     if (toDo.isSet(SOLVE))
//       iterationIF->getProblem()->solve(adaptInfo, false);

//     if (toDo.isSet(SOLVE_RHS))
//       iterationIF->getProblem()->solve(adaptInfo, true);

//     if (toDo.isSet(ESTIMATE)) 
//       iterationIF->getProblem()->estimate(adaptInfo);

    return flag;

  }

}