#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"

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


    /// === Determine to each dof the set of partitions the dof belongs to. ===

    std::map<const DegreeOfFreedom*, std::set<int> > partitionDofs;
    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. ===

    std::vector<const DegreeOfFreedom*> rankDofs;
    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;
	  }
	}
      }
    }


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

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


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

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

    // === Create local and global dofs ordering. ===

    int *gOrder = (int*)(malloc(sizeof(int) * rankDofs.size()));
    int *lOrder = (int*)(malloc(sizeof(int) * rankDofs.size()));

    for (std::vector<const DegreeOfFreedom*>::iterator it = rankDofs.begin();
	 it != rankDofs.end(); ++it) {
      gOrder[nRankDOFs++] = (*it)[0];
    }

    int rstart = 0;
    MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
    rstart -= nRankDOFs;
   
    for (int i = 0; i < nRankDOFs; i++) {
      lOrder[i] = rstart + i;
    }

    AOCreateBasic(PETSC_COMM_WORLD, nRankDOFs, gOrder, lOrder, &applicationOrdering);
    
    free(gOrder);
    free(lOrder);

    /// === 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) {
	int oldDofIndex = (it->first)[0];
	int newDofIndex = 0;
	for (int i = 0; i < static_cast<int>(rankDofs.size()); i++) {
	  if (rankDofs[i] == it->first) {
	    newDofIndex = rstart + i;
	    break;
	  }
	}

	for (std::set<int>::iterator itRanks = partitionDofs[it->first].begin();
	     itRanks != partitionDofs[it->first].end();
	     ++itRanks) {
	  if (*itRanks != mpiRank) {
	    sendNewDofs[*itRanks][oldDofIndex] = newDofIndex;
	  }
	}
      } else {
	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;
      }

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

    /// === 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++) {
	for (std::map<const DegreeOfFreedom*, int>::iterator dofIt = boundaryDofs.begin();
	     dofIt != boundaryDofs.end();
	     ++dofIt) {
	  if ((dofIt->first)[0] == recvBuffers[i][j * 2]) {
	    const_cast<DegreeOfFreedom*>(dofIt->first)[0] = recvBuffers[i][j * 2 + 1];
	    break;
	  }
	}
      }

      delete [] recvBuffers[i];
    }

    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] = rstart + i;
    }

    /// === Create petsc matrix. ===
    int ierr;
    ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
    ierr = MatSetSizes(petscMatrix, rankDofs.size(), rankDofs.size(),
		       partitionDofs.size(), partitionDofs.size());
    ierr = MatSetType(petscMatrix, MATAIJ);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
    ierr = VecSetSizes(petscRhsVec, rankDofs.size(), partitionDofs.size());    
    ierr = VecSetType(petscRhsVec, VECMPI);
  }

  void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
  {
    AODestroy(applicationOrdering);
  }

  void ParallelDomainProblemBase::fillPetscMatrix(DOFMatrix *mat, 
						  DOFVector<double> *vec)
  {
    DOFMatrix::Iterator rowIt(mat, USED_DOFS);
    for (rowIt.reset(); !rowIt.end(); ++rowIt) {
      for (int i = 0; i < static_cast<int>((*rowIt).size()); i++) {
	if ((*rowIt)[i].col >= 0) {
	  MatSetValues(petscMatrix, 1, &i, 1, &((*rowIt)[i].col), &((*rowIt)[i].entry), 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 = dofIt.getDOFIndex();
      double value = *dofIt;

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

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

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

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

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

  }

}