#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,
						       RefinementManager *refineManager)
    : iterationIF(iIF),
      timeIF(tIF),
      feSpace(fe),
      mesh(fe->getMesh()),
      refinementManager(refineManager),
      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. ===

    // 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;
    // Number of DOFs in ranks partition that are owned by the rank.
    int nRankDOFs = 0;
    // Number of DOFs in ranks partition that are at an interior boundary and are
    // owned by other ranks.
    int nOverallDOFs = 0;

    createLocalGlobalNumbering(rankDOFs, boundaryDOFs, nRankDOFs, nOverallDOFs);


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

    createInteriorBoundaryInfo(rankDOFs, boundaryDOFs);


    // === 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++) {
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));

      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
      for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
 	admin.setDOFFree(j, false);

      admin.setUsedSize(mapLocalGlobalDOFs.size() - 1);
      admin.setUsedCount(mapLocalGlobalDOFs.size());
      admin.setFirstHole(mapLocalGlobalDOFs.size());
    }


    /// === Global refinements. ===

    refinementManager->globalRefine(mesh, 1);

    updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);

    exit(0);

    /// === 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(std::vector<const DegreeOfFreedom*>& rankDOFs,
							     std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
  {
    FUNCNAME("ParallelDomainProblemBase::createInteriorBoundaryInfo()");

    // === First, create all the information about the interior boundaries. ===

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

      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) {
	    // We have found an element that is at an interior boundary. 

	    // === Find out, if the boundary part of the element corresponds to the
	    //     rank or to the rank "on the other side" of the interoir boundary. ===

	    const DegreeOfFreedom* boundDOF1 = NULL;
	    const DegreeOfFreedom* boundDOF2 = NULL;
	    
	    switch (i) {
	    case 0:
	      boundDOF1 = element->getDOF(1);
	      boundDOF2 = element->getDOF(2);
	      break;
	    case 1:
	      boundDOF1 = element->getDOF(0);
	      boundDOF2 = element->getDOF(2);
	      break;
	    case 2:
	      boundDOF1 = element->getDOF(0);
	      boundDOF2 = element->getDOF(1);
	      break;
	    default:
	      ERROR_EXIT("Should never happen!\n");
	    }

	    bool isRankDOF1 = (find(rankDOFs.begin(), rankDOFs.end(), boundDOF1) != rankDOFs.end());
	    bool isRankDOF2 = (find(rankDOFs.begin(), rankDOFs.end(), boundDOF2) != rankDOFs.end());
	    bool ranksBoundary = isRankDOF1 || isRankDOF2;

	    /// === And add the part of the interior boundary. ===

	    AtomicBoundary& bound = 
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.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;
	    std::cout << "ADD IN " << mpiRank << ": " << element->getIndex() << " " << elInfo->getNeighbour(i)->getIndex() << std::endl;
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }


    // === Once we have this information, we must care about their order. Eventually ===
    // === all the boundaries have to be in the same order on both ranks (because    ===
    // === each boundary is shared by exactly two ranks).                            ===

    std::vector<int*> sendBuffers;
    std::vector<int*> recvBuffers;

    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end();
	 ++rankIt) {
      int* buffer = new int[rankIt->second.size()];
      for (int i = 0; i < static_cast<int>(rankIt->second.size()); i++)
	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
      sendBuffers.push_back(buffer);
      
      mpiComm.Isend(buffer, rankIt->second.size(), MPI_INT, rankIt->first, 0);
    }

    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end();
	 ++rankIt) {
      int *buffer = new int[rankIt->second.size()];
      recvBuffers.push_back(buffer);
      
      mpiComm.Irecv(buffer, rankIt->second.size(), MPI_INT, rankIt->first, 0);      
    }

    mpiComm.Barrier();

    int i = 0;
    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end();
	 ++rankIt) {

      // === We have received from rank "rankIt->first" the ordered list of element ===
      // === indices. We now have to sort the corresponding list in this rank to    ===
      // === get the same order.                                                    ===
      
      for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
	// If the expected object is not at place, search for it.
	if ((rankIt->second)[j].neighbourObject.el->getIndex() != recvBuffers[i][j]) {
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
	    if ((rankIt->second)[k].neighbourObject.el->getIndex() == recvBuffers[i][j])
	      break;

	  // The element must always be found, because the list is just in another order.
	  TEST_EXIT(k < rankIt->second.size())("Should never happen!\n");

	  // Swap the current with the found element.
	  AtomicBoundary tmpBound = (rankIt->second)[k];
	  (rankIt->second)[k] = (rankIt->second)[j];
	  (rankIt->second)[j] = tmpBound;	
	}
      }

      delete [] recvBuffers[i++];
    }

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


  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(std::vector<const DegreeOfFreedom*>& rankDOFs,
							     std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
							     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;

    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::updateLocalGlobalNumbering(int& nRankDOFs, 
							     int& nOverallDOFs)
  {
    FUNCNAME("ParallelDomainProblemBase::updateLocalGlobalNumbering()");

    std::set<const DegreeOfFreedom*> rankDOFs;
    std::map<const DegreeOfFreedom*, int> boundaryDOFs;

    
    /// === Get all DOFs in ranks partition. ===

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *element = elInfo->getElement();
      
      for (int i = 0; i < 3; i++) 
	rankDOFs.insert(element->getDOF(i));

      elInfo = stack.traverseNext(elInfo);
    }


    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
    // === rankDOFs to boundaryDOFs                                                ===

    for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
	   myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end();
	 ++it) {
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end();
	   ++boundIt) {
	const DegreeOfFreedom *dof1 = NULL;
	const DegreeOfFreedom *dof2 = NULL;

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
	  dof1 = boundIt->rankObject.el->getDOF(1);
	  dof2 = boundIt->rankObject.el->getDOF(2);
	  break;
	case 1:
	  dof1 = boundIt->rankObject.el->getDOF(0);
	  dof2 = boundIt->rankObject.el->getDOF(2);
	  break;
	case 2:
	  dof1 = boundIt->rankObject.el->getDOF(0);
	  dof2 = boundIt->rankObject.el->getDOF(1);
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

	std::vector<const DegreeOfFreedom*> boundDOFs;
	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		   boundDOFs);
      }
    }    
  }


  void ParallelDomainProblemBase::addAllDOFs(Element *el, int ithEdge, 
					     std::vector<const DegreeOfFreedom*>& dofs,
					     bool addVertices)
  {
    FUNCNAME("ParallelDomainProblemBase::addAllDOFs()");

    const DegreeOfFreedom* boundDOF1 = NULL;
    const DegreeOfFreedom* boundDOF2 = NULL;

    if (addVertices) {
      switch (ithEdge) {
      case 0:
	boundDOF1 = el->getDOF(1);
	boundDOF2 = el->getDOF(2);
	break;
      case 1:
	boundDOF1 = el->getDOF(0);
	boundDOF2 = el->getDOF(2);
	break;
      case 2:
	boundDOF1 = el->getDOF(0);
	boundDOF2 = el->getDOF(1);
	break;
      default:
	ERROR_EXIT("Should never happen!\n");
      }

      dofs.push_back(boundDOF1);
    }

    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
	addAllDOFs(el->getSecondChild()->getFirstChild(), 0, dofs, false);
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
	addAllDOFs(el->getSecondChild()->getSecondChild(), 1, dofs, false);
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
	addAllDOFs(el->getFirstChild()->getFirstChild(), 0, dofs, false);
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
	addAllDOFs(el->getFirstChild()->getSecondChild(), 1, dofs, false);
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
	addAllDOFs(el->getFirstChild(), 0, dofs, false);
	dofs.push_back(el->getFirstChild()->getDOF(2));
	addAllDOFs(el->getSecondChild(), 1, dofs, false);
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }

    if (addVertices)
      dofs.push_back(boundDOF2);		   
  }


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

}