ParallelProblem.cc

#include <queue>
#include <time.h>
#include "boost/lexical_cast.hpp"
#include "ParallelProblem.h"
#include "ProblemScal.h"
#include "ProblemVec.h"
#include "ProblemInstat.h"
#include "AdaptInfo.h"
#include "AdaptStationary.h"
#include "ConditionalEstimator.h"
#include "ConditionalMarker.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "DualTraverse.h"
#include "MeshStructure.h"
#include "DOFVector.h"
#include "FiniteElemSpace.h"
#include "RefinementManager.h"
#include "CoarseningManager.h"
#include "Lagrange.h"
#include "ElementFileWriter.h"
#include "MacroWriter.h"
#include "ValueWriter.h"
#include "SystemVector.h"
#include "VtkWriter.h"
#include "mpi.h"

namespace AMDiS {

  using boost::lexical_cast;

  bool elementInPartition(ElInfo *elInfo)
  {
    PartitionElementData *elementData = dynamic_cast<PartitionElementData*>
      (elInfo->getElement()->getElementData(PARTITION_ED));
    if (elementData && elementData->getPartitionStatus() == IN) {
      return true;
    } else {
      return false;
    }
  }

  class MyDualTraverse : public DualTraverse
  {
  public:
    MyDualTraverse(int coarseLevel)
      : coarseLevel_(coarseLevel)
    {}

    bool skipEl1(ElInfo *elInfo)
    {
      PartitionElementData *elementData = dynamic_cast<PartitionElementData*>
	(elInfo->getElement()->getElementData(PARTITION_ED));
      if (elementData) {
	if (elInfo->getElement()->isLeaf() && 
	   elementData->getLevel() < coarseLevel_)
	  return false;
	if(elementData->getLevel() == coarseLevel_)
	  return false;
      }
      return true;
    }
    
    /*bool skipEl2(ElInfo *elInfo)
    {
	    return skipEl1(elInfo);
    }*/
  private:
    int coarseLevel_;
  };

  // =========================================================================
  // ===== class ParallelProblemBase =========================================
  // =========================================================================

  ParallelProblemBase::ParallelProblemBase(std::string name,
					   ProblemIterationInterface *iterationIF,
					   ProblemTimeInterface *timeIF)
    : iterationIF_(iterationIF),
      timeIF_(timeIF),
      debugMode(0),
      debugServerProcess(false)
  {
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();

    GET_PARAMETER(0, name + "->debug mode", "%d", &debugMode);

    groupIsSet=false;
    if (debugMode) {
      MPI::Group group = MPI::COMM_WORLD.Get_group();
      
      int rankSize = mpiSize - 1;
      int *ranks = new int[rankSize];
      for (int i = 0; i < rankSize; i++)
	ranks[i] = i + 1;

      std::cout<<"set amdis-group"<<std::endl;
      amdisGroup = group.Incl(rankSize, ranks);
      groupIsSet=true;
      mpiComm = MPI::COMM_WORLD.Create(amdisGroup);
      
      if (mpiComm != MPI::COMM_NULL) {
	mpiRank = mpiComm.Get_rank();
	mpiSize = mpiComm.Get_size();
	debugServerProcess = false;
      } else {
	debugServerProcess = true;
      }
      
      delete [] ranks;
    } else {
      mpiComm = MPI::COMM_WORLD;
    }
  }

  void ParallelProblemBase::exitParallelization(AdaptInfo *adaptInfo)
  {
    if (!timeIF_) 
      closeTimestep(adaptInfo);
#ifdef DEBUG
    std::cout<<"free amdisGroup"<<std::endl;
#endif
    if(groupIsSet)
    	amdisGroup.Free();
#ifdef DEBUG
    std::cout<<"freeD amdisGroup"<<std::endl;
#endif    
  }

  void ParallelProblemBase::closeTimestep(AdaptInfo *adaptInfo)
  {
    if (mpiSize > 1 && doBuildGlobalSolution(adaptInfo)) {
      if (debugMode && mpiRank == 0) {
	// Send adaptInfo inforamtion to debug server
	double sendTime = adaptInfo->getTime();
	double sendTimestep = adaptInfo->getTimestep();
	MPI::COMM_WORLD.Send(&sendTime, 1, MPI_DOUBLE, 0, 100);
	MPI::COMM_WORLD.Send(&sendTimestep, 1, MPI_DOUBLE, 0, 100);
      }
      synchronizeMeshes(adaptInfo);	
      exchangeRankSolutions(adaptInfo);
      buildGlobalSolution(adaptInfo);
    }
    
    if (timeIF_) 
      timeIF_->closeTimestep(adaptInfo);
  }

  Flag ParallelProblemBase::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    Flag flag;
    
    if (mpiSize > 1 && toDo.isSet(MARK | ADAPT)) {
      flag = iterationIF_->oneIteration(adaptInfo, MARK | ADAPT);
      
      double localWeightSum = setElemWeights(adaptInfo);
      if (doPartitioning(adaptInfo, localWeightSum)) {
	clock_t partitioningStart = clock();

	synchronizeMeshes(adaptInfo);
	partitionMesh(adaptInfo);
	refineOverlap(adaptInfo);
	createOverlap(adaptInfo);
	synchronizeMeshes(adaptInfo);
	exchangeDOFVectors(adaptInfo);
	coarsenOutOfPartition(adaptInfo);
	
	clock_t partitioningEnd = clock();
	partitioningTime = TIME_USED(partitioningStart, 
				     partitioningEnd);
	computationStart = partitioningEnd;
      }
      
      flag |= iterationIF_->oneIteration(adaptInfo, toDo & ~(MARK | ADAPT));
    } else {
      flag = iterationIF_->oneIteration(adaptInfo, toDo);
    }
    
    // synchronize adaption flag
    unsigned long *flagBuffer = new unsigned long[mpiSize];
    
    unsigned long localFlag = flag.getFlags();
    
    mpiComm.Allgather(&localFlag, 1, MPI_UNSIGNED_LONG,
		      flagBuffer, 1, MPI_UNSIGNED_LONG);
    for (int i = 0; i < mpiSize; i++)
      flag.setFlag(flagBuffer[i]);
    
    delete [] flagBuffer;
    
    return flag;
  }

  // =========================================================================
  // ===== class ParallelProblem =============================================
  // =========================================================================

  ParallelProblem::ParallelProblem(std::string name,
				   ProblemIterationInterface *iterationIF,
				   ProblemTimeInterface *timeIF,
				   std::vector<DOFVector<double>*> vectors,
				   Mesh *mesh_,
				   RefinementManager *rm,
				   CoarseningManager *cm)
    : ParallelProblemBase(name, iterationIF, timeIF),
      name_(name),
      mesh(mesh_),
      refinementManager(rm),
      coarseningManager(cm),
      repartitionSteps_(1),
      puEveryTimestep_(false),
      dofVectors_(vectors),
      upperPartThreshold_(1.5),
      lowerPartThreshold_(2.0/3.0),
      globalCoarseGridLevel_(0),
      localCoarseGridLevel_(0),
      globalRefinements_(0),
      adaptiveThresholds_(0),
      thresholdIncFactor_(2.0),
      thresholdDecFactor_(0.5),
      repartTimeFactor_(10.0)
  {
    GET_PARAMETER(0, name_ + "->upper part threshold", "%f", 
		  &upperPartThreshold_);
    GET_PARAMETER(0, name_ + "->lower part threshold", "%f", 
		  &lowerPartThreshold_);
    GET_PARAMETER(0, name_ + "->global coarse grid level", "%d", 
		  &globalCoarseGridLevel_);
    GET_PARAMETER(0, name_ + "->local coarse grid level", "%d", 
		  &localCoarseGridLevel_);
    GET_PARAMETER(0, name_ + "->global refinements", "%d", 
		  &globalRefinements_);


    TEST_EXIT(localCoarseGridLevel_ >= globalCoarseGridLevel_)
      ("local coarse grid level < global coarse grid level\n");

    partitioner = new ParMetisPartitioner(mesh, &mpiComm);

    GET_PARAMETER(0, name_ + "->adaptive thresholds", "%d", 
		  &adaptiveThresholds_);
    GET_PARAMETER(0, name_ + "->threshold inc factor", "%f", 
		  &thresholdIncFactor_);
    GET_PARAMETER(0, name_ + "->threshold dec factor", "%f", 
		  &thresholdDecFactor_);
    GET_PARAMETER(0, name_ + "->repart time factor", "%f", 
		  &repartTimeFactor_);


    TEST_EXIT(lowerPartThreshold_ <= 1.0)("invalid lower part threshold\n");
    TEST_EXIT(upperPartThreshold_ >= 1.0)("invalid upper part threshold\n");

    if (adaptiveThresholds_) {
      TEST_EXIT(thresholdDecFactor_ <= 1.0)("invalid threshold dec factor\n");
      TEST_EXIT(thresholdIncFactor_ >= 1.0)("invalid threshold inc factor\n");
    }
    minUpperTH_ = upperPartThreshold_;
    maxLowerTH_ = lowerPartThreshold_;
  }

  ParallelProblem::~ParallelProblem() 
  {
    delete partitioner;
  }

  bool ParallelProblem::doPartitioning(AdaptInfo *adaptInfo, double localWeightSum) 
  {
    FUNCNAME("ParallelProblem::doPartitioning()");

    double *weightSum = new double[mpiSize];
    int *partArray = new int[mpiSize];
    int part = 0;

    mpiComm.Gather(&localWeightSum, 1, MPI_DOUBLE, weightSum, 1, MPI_DOUBLE, 0);

    if (mpiRank == 0) {

      double average = 0.0;
      for (int i = 0; i < mpiSize; i++)
	average += weightSum[i];

      average /= mpiSize;

      for (int i = 0; i < mpiSize; i++) {
	if ((weightSum[i] / average) > upperPartThreshold_) {
	  part = 1;
	  break;
	}
	if ((weightSum[i] / average) < lowerPartThreshold_) {
	  part = 1;
	  break;
	}
      }

      double computationTime = TIME_USED(computationStart, clock());
      if (adaptiveThresholds_) {

	bool timeOver = ((computationTime / partitioningTime) >= repartTimeFactor_);

	if (part == 1 && !timeOver) {
	  // inc thresholds
	  upperPartThreshold_ *= thresholdIncFactor_;
	  lowerPartThreshold_ /= thresholdIncFactor_;

	  // avoid repartitioning
	  part = 0;
	}
      
	if (part == 0 && timeOver) {
	  // dec thresholds
	  upperPartThreshold_ *= thresholdDecFactor_;
	  lowerPartThreshold_ /= thresholdDecFactor_;

	  upperPartThreshold_ = max(minUpperTH_, upperPartThreshold_);
	  lowerPartThreshold_ = min(maxLowerTH_, lowerPartThreshold_);
	}
      }

      for (int i = 0; i < mpiSize; i++) {
	partArray[i] = part;
      }      
    }

    mpiComm.Scatter(partArray, 1, MPI_INT,
		    &part, 1, MPI_INT, 0);
    
    delete [] weightSum;
    delete [] partArray;

    return (part == 1);
  }

  bool ParallelProblem::doBuildGlobalSolution(AdaptInfo *adaptInfo) 
  {
    return true;
  }

  void ParallelProblem::partitionMesh(AdaptInfo *adaptInfo)
  {
    static bool initial = true;

    if (initial) {
      initial = false;
      partitioner->fillCoarsePartitionVec(&oldPartitionVec);
      partitioner->partition(&elemWeights, INITIAL);
    } else {
      oldPartitionVec = partitionVec;
      partitioner->partition(&elemWeights, ADAPTIVE_REPART, 100.0 /*0.000001*/);
    }    

    partitioner->fillCoarsePartitionVec(&partitionVec);
  }

  void ParallelProblem::refineOverlap(AdaptInfo *adaptInfo)
  {
    int dim = mesh->getDim();
    bool finished = (localCoarseGridLevel_ == 0);

    while (!finished) {
      std::map<DegreeOfFreedom, int> inOut; // 1: in, 2: out, 3: border dof

      // mark in/out/border dofs
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      while (elInfo) {
	Element *element = elInfo->getElement();
	PartitionElementData *partitionData = 
	  dynamic_cast<PartitionElementData*>(elInfo->getElement()->
					      getElementData(PARTITION_ED));

	const DegreeOfFreedom **dofs = element->getDOF();

	if (partitionData->getPartitionStatus() == IN) {
	  for (int i = 0; i < dim + 1; i++) {
	    DegreeOfFreedom dof = dofs[i][0];
	    if (inOut[dof] == 2) 
	      inOut[dof] = 3;
	    if (inOut[dof] == 0) 
	      inOut[dof] = 1;
	  }
	} else {
	  for (int i = 0; i < dim + 1; i++) {
	    DegreeOfFreedom dof = dofs[i][0];
	    if (inOut[dof] == 1) 
	      inOut[dof] = 3;
	    if (inOut[dof] == 0) 
	      inOut[dof] = 2;
	  }
	}

	elInfo = stack.traverseNext(elInfo);
      }

      // refine overlap-border and inner elements
      finished = true;
      bool marked = false;
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      while (elInfo) {
	Element *element = elInfo->getElement();
	PartitionElementData *partitionData = 
	  dynamic_cast<PartitionElementData*>(elInfo->getElement()->getElementData(PARTITION_ED));

	int level = partitionData->getLevel();

	if (level < localCoarseGridLevel_) {
	  if (partitionData->getPartitionStatus() != IN) {
	    const DegreeOfFreedom **dofs = element->getDOF();
	    for (int i = 0; i < dim + 1; i++) {
	      DegreeOfFreedom dof = dofs[i][0];
	      if (inOut[dof] == 3) {
		element->setMark(1);
		marked = true;
		if ((level + 1) < localCoarseGridLevel_) 
		  finished = false;
		break;
	      }
	    }
	  } else {
	    element->setMark(1);
	    marked = true;
	    if ((level + 1) < localCoarseGridLevel_) 
	      finished = false;
	  }
	}

	elInfo = stack.traverseNext(elInfo);
      }
      if (marked) 
	refinementManager->refineMesh(mesh);
    }
  }

  void ParallelProblem::globalRefineOutOfPartition(AdaptInfo *adaptInfo)
  {
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(elInfo->getElement()->
					    getElementData(PARTITION_ED));
      int refinements = globalCoarseGridLevel_ - partitionData->getLevel();
      elInfo->getElement()->setMark(max(0, refinements));
      elInfo = stack.traverseNext(elInfo);
    }

    refinementManager->refineMesh(mesh);
  }

  void ParallelProblem::coarsenOutOfPartition(AdaptInfo *adaptInfo)
  {
    Flag meshCoarsened = 1;    
    while (meshCoarsened.getFlags() != 0) {
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      while (elInfo) {
	Element *element = elInfo->getElement();
	PartitionElementData *partitionData = 
	  dynamic_cast<PartitionElementData*>
	  (element->getElementData(PARTITION_ED));
	if (partitionData->getPartitionStatus() == OUT) {
	  int mark = min(0, -partitionData->getLevel() + globalCoarseGridLevel_);
	  element->setMark(mark);
	}
	elInfo = stack.traverseNext(elInfo);
      }
      meshCoarsened = coarseningManager->coarsenMesh(mesh);
    }
    mpiComm.Barrier();
  }

  void ParallelProblem::exchangeMeshStructureCodes(MeshStructure *structures)
  {
    // every process creates a mesh structure code from its mesh.
    structures[mpiRank].init(mesh);
    const std::vector<unsigned long int>& myCode = structures[mpiRank].getCode();

    // broadcast code sizes
    int *codeSize = new int[mpiSize];
    int tmp = static_cast<int>(myCode.size());
    mpiComm.Allgather(&tmp, 1, MPI_INT, codeSize, 1, MPI_INT);
    if (debugMode) {
      // send code sizes also to debug server
      MPI::COMM_WORLD.Gather(&tmp, 1, MPI_INT, NULL, 1, MPI_INT, 0);
    }

    // broadcast number of elements
    int *elements = new int[mpiSize];
    tmp = structures[mpiRank].getNumElements();
    mpiComm.Allgather(&tmp, 1, MPI_INT, elements, 1, MPI_INT);
    if (debugMode) {
      // send number of elements also to debug server
      MPI::COMM_WORLD.Gather(&tmp, 1, MPI_INT, NULL, 1, MPI_INT, 0);
    }

    // broadcast codes
    int *codeOffset = new int[mpiSize];
    int codeSizeSum = 0;
    for (int rank = 0; rank < mpiSize; rank++) {
      codeOffset[rank] = codeSizeSum;
      codeSizeSum += codeSize[rank];
    }

    unsigned long int *code = new unsigned long int[codeSizeSum];
    unsigned long int *localCode = new unsigned long int[codeSize[mpiRank]];
    unsigned long int *ptr;
    std::vector<unsigned long int>::const_iterator it, end = myCode.end();
  
    for (ptr = localCode, it = myCode.begin(); it != end; ++it, ++ptr)
      *ptr = *it;

    mpiComm.Allgatherv(localCode, codeSize[mpiRank], 
		       MPI_UNSIGNED_LONG, 
		       code, codeSize, codeOffset,
		       MPI_UNSIGNED_LONG);
    if (debugMode) {
      // send codes also to debug server
      MPI::COMM_WORLD.Send(localCode, codeSize[mpiRank],
			   MPI_UNSIGNED_LONG, 0, 100);
    }
    
    for (int rank = 0; rank < mpiSize; rank++) {
      if (rank != mpiRank) {
	std::vector<unsigned long int> remoteCode;
	unsigned long int *ptr;
	unsigned long int *begin = code + codeOffset[rank]; 
	unsigned long int *end = begin + codeSize[rank];
	for (ptr = begin; ptr != end; ++ptr) {
	  remoteCode.push_back(*ptr);
	}
	structures[rank].init(remoteCode, elements[rank]);
      }
    }

    delete [] elements;
    delete [] code;
    delete [] localCode;
    delete [] codeOffset;
    delete [] codeSize;    
  }

  void ParallelProblem::synchronizeMeshes(AdaptInfo *adaptInfo)
  {
    FUNCNAME("ParallelProblem::synchronizeMeshes()");

    MeshStructure *structures = new MeshStructure[mpiSize];

    // build composite mesh structure
    exchangeMeshStructureCodes(structures);

    // merge codes
    for (int rank = 0; rank < mpiSize; rank++)
      if (rank != mpiRank)
	structures[mpiRank].merge(&structures[rank]);
  
    // build finest mesh on the rank partition
    structures[mpiRank].fitMeshToStructure(mesh,
					   refinementManager,
					   true);

    delete [] structures;
  }


  bool ParallelProblem::writeElement(ElInfo *elInfo)
  {
    Element *element = elInfo->getElement();
    PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
      (element->getElementData(PARTITION_ED));
    TEST_EXIT_DBG(partitionData)("no partition data\n");
    PartitionStatus status = partitionData->getPartitionStatus();
    if (status == IN) 
      return true;
    else
      return false;
  }

  void ParallelProblem::exchangeRankSolutions(AdaptInfo *adaptInfo,
					      Mesh *workMesh,
					      std::vector<DOFVector<double>*> rankSolutions)
  {
    FUNCNAME("ParallelProblem::exchangeRankSolutions()");

    ParallelProblem::fillVertexPartitions(localCoarseGridLevel_, 1, true, overlapDistance_);
    overlapDistance_.clear();

    const FiniteElemSpace *feSpace = rankSolutions[0]->getFeSpace();
    int dim = workMesh->getDim();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int numFcts = basFcts->getNumber();
    DegreeOfFreedom *coarseDOFs = new DegreeOfFreedom[numFcts];
    DegreeOfFreedom *fineDOFs = new DegreeOfFreedom[numFcts];
    DOFAdmin *admin = feSpace->getAdmin();

    std::vector<std::vector<DegreeOfFreedom> > sendOrder(mpiSize);
    std::vector<std::vector<DegreeOfFreedom> > recvOrder(mpiSize);

    elementPartitions_.clear();

    int elementPartition = -1;
    Element *coarseElement = NULL;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(workMesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      Element *element = elInfo->getElement();

      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));

      if (partitionData) {
	if (partitionData->getLevel() == 0) {
	  elementPartition = partitionVec[element->getIndex()];
	}

	PartitionStatus status = partitionData->getPartitionStatus();

 	if (status != OUT) {
	  if (partitionData->getLevel() == localCoarseGridLevel_) {
	    basFcts->getLocalIndices(element, admin, coarseDOFs);

	    // collect other partitions element belongs to
	    for (int i = 0; i < dim + 1; i++) {
	      std::set<int>::iterator setBegin = vertexPartitions[coarseDOFs[i]].begin();
	      std::set<int>::iterator setEnd = vertexPartitions[coarseDOFs[i]].end();
	      for (std::set<int>::iterator setIt = setBegin; setIt != setEnd; ++setIt) {
		elementPartitions_[element].insert(*setIt);
	      }	
	    }

	    coarseElement = element;
	  }


	  if (element->isLeaf()) {
	    basFcts->getLocalIndices(element, admin, fineDOFs);

	    for (int i = 0; i < numFcts; i++) {
	      if (status == OVERLAP) {
		// send dofs
		sendOrder[elementPartition].push_back(fineDOFs[i]);
	      } 
	      if (status == IN) {
		// recv dofs
		TEST_EXIT(elementPartition == mpiRank)("???\n");
		std::set<int>::iterator setBegin = elementPartitions_[coarseElement].begin();
		std::set<int>::iterator setEnd = elementPartitions_[coarseElement].end();
		for (std::set<int>::iterator setIt = setBegin; setIt != setEnd; ++setIt) {
		  if (*setIt != mpiRank) {
 		    recvOrder[*setIt].push_back(fineDOFs[i]);
		  }
		}
	      }
	    }
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }

    // create send and recv buffers and fill send buffers
    DOFVector<double> *solution = rankSolutions[mpiRank];

    std::map<int, double*> sendBuffer;
    std::map<int, double*> recvBuffer;
    std::map<int, int> sendBufferSize;
    std::map<int, int> recvBufferSize;

    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	int sendSize = static_cast<int>(sendOrder[partition].size());
	int recvSize = static_cast<int>(recvOrder[partition].size());

	sendBufferSize[partition] = sendSize;	
	recvBufferSize[partition] = recvSize;
	if (sendSize > 0) {
	  sendBuffer[partition] = new double[sendSize];
	  std::vector<DegreeOfFreedom>::iterator dofIt;
	  dofIt = sendOrder[partition].begin();
	  double *bufferIt, *bufferBegin, *bufferEnd;
	  bufferBegin = sendBuffer[partition];
	  bufferEnd = bufferBegin + sendSize;
	  for (bufferIt = bufferBegin; bufferIt < bufferEnd; ++bufferIt, ++dofIt)
	    *bufferIt = (*solution)[*dofIt];
	}
	if (recvSize > 0)
	  recvBuffer[partition] = new double[recvSize];
      }
    }

    // non-blocking sends
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (sendBufferSize[partition] > 0) {
	  mpiComm.Isend(sendBuffer[partition],
			sendBufferSize[partition],
			MPI_DOUBLE,
			partition,
			0);
	}
      }
    }    
   
    // blocking recieves
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (recvBufferSize[partition] > 0) {
	  mpiComm.Recv(recvBuffer[partition],
		       recvBufferSize[partition],
		       MPI_DOUBLE,
		       partition,
		       0);
	}
      }
    }    

    // wait for end of communication
    mpiComm.Barrier();

    // copy values into rank solutions
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	std::vector<DegreeOfFreedom>::iterator dofIt = recvOrder[partition].begin();
	for (int i = 0; i < recvBufferSize[partition]; i++) {
	  (*(rankSolutions[partition]))[*dofIt] = recvBuffer[partition][i];
	  ++dofIt;
	}
      }
    }    

    // free send and recv buffers
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (sendBufferSize[partition] > 0)
	  delete [] sendBuffer[partition];
	if (recvBufferSize[partition] > 0)
	  delete [] recvBuffer[partition];
      }
    }    

    delete [] coarseDOFs;
    delete [] fineDOFs;
  }

  void ParallelProblem::exchangeDOFVector(AdaptInfo *adaptInfo,
					  DOFVector<double> *values)
  {
    partitioner->fillLeafPartitionVec(&oldPartitionVec, &oldPartitionVec);
    partitioner->fillLeafPartitionVec(&partitionVec, &partitionVec);

    // === get send and recieve orders ===
    std::vector<std::vector<DegreeOfFreedom> > sendOrder;
    std::vector<std::vector<DegreeOfFreedom> > recvOrder;
    sendOrder.resize(mpiSize);
    recvOrder.resize(mpiSize);

    const FiniteElemSpace *feSpace = values->getFeSpace();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int numFcts = basFcts->getNumber();
    DegreeOfFreedom *dofs = new DegreeOfFreedom[numFcts];
    DOFAdmin *admin =  feSpace->getAdmin();

    Mesh *mesh = feSpace->getMesh();
    TraverseStack stack;
    ElInfo *elInfo;

    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while(elInfo) {
      Element *element = elInfo->getElement();
      int index = element->getIndex();
      int oldPartition = oldPartitionVec[index];
      int newPartition = partitionVec[index];

      if (oldPartition != newPartition) {
	// get dof indices
	basFcts->getLocalIndices(element, admin, dofs);

	if (oldPartition == mpiRank) {
	  // send element values to new partition
	  for (int i = 0; i < numFcts; i++)
	    sendOrder[newPartition].push_back(dofs[i]);
	}
	if (newPartition == mpiRank) {
	  // recv element values from old partition
	  for (int i = 0; i < numFcts; i++)
	    recvOrder[oldPartition].push_back(dofs[i]);
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    delete [] dofs;

    // === create send and recv buffers and fill send buffers ===
    std::map<int, double*> sendBuffer;
    std::map<int, double*> recvBuffer;
    std::map<int, int> sendBufferSize;
    std::map<int, int> recvBufferSize;

    int partition;
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	int sendSize = static_cast<int>(sendOrder[partition].size());
	int recvSize = static_cast<int>(recvOrder[partition].size());
      
	sendBufferSize[partition] = sendSize;	
	recvBufferSize[partition] = recvSize;
	if (sendSize > 0) {
	  sendBuffer[partition] = new double[sendSize];
	  std::vector<DegreeOfFreedom>::iterator dofIt;
	  dofIt = sendOrder[partition].begin();
	  double *bufferIt, *bufferBegin, *bufferEnd;
	  bufferBegin = sendBuffer[partition];
	  bufferEnd = bufferBegin + sendSize;
	  for (bufferIt = bufferBegin; bufferIt < bufferEnd; ++bufferIt, ++dofIt)
	    *bufferIt = (*values)[*dofIt];
	}
	if (recvSize > 0)
	  recvBuffer[partition] = new double[recvSize];
      }
    }

    // === non-blocking sends ===
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (sendBufferSize[partition] > 0) {
	  mpiComm.Isend(sendBuffer[partition],
			  sendBufferSize[partition],
			  MPI_DOUBLE,
			  partition,
			  0);
	}
      }
    }    
    
    // === blocking receives ===
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (recvBufferSize[partition] > 0) {
	  mpiComm.Recv(recvBuffer[partition],
			 recvBufferSize[partition],
			 MPI_DOUBLE,
			 partition,
			 0);
	}
      }
    }    

    // === wait for end of MPI communication ===
    mpiComm.Barrier();

    // === copy received values into DOFVector ===
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	std::vector<DegreeOfFreedom>::iterator dofIt = recvOrder[partition].begin();
	for (int i = 0; i < recvBufferSize[partition]; i++) {
	  (*values)[*dofIt] = recvBuffer[partition][i];
	  ++dofIt;
	}
      }
    }    

    // === free send and receive buffers ===
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (sendBufferSize[partition] > 0)
	  delete [] sendBuffer[partition];
	if (recvBufferSize[partition] > 0)
	  delete [] recvBuffer[partition];
      }
    }
  }

  void ParallelProblem::buildGlobalSolution(AdaptInfo *adaptInfo,
					    std::vector<DOFVector<double>*> rankSolutions,
					    DOFVector<double> *globalSolution)
  {
    FUNCNAME("ParallelProblem::buildGlobalSolution()");

    const FiniteElemSpace *feSpace = globalSolution->getFeSpace();
    int dim = mesh->getDim();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int numFcts = basFcts->getNumber();
    DegreeOfFreedom *coarseDOFs = new DegreeOfFreedom[numFcts];
    DegreeOfFreedom *fineDOFs = new DegreeOfFreedom[numFcts];
    DOFAdmin *admin = feSpace->getAdmin();
    Lagrange *linearFunctions = Lagrange::getLagrange(dim, 1);

    MSG("Building global solution\n");

    // compute w[DOF][partition]->value
    std::map<DegreeOfFreedom, std::map<int, double> > w;
    std::map<DegreeOfFreedom, std::map<int, double> >::iterator wIt, wBegin, wEnd;

    std::map<DegreeOfFreedom, double> sumW;

    Element *lastCoarseElement = NULL;

    WorldVector<double> worldCoord;
    DimVec<double> baryCoord(dim, NO_INIT);

    std::set<int>::iterator partIt, partBegin, partEnd;

    std::map<DegreeOfFreedom, bool> visited;

    MyDualTraverse dualTraverse(localCoarseGridLevel_);
    ElInfo *elInfo1, *elInfo2;
    ElInfo *large, *small;

    bool cont = dualTraverse.traverseFirst(mesh, mesh,
					   -1, -1,
					   Mesh::CALL_EVERY_EL_PREORDER |
					   Mesh::FILL_COORDS | 
					   Mesh::FILL_DET,
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_COORDS,
					   &elInfo1, &elInfo2,
					   &small, &large);

    while (cont) {
      Element *element1 = elInfo1->getElement();
      Element *element2 = elInfo2->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>
	(element1->getElementData(PARTITION_ED));
#ifdef DEBUG
	if(element1->getElementData()==NULL)
		std::cout<<"even no elementData, id:"<<element1->getIndex()<<std::endl;
	if(partitionData==NULL) std::cout<<"elem_id:"<<element1->getIndex()<<std::endl;
#endif
	TEST_EXIT_DBG(partitionData!=NULL)("PARTITION_ED data not found\n");
      if (partitionData->getPartitionStatus() == IN) {

	// get coarse dofs
	if (element1 != lastCoarseElement) {
	  basFcts->getLocalIndices(element1, admin, coarseDOFs);
	  lastCoarseElement = element1;
	}
      
	if (elementPartitions_[element1].size() > 1) {
	  // get fine dofs
	  basFcts->getLocalIndices(element2, admin, fineDOFs);
      
	  // for all fine DOFs
	  for (int i = 0; i < numFcts; i++) {
	    if (!visited[fineDOFs[i]]) {
	      visited[fineDOFs[i]] = true;

	      elInfo2->coordToWorld(*(basFcts->getCoords(i)), worldCoord);
	      elInfo1->worldToCoord(worldCoord, &baryCoord);

	      // for all coarse vertex DOFs
	      for (int j = 0; j < dim + 1; j++) {
		partBegin = vertexPartitions[coarseDOFs[j]].begin();
		partEnd = vertexPartitions[coarseDOFs[j]].end();
		for (partIt = partBegin; partIt != partEnd; ++partIt) {
		  int partition = *partIt/* - 1*/;
		  double val = (*(linearFunctions->getPhi(j)))(baryCoord);
		  w[fineDOFs[i]][partition] += val;

		  sumW[fineDOFs[i]] += val;
		}
	      }
	    }
	  }
	}
      }
      cont = dualTraverse.traverseNext(&elInfo1, &elInfo2,
				       &small, &large);
    }

    delete [] coarseDOFs;
    delete [] fineDOFs;

    MSG("PU ...\n");

    wBegin = w.begin();
    wEnd = w.end();

    for (wIt = wBegin; wIt != wEnd; ++wIt) {
      DegreeOfFreedom dof = wIt->first;
      (*globalSolution)[dof] = 0.0;
    }
    
    for (wIt = wBegin; wIt != wEnd; ++wIt) {
      DegreeOfFreedom dof = wIt->first;
      std::map<int, double>::iterator partIt, partBegin, partEnd;
      partBegin = wIt->second.begin();
      partEnd = wIt->second.end();
    
      for (partIt = partBegin; partIt != partEnd; ++partIt) {
	int partition = partIt->first;
	double wDOF = partIt->second;
	(*globalSolution)[dof] += wDOF / sumW[dof] * (*(rankSolutions[partition]))[dof];
      }
    }
  }

  double ParallelProblem::errors2map(std::map<int, double> &errVec, 
				     int comp,
				     bool add)
  {
    int elNum = -1;
    double totalErr, error;

    if (!add) 
      errVec.clear();

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 
					 -1,
					 Mesh::CALL_EVERY_EL_PREORDER);
    totalErr = 0.0;
    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()) {
	  error = element->getEstimation(comp);
	  errVec[elNum] += error;
	  totalErr += error;
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }

    return totalErr;
  }

  double ParallelProblem::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 ParallelProblem::globalRefinements() 
  {
    if (globalRefinements_ <= 0) 
      return;

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

      if (partitionData && partitionData->getPartitionStatus() == IN) {
	element->setMark(globalRefinements_);
      }
      
      elInfo = stack.traverseNext(elInfo);
    }

    refinementManager->refineMesh(mesh);
  }


  void ParallelProblem::createOverlap(int level, int overlap, bool openOverlap,
				      std::map<Element*, int> &overlapDistance)
  {
    int dim = mesh->getDim();

    // === create dual graph (one common node) and prepare breadth-first search ===
    std::map<DegreeOfFreedom, std::vector<Element*> > vertexElements;
    std::queue<Element*> overlapQueue;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
      
      if (partitionData) {
	if (partitionData->getLevel() == level) {
	  for (int i = 0; i < dim + 1; i++) {	    
	    vertexElements[element->getDOF(i, 0)].push_back(element);
	  }

	  if(partitionData->getPartitionStatus() == IN) {
	    overlapDistance[element] = 0;
	    overlapQueue.push(element);
	  } else {
	    overlapDistance[element] = -1; // still unknown
	  }
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }

    // === breadth-first search on dual graph ===
    std::vector<Element*>::iterator it, itBegin, itEnd;

    while (!overlapQueue.empty()) {
      // get first element in queue
      Element *element = overlapQueue.front();
      overlapQueue.pop();

      // get distance
      int distance = overlapDistance[element];
      TEST_EXIT(distance >= 0)("invalid distance\n");

      if (distance >= overlap) 
	continue;
      
      // get all adjacent elements
      for (int i = 0; i < dim + 1; i++) {
	itBegin = vertexElements[element->getDOF(i, 0)].begin();
	itEnd = vertexElements[element->getDOF(i, 0)].end();
	for (it = itBegin; it != itEnd; ++it) {
	  if (overlapDistance[*it] == -1) {
	    // set distance for new member
	    overlapDistance[*it] = distance + 1;
	    // push neighbour to queue
	    overlapQueue.push(*it);
	    // mark as overlap on AMDiS mesh	    
	    PartitionElementData *partitionData = 
	      dynamic_cast<PartitionElementData*>((*it)->getElementData(PARTITION_ED));
	    TEST_EXIT(partitionData)("no partition data\n");
	    partitionData->setPartitionStatus(OVERLAP);
	    partitionData->descend(*it);
	  }
	}
      }
    }
  }

  void ParallelProblem::fillVertexPartitions(int level, int overlap, 
					     bool openOverlap,
					     std::map<Element*, int> &overlapDistance)
  {
    int dim = mesh->getDim();

    // clear partition dof vector
    vertexPartitions.clear();

    // first: partition elements ...
    int partition;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
      if (partitionData) {
	if (partitionData->getLevel() == 0) {
	  partition = partitionVec[element->getIndex()];
	}
	if (partitionData->getLevel() == level) {
	  for (int i = 0; i < dim + 1; i++) {
	    vertexPartitions[element->getDOF(i, 0)].insert(partition);
	  }
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }

    if (overlap > 1 || openOverlap == false) {
      // exchange mesh structure codes
      MeshStructure *structures = new MeshStructure[mpiSize];
      exchangeMeshStructureCodes(structures);

      // merge codes
      for (int rank = 0; rank < mpiSize; rank++) {
	if (rank != mpiRank) {
	  structures[mpiRank].merge(&structures[rank]);
	}
      }
    
      MeshStructure &compositeStructure = structures[mpiRank];
      compositeStructure.reset();

      // get composite indices of local overlap elements
      std::map<int, Element*> indexElement;
      std::vector<int> innerOverlapElements; // not at open overlap boundary

      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
      while (elInfo) {
	Element *element = elInfo->getElement();
	PartitionElementData *partitionData = 
	  dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
	
	if (partitionData && partitionData->getLevel() == level) {
	  int compositeI