ParallelDomainProblem.cc

#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 "ElementDofIterator.h"

#include "petscksp.h"

namespace AMDiS {

  PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
  {
    if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
      std::cout << "  Iteration " << iter << ": " << rnorm << "\n";

    return 0;
  }

  ParallelDomainBase::ParallelDomainBase(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 ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
  {
    if (mpiSize <= 1)
      return;

    // Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
    // only for macro meshes, so it will not work yet if the mesh is already refined
    // in some way.
    testForMacroMesh();

    // 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;
    // Number of DOFs in ranks partition that are owned by the rank.
    int nRankDOFs = 0;
    // Number of all DOFs in the macro mesh.
    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());
      admin.setUsedCount(mapLocalGlobalDOFs.size());
      admin.setFirstHole(mapLocalGlobalDOFs.size());
    }

    exit(0);

    // === Global refinements. ===

    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);

      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
    }


// #if (DEBUG != 0)
//     testInteriorBoundary();
// #endif

    // === 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 ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
  {
  }


  void ParallelDomainBase::testForMacroMesh()
  {
    FUNCNAME("ParallelDomainBase::testForMacroMesh()");

    int nMacroElements = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      TEST_EXIT(elInfo->getLevel() == 0)
	("Mesh is already refined! This does nor work with parallelization!\n");
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

    TEST_EXIT(nMacroElements >= mpiSize)
      ("The mesh has less macro elements than number of mpi processes!\n");
  }


  void ParallelDomainBase::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 ParallelDomainBase::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, myKSPMonitor, 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());

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

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

      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
    }

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

      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
    }


    MPI::Request::Waitall(requestCounter, request);

    
    i = 0;
    for (std::map<int, std::vector<const 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][0]] = recvBuffers[i][j];

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


  double ParallelDomainBase::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_DBG(elNum != -1)("invalid element number\n");
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }


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

      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;

    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

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

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


    MPI::Request::Waitall(requestCounter, request);


    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_DBG(k < static_cast<int>(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 ParallelDomainBase::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 ParallelDomainBase::createLocalGlobalNumbering(
                                 std::vector<const DegreeOfFreedom*>& rankDOFs,
				 std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
				 int& nRankDOFs, 
				 int& nOverallDOFs)
  {
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");

    // === 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;
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
    rstart -= nRankDOFs;

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

    // Stores to each rank a map from DOF pointers (which are all owned by the rank
    // and lie on an interior boundary) to new global DOF indices.
    std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> > sendNewDofs;

    // Maps to each rank the number of new DOF indices this rank will receive from.
    // All these DOFs are at an interior boundary on this rank, but are owned by
    // another rank.
    std::map<int, int> 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.

	// search for new dof index in ranks partition for this boundary dof
	DegreeOfFreedom newDofIndex = 0;
	for (int i = 0; i < nRankDOFs; 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][it->first] = 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.
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
      }
    }


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

    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

    sendDofs.clear();
    recvDofs.clear();
    
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

    int i = 0;
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];

      int c = 0;
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
	   dofIt != sendIt->second.end();
	   ++dofIt) {
	sendBuffers[i][c++] = (dofIt->first)[0];
	sendBuffers[i][c++] = dofIt->second;

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

      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
    }

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

      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
    }


    MPI::Request::Waitall(requestCounter, request);

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

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


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

    mapLocalGlobalDOFs.clear();
    mapGlobalLocalDOFs.clear();
    isRankDOF.clear();

    for (int i = 0; i < nRankDOFs; 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, int>::iterator recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

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

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

	bool found = false;

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

	  if ((dofIt->first)[0] == oldDof) {
	    recvDofs[recvIt->first].push_back(dofIt->first);
	    const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newLocalDof;
	    mapLocalGlobalDOFs[newLocalDof] = newGlobalDof;
	    mapGlobalLocalDOFs[newGlobalDof] = newLocalDof;
	    isRankDOF[newLocalDof] = false;
	    found = true;
	    break;
	  }
	}

	TEST_EXIT_DBG(found)("Should not happen!\n");
      }

      delete [] recvBuffers[i];
    }
  }


  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, 
						      int& nOverallDOFs)
  {
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");

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

    std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDofs = sendDofs;
    std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDofs = recvDofs;

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

	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
	  ("Should never happen!\n");
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
	  ("Should never happen!\n");

	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];

	std::vector<const DegreeOfFreedom*> boundDOFs;
	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		   boundDOFs);
	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
	  sendNewDofs[it->first].push_back(boundDOFs[i]);
	}
      }
    }    


    for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
	   otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.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");
	}

	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
	  ("Should never happen!\n");
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
	  ("Should never happen!\n");

	rankDOFs.erase(dof1);
	rankDOFs.erase(dof2);
	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];

	std::vector<const DegreeOfFreedom*> boundDOFs;
	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		   boundDOFs);
	for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) {
	//	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  rankDOFs.erase(boundDOFs[i]);
	  newBoundaryDOFs[boundDOFs[i]] = it->first;
	  recvNewDofs[it->first].push_back(boundDOFs[i]);
	}
      }

    }

    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
    rstart -= nRankDOFs;


    // === Calculate number of overall DOFs of all partitions. ===

    mpiComm.Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);


    // === Create new local DOF index numbering. ===

    mapLocalGlobalDOFs.clear();
    mapGlobalLocalDOFs.clear();
    isRankDOF.clear();

    int i = 0;
    for (std::set<const DegreeOfFreedom*>::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end();
	 ++dofIt, i++) {
      const_cast<DegreeOfFreedom*>(*dofIt)[0] = i;
      mapLocalGlobalDOFs[i] = rstart + i;
      mapGlobalLocalDOFs[rstart + i] = i;
      isRankDOF[i] = true;
    }


    // === Send new DOF indices. ===

    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

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

      request[requestCounter++] = 
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
    }

    i = 0;
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
      int nRecvDofs = recvIt->second.size();
      recvBuffers[i] = new int[nRecvDofs];
	
      request[requestCounter++] = 
	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
    }


    MPI::Request::Waitall(requestCounter, request);


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

    i = 0;
    int newDofIndex = nRankDOFs;
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end(); ++recvIt) {
      
      int j = 0;
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = recvIt->second.begin();
	   dofIt != recvIt->second.end();
	   ++dofIt) {
	const_cast<DegreeOfFreedom*>(*dofIt)[0] = newDofIndex;
	mapLocalGlobalDOFs[newDofIndex] = recvBuffers[i][j];
	mapGlobalLocalDOFs[recvBuffers[i][j]] = newDofIndex;
	isRankDOF[newDofIndex] = false;
	newDofIndex++;
	j++;
      }

      delete [] recvBuffers[i++];
    }

    // === Update list of dofs that must be communicated for solution exchange. ===

    sendDofs = sendNewDofs;
    recvDofs = recvNewDofs;
  }


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

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


  void ParallelDomainBase::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. ===

    ElementDofIterator elDofIt(feSpace->getAdmin());

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

      elDofIt.reset(element);
      do {
	// Determine to each dof the partition(s) it corresponds to.
	partitionDOFs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
      } while(elDofIt.next());

      elInfo = stack.traverseNext(elInfo);
    }


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

    // iterate over all DOFs
    for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDOFs.begin();
	 it != partitionDOFs.end();
	 ++it) {

      // iterate over all partition the current DOF is part of.
      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;
	}

      }
    }
  }


  void ParallelDomainBase::testInteriorBoundary()
  {
    FUNCNAME("ParallelDomainBase::testInteriorBoundary()");

    // Maps to each neighbour rank an array of WorldVectors. This array contain the 
    // coordinates of all dofs this rank shares on the interior boundary with the 
    // neighbour rank. A rank sends the coordinates to another rank, if it owns the
    // boundarys dof.
    std::map<int, std::vector<WorldVector<double> > > sendCoords;
    // A rank receives all boundary dofs that are at its interior boundaries but are
    // not owned by the rank. This map stores for each rank the the coordinates of dofs
    // this rank expectes to receive from.
    std::map<int, std::vector<WorldVector<double> > > recvCoords;

    WorldVector<double> coords;
   
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator it = 
	   sendDofs.begin();
	 it != sendDofs.end(); ++it) {
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
	TEST_EXIT(b)("Cannot find DOF in mesh!\n");
	sendCoords[it->first].push_back(coords);
      }
    }

    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator it = 
	   recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
	TEST_EXIT(b)("Cannot find DOF in mesh!\n");
	recvCoords[it->first].push_back(coords);
      }
    }

    for (std::map<int, std::vector<WorldVector<double> > >::iterator it 
	   = sendCoords.begin();
	 it != sendCoords.end(); ++it) {
    }

    for (std::map<int, std::vector<WorldVector<double> > >::iterator it = recvCoords.begin();
	 it != recvCoords.end(); ++it) {
    }
  }


  ParallelDomainScal::ParallelDomainScal(const std::string& name,
					 ProblemScal *problem,
					 ProblemInstatScal *problemInstat)
    : ParallelDomainBase(name, 
			 problem, 
			 problemInstat, 
			 problem->getFESpace(),
			 problem->getRefinementManager()),
      probScal(problem)
  {
    info = problem->getInfo();
  }

  void ParallelDomainScal::initParallelization(AdaptInfo *adaptInfo)
  {
    FUNCNAME("ParallelDomainScal::initParallelization()");

    ParallelDomainBase::initParallelization(adaptInfo);

    DOFMatrix* m = probScal->getSystemMatrix();

    TEST_EXIT(m)("No DOF Matrix!\n");

    m->setIsRankDOF(isRankDOF);
  }

  void ParallelDomainScal::solve()
  {
    FUNCNAME("ParallelDomainScal::solve()");

#ifdef _OPENMP
    double wtime = omp_get_wtime();
#endif
    clock_t first = clock();

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

#ifdef _OPENMP
    INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
		   TIME_USED(first, clock()),
		   omp_get_wtime() - wtime);
#else
    INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
		   TIME_USED(first, clock()));
#endif
  }

  Flag ParallelDomainScal::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    FUNCNAME("ParallelDomainScal::oneIteration()");

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

    if (toDo.isSet(SOLVE))
      solve();

    if (toDo.isSet(SOLVE_RHS))
      ERROR_EXIT("Not yet implemented!\n");

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

    return flag;
  }

}