#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 "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 "mpi.h"
#include <queue>

namespace AMDiS {

  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;
    };
  private:
    int coarseLevel_;
  };

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

  ParallelProblemBase::ParallelProblemBase(ProblemIterationInterface *iterationIF,
					   ProblemTimeInterface *timeIF)
    : iterationIF_(iterationIF),
      timeIF_(timeIF)
  {
    mpiRank_ = MPI::COMM_WORLD.Get_rank();
    mpiSize_ = MPI::COMM_WORLD.Get_size();
  }

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

  ParallelProblem::ParallelProblem(const std::string& name,
				   ProblemIterationInterface *iterationIF,
				   ProblemTimeInterface *timeIF,
				   std::vector<DOFVector<double>*> vectors,
				   Mesh *mesh,
				   RefinementManager *rm,
				   CoarseningManager *cm)
    : ParallelProblemBase(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_);

    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 = GET_MEMORY(double, mpiSize_);
    int *partArray = GET_MEMORY(int, mpiSize_);
    int part;

    //weightSum[mpiRank_] = localWeightSum;

    MPI::COMM_WORLD.Gather(&localWeightSum, 1, MPI_DOUBLE,
			   weightSum, 1, MPI_DOUBLE, 0);

    if(mpiRank_ == 0) {

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

      average /= mpiSize_;

      //MSG("average weight: %f\n", average);

      part = 0;
      for(i = 0; i < mpiSize_; i++) {
	//MSG("weight sum %d: %f\n", i, weightSum[i]);

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

      //MSG("upper threshold initial: %f\n", upperPartThreshold_);
      //MSG("lower threshold initial: %f\n", lowerPartThreshold_);
      //MSG("part initial: %d\n", part);

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

      //MSG("comp time: %f\n", computationTime);
      //MSG("part time: %f\n", partitioningTime);
      //MSG("time quotient: %f\n", computationTime/partitioningTime);
      //MSG("upper threshold final: %f\n", upperPartThreshold_);
      //MSG("lower threshold final: %f\n", lowerPartThreshold_);
      //MSG("part final: %d\n", part);

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

    MPI::COMM_WORLD.Scatter(partArray, 1, MPI_INT,
			    &part, 1, MPI_INT, 0);

    
    //MSG("rank %d: part: %d\n", mpiRank_, part);


    FREE_MEMORY(weightSum, double, mpiSize_);
    FREE_MEMORY(partArray, int, mpiSize_);

    return (part == 1);
  }

  bool ParallelProblem::doBuildGlobalSolution(AdaptInfo *adaptInfo) {
    return true;
//       (puEveryTimestep_ || !timeIF_ ||
// 	    (adaptInfo->getTimestepNumber() % repartitionSteps_ == 0) ||
// 	    adaptInfo->getTime() >= adaptInfo->getEndTime());
  }

  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 i, dim = mesh_->getDim();

    bool finished = (localCoarseGridLevel_ == 0);

    //for(j = globalCoarseGridLevel_; j < localCoarseGridLevel_; j++) {
    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(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(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(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_);
	  //int mark = -partitionData->getLevel();
	  element->setMark(mark);
	}
	elInfo = stack.traverseNext(elInfo);
      }
      meshCoarsened = coarseningManager_->coarsenMesh(mesh_);
    }
    MPI::COMM_WORLD.Barrier();
  }

  void ParallelProblem::exchangeMeshStructureCodes(MeshStructure *structures)
  {
    structures[mpiRank_].init(mesh_);
    const std::vector<unsigned long int>& myCode = structures[mpiRank_].getCode();
    int rank;


    // broadcast code sizes
    int *codeSize = GET_MEMORY(int, mpiSize_);
    int tmp = static_cast<int>(myCode.size());

    MPI::COMM_WORLD.Allgather(&tmp, 1, MPI_INT,
			      codeSize, 1, MPI_INT);
    
    // broadcast number of elements
    int *elements = GET_MEMORY(int, mpiSize_);
    tmp = structures[mpiRank_].getNumElements();
    MPI::COMM_WORLD.Allgather(&tmp, 1, MPI_INT,
			      elements, 1, MPI_INT);

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

    unsigned long int *code = GET_MEMORY(unsigned long int, codeSizeSum);
    unsigned long int *localCode = GET_MEMORY(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;
    }
  
    MPI::COMM_WORLD.Allgatherv(localCode, codeSize[mpiRank_], 
			       MPI_UNSIGNED_LONG, 
			       code, codeSize, codeOffset,
			       MPI_UNSIGNED_LONG);
    
    for(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]);
      }
    }

    // free memory
    FREE_MEMORY(elements, int, mpiSize_);
    FREE_MEMORY(code, unsigned long int, codeSizeSum);
    FREE_MEMORY(localCode, unsigned long int, codeSize[mpiRank_]);
    FREE_MEMORY(codeOffset, int, mpiSize_);
    FREE_MEMORY(codeSize, int, mpiSize_);
    
  }

  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 ===
    structures[mpiRank_].fitMeshToStructure(mesh_,
					    refinementManager_,
					    true);

    DELETE [] structures;

#if 0
    // === count partition elements (only for debug) ===
    int partitionElements = 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 && partitionData->getPartitionStatus() == IN) {
	partitionElements++;
      }
      elInfo = stack.traverseNext(elInfo);
    }
    MSG("rank %d partition elements: %d\n", mpiRank_, partitionElements);
#endif
  }


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

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

    int level = localCoarseGridLevel_, overlap = 1;
    bool openOverlap = true;
    ParallelProblem::fillVertexPartitions(level, overlap, openOverlap, 
 					  overlapDistance_);
    overlapDistance_.clear();


    const FiniteElemSpace *feSpace = rankSolutions[0]->getFESpace();
    TEST_EXIT(feSpace->getMesh() == mesh_)("invalid mesh\n");
    int dim = mesh_->getDim();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int numFcts = basFcts->getNumber();
    DegreeOfFreedom *coarseDOFs = GET_MEMORY(DegreeOfFreedom, numFcts);
    DegreeOfFreedom *fineDOFs = GET_MEMORY(DegreeOfFreedom, numFcts);
    DOFAdmin *admin =  feSpace->getAdmin();
    int partition;

    std::vector<std::vector<DegreeOfFreedom> > sendOrder;
    std::vector<std::vector<DegreeOfFreedom> > recvOrder;
    sendOrder.resize(mpiSize_);
    recvOrder.resize(mpiSize_);

    std::set<int>::iterator setIt, setBegin, setEnd;

    elementPartitions_.clear();

    int elementPartition = -1;
    Element *coarseElement = NULL;
    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) {
	  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++) {
	      setBegin = vertexPartitions_[coarseDOFs[i]].begin();
	      setEnd = vertexPartitions_[coarseDOFs[i]].end();
	      for (setIt = setBegin; setIt != setEnd; ++setIt) {
		elementPartitions_[element].insert(*setIt/* - 1*/);
	      }	
	    }

	    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");
		setBegin = elementPartitions_[coarseElement].begin();
		setEnd = elementPartitions_[coarseElement].end();
		for (setIt = setBegin; setIt != setEnd; ++setIt) {
		  if (*setIt != mpiRank_) {
 		    recvOrder[*setIt].push_back(fineDOFs[i]);
		  }
		}
	      }
	    }
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }


#if 0
    MyDualTraverse dualTraverse(localCoarseGridLevel_);//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) {
      TEST_EXIT(elInfo1 == large && elInfo2 == small)
	("error in dual traverse\n");

      Element *element1 = elInfo1->getElement();
      Element *element2 = elInfo2->getElement();
    
      // get partition status of element2
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element2->getElementData(PARTITION_ED));
      TEST_EXIT(partitionData)("no partition data\n");
      PartitionStatus status = partitionData->getPartitionStatus();

      if(status != OUT) {
	basFcts->getLocalIndices(element1,
				 admin,
				 coarseDOFs);
	basFcts->getLocalIndices(element2,
				 admin,
				 fineDOFs);

	// collect other partitions element2 belongs to
	for(i = 0; i < dim + 1; i++) {
	  setBegin = vertexPartitions_[coarseDOFs[i]].begin();
	  setEnd = vertexPartitions_[coarseDOFs[i]].end();
	  for(setIt = setBegin; setIt != setEnd; ++setIt) {
	    elementPartitions_[element1].insert(*setIt/* - 1*/);
	  }	
	}
      
	int elementPartition = partitionVec_[element1->getIndex()];

	// collect send- and recv-dofs for every rank 
	WorldVector<double> worldCoords;
	for(i = 0; i < numFcts; i++) {
	  elInfo2->coordToWorld(*(basFcts->getCoords(i)), &worldCoords);
	  if (status == OVERLAP) {
	    // send dofs
	    //sortedSendDOFs[elementPartition][worldCoords] = fineDOFs[i];
	    sendOrder[elementPartition].push_back(fineDOFs[i]);
	  } 
	  if (status == IN) {
	    // recv dofs
	    TEST_EXIT(elementPartition == mpiRank_)("???\n");
	    setBegin = elementPartitions_[element1].begin();
	    setEnd = elementPartitions_[element1].end();
	    for(setIt = setBegin; setIt != setEnd; ++setIt) {
	      //sortedRecvDOFs[*setIt][worldCoords] = fineDOFs[i];
	      if(*setIt != mpiRank_) {
		recvOrder[*setIt].push_back(fineDOFs[i]);
	      }
	    }
	  }
	}
      }
      
      cont = dualTraverse.traverseNext(&elInfo1, &elInfo2,
				       &small, &large);
    }
#endif

    // 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 (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] = GET_MEMORY(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] = GET_MEMORY(double, recvSize);
      }
    }

    // non-blocking sends
    for(partition = 0; partition < mpiSize_; partition++) {
      if(partition != mpiRank_) {
	if(sendBufferSize[partition] > 0) {
	  MPI::COMM_WORLD.Isend(sendBuffer[partition],
				sendBufferSize[partition],
				MPI_DOUBLE,
				partition,
				0);
	}
      }
    }    
  
    // blocking recieves
    for(partition = 0; partition < mpiSize_; partition++) {
      if(partition != mpiRank_) {
	if(recvBufferSize[partition] > 0) {
	  MPI::COMM_WORLD.Recv(recvBuffer[partition],
			       recvBufferSize[partition],
			       MPI_DOUBLE,
			       partition,
			       0);
	}
      }
    }    

    // wait for end of communication
    MPI::COMM_WORLD.Barrier();

    // copy values into rank solutions
    for (partition = 0; partition < mpiSize_; partition++) {
      if (partition != mpiRank_) {
	std::vector<DegreeOfFreedom>::iterator dofIt = recvOrder[partition].begin();
	for (i = 0; i < recvBufferSize[partition]; i++) {
	  (*(rankSolutions[partition]))[*dofIt] = 
	    recvBuffer[partition][i];
	  ++dofIt;
	}
      }
    }    
    
    // free send and recv buffers
    for (partition = 0; partition < mpiSize_; partition++) {
      if (partition != mpiRank_) {
	if (sendBufferSize[partition] > 0)
	  FREE_MEMORY(sendBuffer[partition], 
		      double,
		      sendBufferSize[partition]);
	if(recvBufferSize[partition] > 0)
	  FREE_MEMORY(recvBuffer[partition], 
		      double,
		      recvBufferSize[partition]);
      }
    }    
    FREE_MEMORY(coarseDOFs, DegreeOfFreedom, numFcts);
    FREE_MEMORY(fineDOFs, DegreeOfFreedom, numFcts);
  }

#if 0
  void ParallelProblem::exchangeElementMarkings(AdaptInfo *adaptInfo)
  {
    FUNCNAME("ParallelProblem::exchangeElementMarkings()");
    int i;

    partitioner_->fillLeafPartitionVec(&oldPartitionVec_, &oldPartitionVec_);
    partitioner_->fillLeafPartitionVec(&partitionVec_, &partitionVec_);

    //std::map<int, Element*> elementWithIndex;

    // === get send and recieve orders ===
    std::vector<std::vector<Element*> > recvOrder;
    std::vector<std::vector<int> > markings;

    markings.resize(mpiSize_);
    recvOrder.resize(mpiSize_);

    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];
      int mark = element->getMark();

//       if(mark!=0)
// 	MSG("rank %d: index %d  mark %d\n",
// 	    mpiRank_, index,  mark);

      if(oldPartition != newPartition) {
	if(oldPartition == mpiRank_) {
	  markings[newPartition].push_back(mark);
// 	  MSG("rank %d: index %d   %d -> %d   mark %d\n",
// 	      mpiRank_, index, oldPartition, newPartition, mark);
	}
	if(newPartition == mpiRank_) {
	  recvOrder[oldPartition].push_back(element);
// 	  MSG("rank %d: index %d   %d <- %d   mark %d\n",
// 	      mpiRank_, index, oldPartition, newPartition, mark);
	  //elementWithIndex[index] = element;
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }

    // === create send and recv buffers and fill send buffers ===
    std::map<int, int*> sendBuffer;
    std::map<int, int*> 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>(markings[partition].size());
	int recvSize = static_cast<int>(recvOrder[partition].size());
      
	sendBufferSize[partition] = sendSize;	
	recvBufferSize[partition] = recvSize;
	if(sendSize > 0) {
	  sendBuffer[partition] = GET_MEMORY(int, sendSize);
	  std::vector<int>::iterator it;
	  it = markings[partition].begin();
	  int *bufferIt, *bufferBegin, *bufferEnd;
	  bufferBegin = sendBuffer[partition];
	  bufferEnd = bufferBegin + sendSize;
	  for(bufferIt = bufferBegin; 
	      bufferIt < bufferEnd; 
	      ++bufferIt, ++it) 
	  {
	    *bufferIt = *it;
	  }
	}
	if(recvSize > 0)
	  recvBuffer[partition] = GET_MEMORY(int, recvSize);
      }
    }

    // === non-blocking sends ===
    for(partition = 0; partition < mpiSize_; partition++) {
      if(partition != mpiRank_) {
	if(sendBufferSize[partition] > 0) {
	  MPI::COMM_WORLD.Isend(sendBuffer[partition],
				sendBufferSize[partition],
				MPI_INT,
				partition,
				0);
	}
      }
    }    
    
    // === blocking receives ===
    for(partition = 0; partition < mpiSize_; partition++) {
      if(partition != mpiRank_) {
	if(recvBufferSize[partition] > 0) {
	  MPI::COMM_WORLD.Recv(recvBuffer[partition],
			       recvBufferSize[partition],
			       MPI_INT,
			       partition,
			       0);
	}
      }
    }    

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

    // === copy received values into elements ===
    for(partition = 0; partition < mpiSize_; partition++) {
      if(partition != mpiRank_) {
	std::vector<Element*>::iterator it = recvOrder[partition].begin();
	for(i = 0; i < recvBufferSize[partition]; i++) {
	  (*it)->setMark(recvBuffer[partition][i]);
// 	  MSG(" *** rank %d, index %d from %d mark %d\n",
// 	      mpiRank_, (*it)->getIndex(), i, recvBuffer[partition][i]);
	  ++it;
	}
      }
    }    

    // === free send and receive buffers ===
    for(partition = 0; partition < mpiSize_; partition++) {
      if(partition != mpiRank_) {
	if(sendBufferSize[partition] > 0)
	  FREE_MEMORY(sendBuffer[partition], 
		      int,
		      sendBufferSize[partition]);
	if(recvBufferSize[partition] > 0)
	  FREE_MEMORY(recvBuffer[partition], 
		      int,
		      recvBufferSize[partition]);
      }
    }
  }
#endif

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

    int i;
    const FiniteElemSpace *feSpace = values->getFESpace();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int numFcts = basFcts->getNumber();
    DegreeOfFreedom *dofs = GET_MEMORY(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_) {
	  for(i = 0; i < numFcts; i++) {
	    // send element values to new partition
	    sendOrder[newPartition].push_back(dofs[i]);
	  }
	}
	if(newPartition == mpiRank_) {
	  for(i = 0; i < numFcts; i++) {
	    // recv element values from old partition
	    recvOrder[oldPartition].push_back(dofs[i]);
	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    FREE_MEMORY(dofs, DegreeOfFreedom, numFcts);

    // === 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] = GET_MEMORY(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] = GET_MEMORY(double, recvSize);
      }
    }

//     MSG("rank %d: send %d %d %d %d    recv %d %d %d %d\n", 
// 	mpiRank_,
// 	sendBufferSize[0],
// 	sendBufferSize[1],
// 	sendBufferSize[2],
// 	sendBufferSize[3],
// 	recvBufferSize[0],
// 	recvBufferSize[1],
// 	recvBufferSize[2],
// 	recvBufferSize[3]);

    // === non-blocking sends ===
    for(partition = 0; partition < mpiSize_; partition++) {
      if(partition != mpiRank_) {
	if(sendBufferSize[partition] > 0) {
	  MPI::COMM_WORLD.Isend(sendBuffer[partition],
				sendBufferSize[partition],
				MPI_DOUBLE,
				partition,
				0);
	}
      }
    }    
    
    // === blocking receives ===
    for(partition = 0; partition < mpiSize_; partition++) {
      if(partition != mpiRank_) {
	if(recvBufferSize[partition] > 0) {
	  MPI::COMM_WORLD.Recv(recvBuffer[partition],
			       recvBufferSize[partition],
			       MPI_DOUBLE,
			       partition,
			       0);
	}
      }
    }    

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

    // === copy received values into DOFVector ===
    for(partition = 0; partition < mpiSize_; partition++) {
      if(partition != mpiRank_) {
	std::vector<DegreeOfFreedom>::iterator dofIt = recvOrder[partition].begin();
	for(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)
	  FREE_MEMORY(sendBuffer[partition], 
		      double,
		      sendBufferSize[partition]);
	if(recvBufferSize[partition] > 0)
	  FREE_MEMORY(recvBuffer[partition], 
		      double,
		      recvBufferSize[partition]);
      }
    }
  }

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

    const FiniteElemSpace *feSpace = globalSolution->getFESpace();
    int i, dim = mesh_->getDim();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int numFcts = basFcts->getNumber();
    DegreeOfFreedom *coarseDOFs = GET_MEMORY(DegreeOfFreedom, numFcts);
    DegreeOfFreedom *fineDOFs = GET_MEMORY(DegreeOfFreedom, numFcts);
    DOFAdmin *admin =  feSpace->getAdmin();

    Lagrange *linearFunctions = Lagrange::getLagrange(dim, 1);

    MSG("Building global solution\n");

    int j;

//     DOFVector<double> wtest0(feSpace, "wtest0");
//     wtest0.set(0.0);
//     DOFVector<double> wtest1(feSpace, "wtest1");
//     wtest1.set(0.0);
//     DOFVector<double> wtest2(feSpace, "wtest2");
//     wtest2.set(0.0);
//     DOFVector<double> wtest3(feSpace, "wtest3");
//     wtest3.set(0.0);

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


#if 0    
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh_, -1, 
					 Mesh::CALL_EVERY_EL_PREORDER |
					 Mesh::FILL_COORDS);
    while(elInfo) {
      Element *element = elInfo->getElement();

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

      if(partitionData) {
	PartitionStatus status = partitionData->getPartitionStatus();

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

	  if(element->isLeaf()) {
	    if(elementPartitions_[element].size() > 1) {
	      // get fine dofs
	      basFcts->getLocalIndices(element, admin, fineDOFs);
      
	      // for all fine DOFs
	      for(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(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;
		    }
		  }
		}
	      }
	    }
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }
#endif

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

    bool cont;

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

// 		  if(partition == 0) {
// 		    wtest0[fineDOFs[i]] += val;
// 		  }
// 		  if(partition == 1) {
// 		    wtest1[fineDOFs[i]] += val;
// 		  }
// 		  if(partition == 2) {
// 		    wtest2[fineDOFs[i]] += val;
// 		  }
// 		  if(partition == 3) {
// 		    wtest3[fineDOFs[i]] += val;
// 		  }

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


//     if(mpiRank_ == 0) {
//       MacroWriter::writeMacro(feSpace, "output/wtest0.mesh");
//       ValueWriter::writeValues(&wtest0, "output/wtest0.dat");
//       MacroWriter::writeMacro(feSpace, "output/wtest1.mesh");
//       ValueWriter::writeValues(&wtest1, "output/wtest1.dat");
//       MacroWriter::writeMacro(feSpace, "output/wtest2.mesh");
//       ValueWriter::writeValues(&wtest2, "output/wtest2.dat");
//       MacroWriter::writeMacro(feSpace, "output/wtest3.mesh");
//       ValueWriter::writeValues(&wtest3, "output/wtest3.dat");
//     }

    FREE_MEMORY(coarseDOFs, DegreeOfFreedom, numFcts);
    FREE_MEMORY(fineDOFs, DegreeOfFreedom, numFcts);

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

//     double total_sum = 0.0;
//     MPI::COMM_WORLD.Allreduce(&totalErr,
// 			      &total_sum,
// 			      1,
// 			      MPI_DOUBLE,
// 			      MPI_SUM);
//     total_sum = sqrt(total_sum);

//     MSG("error rank %d = %e (total %e)\n", mpiRank_, totalErr, total_sum);

    return totalErr;
  }

//   void ParallelProblem::writeRankMacroAndValues(DOFVector<double> *vec, 
// 						const char *name, 
// 						double time)
//   {
//     char number[3];
//     sprintf(number, "%d", MPI::COMM_WORLD.Get_rank());
    
//     std::string macroFile(name);
//     macroFile += number;
//     macroFile += ".mesh";

//     ConditionalMacroWriter::writeMacro(vec->getFESpace(),
// 				       const_cast<char*>(macroFile.c_str()),
// 				       time,
// 				       -1,
// 				       Mesh::CALL_LEAF_EL,
// 				       &writeElement);

//     std::string valueFile(name);
//     valueFile += number;
//     valueFile += ".dat";

//     ConditionalValueWriter::writeValues(vec,
// 					const_cast<char*>(valueFile.c_str()),
// 					time,
// 					-1,
// 					Mesh::CALL_LEAF_EL,
// 					&writeElement);
//   }


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

    TraverseStack stack;
    ElInfo *elInfo;

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

    // first: partition elements ...
    int index, partition;
    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) {
	  index = element->getIndex();
	  partition = partitionVec_[index];
	}
	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
      int rank;
      for(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 compositeIndex = compositeStructure.getCurrentElement();
	  indexElement[compositeIndex] = element;
	  if(partitionData->getPartitionStatus() == OVERLAP) {
	    int distance = overlapDistance[element];
	    if(distance < overlap || !openOverlap) {
	      innerOverlapElements.push_back(compositeIndex);
	    } 
	  }
	}
	
	if(element->isLeaf()) {
	  compositeStructure.skipBranch();
	} else {
	  compositeStructure.nextElement();
	}
	elInfo = stack.traverseNext(elInfo);
      }
    
      // === exchange 'inner' overlap elements ===

      // exchange number of overlap elements
      int *numOverlapElements = GET_MEMORY(int, mpiSize_);
      int tmp = static_cast<int>(innerOverlapElements.size());
      MPI::COMM_WORLD.Allgather(&tmp, 1, MPI_INT,
			      numOverlapElements, 1, MPI_INT);
      
      // exchange overlap elements
      int *offset = GET_MEMORY(int, mpiSize_);
      int sum = 0;
      for(rank = 0; rank < mpiSize_; rank++) {
	offset[rank] = sum;
	sum += numOverlapElements[rank];
      }

      int *recvBuffer = GET_MEMORY(int, sum);
      int *sendBuffer = GET_MEMORY(int, numOverlapElements[mpiRank_]);
   
      int *ptr;
      std::vector<int>::iterator elemIt, elemEnd = innerOverlapElements.end();
  
      for(ptr = sendBuffer, elemIt = innerOverlapElements.begin();
	  elemIt != elemEnd; 
	  ++elemIt, ++ptr) 
      {
	*ptr = *elemIt;
      }
  
      MPI::COMM_WORLD.Allgatherv(sendBuffer, numOverlapElements[mpiRank_], 
				 MPI_INT, 
				 recvBuffer, numOverlapElements, offset,
				 MPI_INT);



      // fill vertexPartitions for 'inner' overlap elements
      int el;
      for(rank = 0; rank < mpiSize_; rank++) {
	int numElements = numOverlapElements[rank];
	//       MSG("rank %d overlap elements with %d: %d\n", 
	// 	  mpiRank_, rank, numOverlapElements[rank]);
	int *elements = recvBuffer + offset[rank];
	for(el = 0; el < numElements; el++) {
	  Element *element = indexElement[elements[el]];
	  for(i = 0; i < dim + 1; i++) {
	    vertexPartitions_[element->getDOF(i, 0)].insert(rank/* + 1*/);
	  }
	}
      }

      // free memory
      DELETE [] structures;
      FREE_MEMORY(recvBuffer, int, sum);
      FREE_MEMORY(sendBuffer, int, numOverlapElements[mpiRank_]);
      FREE_MEMORY(numOverlapElements, int, mpiSize_);
    }
  }

  void ParallelProblem::createOverlap(AdaptInfo *adaptInfo)
  {
    int level = localCoarseGridLevel_, overlap = 1;
    bool openOverlap = true;

    ParallelProblem::createOverlap(level, overlap, openOverlap, overlapDistance_);

//     ParallelProblem::fillVertexPartitions(level, overlap, openOverlap, 
//  					  overlapDistance_);
  }

  // =========================================================================
  // ===== class ParallelProblemScal =========================================
  // =========================================================================
  
  ParallelProblemScal::ParallelProblemScal(const std::string& name,
					   ProblemScal *problem,
					   ProblemInstatScal *problemInstat,
					   std::vector<DOFVector<double>*> vectors)
    : ParallelProblem(name,
		      problem,
		      problemInstat,
		      vectors,
		      problem->getMesh(),
		      problem->getRefinementManager(),
		      problem->getCoarseningManager()),
      problem_(problem),
      problemInstat_(problemInstat),
      oldEstimator_(NULL),
      oldMarker_(NULL),
      rankSolution_(mpiSize_),
      usersEstimator_(NULL),
      usersMarker_(NULL)
  {
    int i;
    for(i = 0; i < mpiSize_; i++) {
      rankSolution_[i] = NEW DOFVector<double>(problem_->getFESpace(), "rank solution");
    }

    //vertexPartitions_ = NEW std::map<DegreeOfFreedom, std::set<int> >(problem_->getFESpace(), 
    //						  "partition dof"); 

    if(problemInstat_) {
      dofVectors_.push_back(problemInstat_->getOldSolution());
    }
  }

  ParallelProblemScal::~ParallelProblemScal()
  {
    int i;
    for(i = 0; i < mpiSize_; i++) {
      DELETE rankSolution_[i];
    }
    //DELETE vertexPartitions_;
  }


  void ParallelProblemScal::initParallelization(AdaptInfo *adaptInfo)
  {
    if (mpiSize_ > 1) {
      clock_t partitioningStart = clock();

      partitioner_->createPartitionData();
      setElemWeights(adaptInfo);
      partitionMesh(adaptInfo);

      globalRefineOutOfPartition(adaptInfo);

      refineOverlap(adaptInfo);
      createOverlap(adaptInfo);

      globalRefinements();

      oldEstimator_ = problem_->getEstimator();
      oldMarker_ = problem_->getMarker();

      ConditionalEstimator *condEstimator = NULL;
      if (usersEstimator_) {
	problem_->setEstimator(usersEstimator_);
      } else {
	condEstimator = NEW ConditionalEstimator(oldEstimator_);
	problem_->setEstimator(condEstimator);
      }

      if (usersMarker_) {
	problem_->setMarker(usersMarker_);
      } else {
	ConditionalMarker *newMarker = NEW ConditionalMarker("parallel marker",
							     -1,
							     oldMarker_,
							     globalCoarseGridLevel_,
							     localCoarseGridLevel_);
	problem_->setMarker(newMarker);
      }

      if (mpiRank_ == 0) {
	clock_t partitioningEnd = clock();
	partitioningTime = TIME_USED(partitioningStart, 
				     partitioningEnd);
	computationStart = partitioningEnd;
      }

      // modify file writers
      char number[10];
      sprintf(number, "%d", MPI::COMM_WORLD.Get_rank());
      ::std::list<FileWriterInterface*> fileWriters = problem_->getFileWriterList();
      ::std::list<FileWriterInterface*>::iterator fwIt, fwBegin, fwEnd;
      fwBegin = fileWriters.begin();
      fwEnd = fileWriters.end();
      for(fwIt = fwBegin; fwIt != fwEnd; ++fwIt) {
	(*fwIt)->setFilename((*fwIt)->getFilename() + "_proc" + 
			     ::std::string(number) + "_");
	(*fwIt)->setTraverseProperties(-1, 0, elementInPartition);
      }
    }
  }

  void ParallelProblemScal::exitParallelization(AdaptInfo *adaptInfo)
  {
    if (mpiSize_ > 1) {
      ParallelProblem::exitParallelization(adaptInfo);

      if (!timeIF_) problem_->writeFiles(adaptInfo, true);

      partitioner_->deletePartitionData();

      if (!usersEstimator_) DELETE problem_->getEstimator();
      if (!usersMarker_) DELETE problem_->getMarker();

      problem_->setEstimator(oldEstimator_);
      problem_->setMarker(oldMarker_);    
    }
  }

  int ParallelProblemScal::getNumProblems() { return 1; }

  ProblemStatBase *ParallelProblemScal::getProblem(int number) 
  { 
    return problem_; 
  }

  ProblemStatBase *ParallelProblemScal::getProblem(const std::string& name)
  { 
    return problem_;
  }


  void ParallelProblemScal::exchangeDOFVectors(AdaptInfo *adaptInfo)
  {
    std::vector<DOFVector<double>*>::iterator it, itBegin, itEnd;
    itBegin = dofVectors_.begin();
    itEnd = dofVectors_.end();
    for(it = itBegin; it != itEnd; ++it) {
      ParallelProblem::exchangeDOFVector(adaptInfo, *it);
    }
  }

  void ParallelProblemScal::exchangeRankSolutions(AdaptInfo *adaptInfo)
  {
    rankSolution_[mpiRank_]->copy(*(problem_->getSolution()));
    ParallelProblem::exchangeRankSolutions(adaptInfo,
					   rankSolution_);
  }

  void ParallelProblemScal::buildGlobalSolution(AdaptInfo *adaptInfo)
  {
    ParallelProblem::buildGlobalSolution(adaptInfo,
					 rankSolution_,
					 problem_->getSolution());

//     ParallelProblem::writeRankMacroAndValues(problem_->getSolution(),
// 					     "output/debug_rank");


//    {      
//       char number[3];
//       sprintf(number, "%d", MPI::COMM_WORLD.Get_rank());
//       std::string meshfile = "output/debug_" + std::string(number) + ".mesh";
//       std::string datfile = "output/debug_" + std::string(number) + ".dat";

//       MacroWriter::writeMacro(problem_->getFESpace(), (char*)meshfile.c_str());
//       ValueWriter::writeValues(problem_->getSolution(), (char*)datfile.c_str());

//       if(mpiRank_ == 0) {
// 	std::map<int, int> finePartitionVec;
// 	partitioner_->fillLeafPartitionVec(&partitionVec_,
// 					   &finePartitionVec);
//  	ElementFileWriter elemFW("partition",
//  				 problem_->getMesh(),
//  				 problem_->getFESpace(),
//  				 finePartitionVec);
//  	elemFW.writeFiles(NULL, true);
	
// 	bool wait=true;
// 	while(wait) {}
//       }

//     }
  }

  void ParallelProblemScal::coarsenOutOfPartition(AdaptInfo *adaptInfo)
  {
    bool coarsenAllowed = adaptInfo->isCoarseningAllowed(0);
    adaptInfo->allowCoarsening(true, 0);
    ParallelProblem::coarsenOutOfPartition(adaptInfo);
    adaptInfo->allowCoarsening(coarsenAllowed, 0);

//     int level = localCoarseGridLevel_, overlap = 1;
//     bool openOverlap = true;

//     problem_->getMesh()->dofCompress();

//     ParallelProblem::fillVertexPartitions(level, overlap, openOverlap, 
// 					  overlapDistance_);

//     overlapDistance_.clear();
//     if(mpiRank_ == 0) { 
//       problem_->writeFiles(adaptInfo, true);
//       bool wait=true;
//       while(wait) {}
//     }

//     ParallelProblem::writeRankMacroAndValues(problem_->getSolution(),
// 					     "output/debug_rank",
// 					     0.0);

//     if(mpiRank_ == 0) {
//       std::map<int, int> finePartitionVec;
//       partitioner_->fillLeafPartitionVec(&partitionVec_,
// 					 &finePartitionVec);
//       ElementFileWriter elemFW("partition",
// 			       problem_->getMesh(),
// 			       problem_->getFESpace(),
// 			       finePartitionVec);
//       elemFW.writeFiles(NULL, true);
//       bool wait = true;
//       while(wait) {}
//     }
  }


  // =========================================================================
  // ===== class ParallelProblemVec ==========================================
  // =========================================================================
  
  ParallelProblemVec::ParallelProblemVec(const std::string& name,
					 ProblemVec *problem,
					 ProblemInstatVec *problemInstat,
					 std::vector<DOFVector<double>*> vectors)
    : ParallelProblem(name,
		      problem,
		      problemInstat,
		      vectors,
		      problem->getMesh(0),
		      problem->getRefinementManager(0),
		      problem->getCoarseningManager(0)),
      problem_(problem),
      problemInstat_(problemInstat)
  {
    numComponents_ = problem_->getNumComponents();

    std::vector<FiniteElemSpace*> feSpaces(numComponents_);
    int i, j;
    for(i = 0; i < numComponents_; i++) {
      feSpaces[i] = problem_->getFESpace(i);
    }
    rankSolution_.resize(mpiSize_);
    for(i = 0; i < mpiSize_; i++) {
      rankSolution_[i] = NEW SystemVector("rank solution", feSpaces, numComponents_);
      for(j = 0; j < numComponents_; j++) {
	rankSolution_[i]->setDOFVector(j, NEW DOFVector<double>(feSpaces[j], 
								"rank solution"));
      }
    }
    if(problemInstat_) {
      for(i = 0; i < numComponents_; i++) {
	dofVectors_.push_back(problemInstat_->getOldSolution()->getDOFVector(i));
      }
    }
  }

  ParallelProblemVec::~ParallelProblemVec()
  {
    int i, j;
    for(i = 0; i < mpiSize_; i++) {
      for(j = 0; j < numComponents_; j++) {
	DELETE rankSolution_[i]->getDOFVector(j);
      }
      DELETE rankSolution_[i];
    }
  }

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

    if(mpiSize_ > 1) {
      int i;
      partitioner_->createPartitionData();
      setElemWeights(adaptInfo);
      partitionMesh(adaptInfo);

      refineOverlap(adaptInfo);
      createOverlap(adaptInfo);

      oldEstimator_ = problem_->getEstimator();
      oldMarker_ = problem_->getMarker();

      std::vector<ConditionalEstimator*> condEstimator;
      if(static_cast<int>(usersEstimator_.size()) == numComponents_) {
	problem_->setEstimator(usersEstimator_);
      } else {
	for(i = 0; i < numComponents_; i++) {
	  condEstimator.push_back(NEW ConditionalEstimator(oldEstimator_[i]));
	  problem_->setEstimator(condEstimator[i], i);
	}
      }

      if(static_cast<int>(usersMarker_.size()) == numComponents_) {
	for(i = 0; i < numComponents_; i++) {
	  problem_->setMarker(usersMarker_[i], i);
	}
      } else {
	TEST_EXIT(static_cast<int>(condEstimator.size()) == numComponents_)
	  ("use conditional marker only together with conditional estimator\n");
	for(i = 0; i < numComponents_; i++) {
	  ConditionalMarker *newMarker = 
	    NEW ConditionalMarker("parallel marker",
				  i,
				  oldMarker_[i],
				  globalCoarseGridLevel_,
				  localCoarseGridLevel_);
	  problem_->setMarker(newMarker, i);
	}
      }

      // modify file writers
      char number[10];
      sprintf(number, "%d", MPI::COMM_WORLD.Get_rank());
      ::std::list<FileWriterInterface*> fileWriters = problem_->getFileWriterList();
      ::std::list<FileWriterInterface*>::iterator fwIt, fwBegin, fwEnd;
      fwBegin = fileWriters.begin();
      fwEnd = fileWriters.end();
      for(fwIt = fwBegin; fwIt != fwEnd; ++fwIt) {
	(*fwIt)->setFilename((*fwIt)->getFilename() + "_proc" + 
			     ::std::string(number) + "_");
	(*fwIt)->setTraverseProperties(-1, 0, elementInPartition);
      }    
    }
  }

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

    if(mpiSize_ > 1) {
      ParallelProblem::exitParallelization(adaptInfo);

      if(!timeIF_) problem_->writeFiles(adaptInfo, true);

      partitioner_->deletePartitionData();

      int i;
      for(i = 0; i < numComponents_; i++) {
	if(static_cast<int>(usersEstimator_.size()) == numComponents_) 
	  DELETE problem_->getEstimator(i);
	if(static_cast<int>(usersMarker_.size()) == numComponents_) 
	  DELETE problem_->getMarker(i);

	problem_->setEstimator(oldEstimator_[i], i);
	problem_->setMarker(oldMarker_[i], i);    
      }

      usersEstimator_.resize(0);
      usersMarker_.resize(0);
    }
  }

  void ParallelProblemVec::exchangeDOFVectors(AdaptInfo *adaptInfo) 
  {
    std::vector<DOFVector<double>*>::iterator it, itBegin, itEnd;
    itBegin = dofVectors_.begin();
    itEnd = dofVectors_.end();
    for(it = itBegin; it != itEnd; ++it) {
      ParallelProblem::exchangeDOFVector(adaptInfo, *it);
    }
  }

  void ParallelProblemVec::exchangeRankSolutions(AdaptInfo *adaptInfo) 
  {
    rankSolution_[mpiRank_]->copy(*(problem_->getSolution()));

    std::vector<DOFVector<double>*> rankSol(mpiSize_);

    for (int i = 0; i < numComponents_; i++) {
      for (int j = 0; j < mpiSize_; j++) {
	rankSol[j] = rankSolution_[j]->getDOFVector(i);
      }

      ParallelProblem::exchangeRankSolutions(adaptInfo,
					     rankSol);
    }
  }

  void ParallelProblemVec::buildGlobalSolution(AdaptInfo *adaptInfo) 
  {
    std::vector<DOFVector<double>*> rankSol(mpiSize_);

    int i, j;
    for(i = 0; i < numComponents_; i++) {
      for(j = 0; j < mpiSize_; j++) {
	rankSol[j] = rankSolution_[j]->getDOFVector(i);
      }

      ParallelProblem::buildGlobalSolution(adaptInfo,
					   rankSol,
					   problem_->getSolution()->getDOFVector(i));
    }
  }

  int ParallelProblemVec::getNumProblems() { return 1; }

  ProblemStatBase *ParallelProblemVec::getProblem(int number) 
  { 
    return problem_; 
  }

  ProblemStatBase *ParallelProblemVec::getProblem(const std::string& name)
  { 
    return problem_;
  }

  void ParallelProblemVec::coarsenOutOfPartition(AdaptInfo *adaptInfo)
  {
    std::vector<bool> coarsenAllowed(numComponents_);
    int i;
    for(i = 0; i < numComponents_; i++) {
      coarsenAllowed[i] = adaptInfo->isCoarseningAllowed(i);
      adaptInfo->allowCoarsening(true, i);
    }

    ParallelProblem::coarsenOutOfPartition(adaptInfo);
    for(i = 0; i < numComponents_; i++) {
      adaptInfo->allowCoarsening(coarsenAllowed[i], i);
    }
  }

}