ParallelProblem.cc 48.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#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)
  {
33
34
35
36
37
38
39
    PartitionElementData *elementData = dynamic_cast<PartitionElementData*>
      (elInfo->getElement()->getElementData(PARTITION_ED));
    if (elementData && elementData->getPartitionStatus() == IN) {
      return true;
    } else {
      return false;
    }
40
41
  }

42
43
44
45
46
47
48
49
50
51
52
  bool elementInPartitionDbg(ElInfo *elInfo) 
  {
    // In debug mode, the first partition has to write the whole domain.
    if (MPI::COMM_WORLD.Get_rank() == 0) {
      return true;
    }

    return elementInPartition(elInfo);
  }
      

53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
  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)
  {
86
87
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
  }

  // =========================================================================
  // ===== 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),
103
104
105
      mesh(mesh),
      refinementManager(rm),
      coarseningManager(cm),
106
107
108
109
110
111
112
113
114
115
116
      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),
117
118
      repartTimeFactor_(10.0),
      debugMode(0)
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
  {
    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");

135
    partitioner_ = NEW ParMetisPartitioner(mesh);
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

    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_;
156
157
158
159
160
161

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

    if (debugMode) {
      dbgMesh = NEW Mesh(mesh->getName(), mesh->getDim());
    }
162
163
164
165
166
  }

  ParallelProblem::~ParallelProblem() 
  {
    DELETE partitioner_;
167
168
169
170

    if (debugMode) {
      DELETE dbgMesh;
    }
171
172
173
174
175
176
  }

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

177
178
179
    double *weightSum = GET_MEMORY(double, mpiSize);
    int *partArray = GET_MEMORY(int, mpiSize);
    int part = 0;
180
181
182
183

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

184
    if (mpiRank == 0) {
185
186

      double average = 0.0;
187
      for (int i = 0; i < mpiSize; i++) {
188
189
190
	average += weightSum[i];
      }

191
      average /= mpiSize;
192

193
194
      for (int i = 0; i < mpiSize; i++) {
	if ((weightSum[i] / average) > upperPartThreshold_) {
195
196
197
	  part = 1;
	  break;
	}
198
	if ((weightSum[i] / average) < lowerPartThreshold_) {
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
	  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_);
	}
      }

228
      for (int i = 0; i < mpiSize; i++) {
229
230
231
232
233
234
235
	partArray[i] = part;
      }      
    }

    MPI::COMM_WORLD.Scatter(partArray, 1, MPI_INT,
			    &part, 1, MPI_INT, 0);
    
236
237
    FREE_MEMORY(weightSum, double, mpiSize);
    FREE_MEMORY(partArray, int, mpiSize);
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262

    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)
  {
263
    int i, dim = mesh->getDim();
264
265
266
267
268
269
270
271
272

    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;
273
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
      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;
301
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
      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);
      }
330
      if(marked) refinementManager->refineMesh(mesh);
331
332
333
334
335
336
    }
  }

  void ParallelProblem::globalRefineOutOfPartition(AdaptInfo *adaptInfo)
  {
    TraverseStack stack;
337
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
338
339
340
341
342
343
344
345
    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);
    }

346
    refinementManager->refineMesh(mesh);
347
348
349
350
351
352
353
  }

  void ParallelProblem::coarsenOutOfPartition(AdaptInfo *adaptInfo)
  {
    Flag meshCoarsened = 1;    
    while(meshCoarsened.getFlags() != 0) {
      TraverseStack stack;
354
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
355
356
357
358
359
360
361
362
363
364
365
366
      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);
      }
367
      meshCoarsened = coarseningManager->coarsenMesh(mesh);
368
369
370
371
372
373
    }
    MPI::COMM_WORLD.Barrier();
  }

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

    // broadcast code sizes
379
    int *codeSize = GET_MEMORY(int, mpiSize);
380
381
382
383
384
385
    int tmp = static_cast<int>(myCode.size());

    MPI::COMM_WORLD.Allgather(&tmp, 1, MPI_INT,
			      codeSize, 1, MPI_INT);
    
    // broadcast number of elements
386
387
    int *elements = GET_MEMORY(int, mpiSize);
    tmp = structures[mpiRank].getNumElements();
388
389
390
391
    MPI::COMM_WORLD.Allgather(&tmp, 1, MPI_INT,
			      elements, 1, MPI_INT);

    // broadcast codes
392
    int *codeOffset = GET_MEMORY(int, mpiSize);
393
    int codeSizeSum = 0;
394
    for (int rank = 0; rank < mpiSize; rank++) {
395
396
397
398
399
      codeOffset[rank] = codeSizeSum;
      codeSizeSum += codeSize[rank];
    }

    unsigned long int *code = GET_MEMORY(unsigned long int, codeSizeSum);
400
    unsigned long int *localCode = GET_MEMORY(unsigned long int, codeSize[mpiRank]);
401
402
403
404
   
    unsigned long int *ptr;
    std::vector<unsigned long int>::const_iterator it, end = myCode.end();
  
405
406
407
    for (ptr = localCode, it = myCode.begin();
	 it != end; 
	 ++it, ++ptr) 
408
409
410
411
    {
      *ptr = *it;
    }
  
412
    MPI::COMM_WORLD.Allgatherv(localCode, codeSize[mpiRank], 
413
414
415
416
			       MPI_UNSIGNED_LONG, 
			       code, codeSize, codeOffset,
			       MPI_UNSIGNED_LONG);
    
417
418
    for (int rank = 0; rank < mpiSize; rank++) {
      if (rank != mpiRank) {
419
420
421
422
	std::vector<unsigned long int> remoteCode;
	unsigned long int *ptr;
	unsigned long int *begin = code + codeOffset[rank]; 
	unsigned long int *end = begin + codeSize[rank];
423
	for (ptr = begin; ptr != end; ++ptr) {
424
425
426
427
428
429
430
	  remoteCode.push_back(*ptr);
	}
	structures[rank].init(remoteCode, elements[rank]);
      }
    }

    // free memory
431
    FREE_MEMORY(elements, int, mpiSize);
432
    FREE_MEMORY(code, unsigned long int, codeSizeSum);
433
434
435
    FREE_MEMORY(localCode, unsigned long int, codeSize[mpiRank]);
    FREE_MEMORY(codeOffset, int, mpiSize);
    FREE_MEMORY(codeSize, int, mpiSize);
436
437
438
439
440
441
442
    
  }

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

443
    MeshStructure *structures = NEW MeshStructure[mpiSize];
444

445
    // build composite mesh structure
446
447
448
    exchangeMeshStructureCodes(structures);

    // merge codes
449
450
451
    for (int rank = 0; rank < mpiSize; rank++) {
      if (rank != mpiRank) {
	structures[mpiRank].merge(&structures[rank]);
452
453
      }
    }
454
455
456
457
458
459
460
461
462
463
464
465
466
  
    // build finest mesh on the rank partition
    structures[mpiRank].fitMeshToStructure(mesh,
					   refinementManager,
					   true);

    // In debug mode, process 0 builds the global solution mesh.
    if (debugMode) {
      dbgMesh = mesh;
      if (mpiRank == 0) {
	structures[mpiRank].fitMeshToStructure(dbgMesh,
					       refinementManager,
					       true, true);
467
468
      }
    }
469
470

    DELETE [] structures;
471
472
473
474
475
476
477
478
  }


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

  void ParallelProblem::exchangeRankSolutions(AdaptInfo *adaptInfo,
488
					      Mesh *workMesh,
489
490
491
492
493
494
495
496
497
498
499
500
					      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();
501
    int dim = workMesh->getDim();
502
503
504
505
506
507
    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();

508
509
    std::vector<std::vector<DegreeOfFreedom> > sendOrder(mpiSize);
    std::vector<std::vector<DegreeOfFreedom> > recvOrder(mpiSize);
510
511
512
513
514
515

    elementPartitions_.clear();

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

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

523
524
      if (partitionData) {
	if (partitionData->getLevel() == 0) {
525
526
527
528
529
	  elementPartition = partitionVec_[element->getIndex()];
	}

	PartitionStatus status = partitionData->getPartitionStatus();

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

	    // collect other partitions element belongs to
535
	    for (int i = 0; i < dim + 1; i++) {
536
537
538
539
	      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);
540
541
542
543
544
545
546
	      }	
	    }

	    coarseElement = element;
	  }


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

550
551
	    for (int i = 0; i < numFcts; i++) {
	      if (status == OVERLAP) {
552
553
554
		// send dofs
		sendOrder[elementPartition].push_back(fineDOFs[i]);
	      } 
555
	      if (status == IN) {
556
		// recv dofs
557
558
559
560
561
		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) {
562
563
564
565
566
567
568
569
570
571
572
573
574
 		    recvOrder[*setIt].push_back(fineDOFs[i]);
		  }
		}
	      }
	    }
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }

    // create send and recv buffers and fill send buffers
575
    DOFVector<double> *solution = rankSolutions[mpiRank];
576
577
578
579
580
581

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

582
583
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
	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
609
610
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
611
612
613
614
615
616
617
618
619
	if(sendBufferSize[partition] > 0) {
	  MPI::COMM_WORLD.Isend(sendBuffer[partition],
				sendBufferSize[partition],
				MPI_DOUBLE,
				partition,
				0);
	}
      }
    }    
620
   
621
    // blocking recieves
622
623
624
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (recvBufferSize[partition] > 0) {
625
626
627
628
629
630
631
632
633
634
635
636
637
	  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
638
639
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
640
	std::vector<DegreeOfFreedom>::iterator dofIt = recvOrder[partition].begin();
641
642
	for (int i = 0; i < recvBufferSize[partition]; i++) {
	  (*(rankSolutions[partition]))[*dofIt] = recvBuffer[partition][i];
643
644
645
646
	  ++dofIt;
	}
      }
    }    
647

648
    // free send and recv buffers
649
650
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
651
	if (sendBufferSize[partition] > 0)
652
653
654
655
656
657
658
659
660
	  FREE_MEMORY(sendBuffer[partition], 
		      double,
		      sendBufferSize[partition]);
	if(recvBufferSize[partition] > 0)
	  FREE_MEMORY(recvBuffer[partition], 
		      double,
		      recvBufferSize[partition]);
      }
    }    
661

662
663
664
665
666
667
668
669
670
671
672
673
674
    FREE_MEMORY(coarseDOFs, DegreeOfFreedom, numFcts);
    FREE_MEMORY(fineDOFs, DegreeOfFreedom, numFcts);
  }

  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;
675
676
    sendOrder.resize(mpiSize);
    recvOrder.resize(mpiSize);
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699

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

700
	if(oldPartition == mpiRank) {
701
702
703
704
705
	  for(i = 0; i < numFcts; i++) {
	    // send element values to new partition
	    sendOrder[newPartition].push_back(dofs[i]);
	  }
	}
706
	if(newPartition == mpiRank) {
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
	  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;
726
727
    for(partition = 0; partition < mpiSize; partition++) {
      if(partition != mpiRank) {
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
	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", 
753
// 	mpiRank,
754
755
756
757
758
759
760
761
762
763
// 	sendBufferSize[0],
// 	sendBufferSize[1],
// 	sendBufferSize[2],
// 	sendBufferSize[3],
// 	recvBufferSize[0],
// 	recvBufferSize[1],
// 	recvBufferSize[2],
// 	recvBufferSize[3]);

    // === non-blocking sends ===
764
765
    for(partition = 0; partition < mpiSize; partition++) {
      if(partition != mpiRank) {
766
767
768
769
770
771
772
773
774
775
776
	if(sendBufferSize[partition] > 0) {
	  MPI::COMM_WORLD.Isend(sendBuffer[partition],
				sendBufferSize[partition],
				MPI_DOUBLE,
				partition,
				0);
	}
      }
    }    
    
    // === blocking receives ===
777
778
    for(partition = 0; partition < mpiSize; partition++) {
      if(partition != mpiRank) {
779
780
781
782
783
784
785
786
787
788
789
790
791
792
	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 ===
793
794
    for(partition = 0; partition < mpiSize; partition++) {
      if(partition != mpiRank) {
795
796
797
798
799
800
801
802
803
	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 ===
804
805
    for(partition = 0; partition < mpiSize; partition++) {
      if(partition != mpiRank) {
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
	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();
825
    int dim = mesh->getDim();
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
    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");

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


852
    MyDualTraverse dualTraverse(localCoarseGridLevel_);
853
854
855
    ElInfo *elInfo1, *elInfo2;
    ElInfo *large, *small;

856
857
858
859
860
861
862
863
864
    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);
865
866
867
868
869
870
871

    while (cont) {
      Element *element1 = elInfo1->getElement();
      Element *element2 = elInfo2->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>
	(element1->getElementData(PARTITION_ED));
872
      if (partitionData->getPartitionStatus() == IN) {
873
874

	// get coarse dofs
875
	if (element1 != lastCoarseElement) {
876
877
878
879
	  basFcts->getLocalIndices(element1, admin, coarseDOFs);
	  lastCoarseElement = element1;
	}
      
880
	if (elementPartitions_[element1].size() > 1) {
881
882
883
884
	  // get fine dofs
	  basFcts->getLocalIndices(element2, admin, fineDOFs);
      
	  // for all fine DOFs
885
886
	  for (int i = 0; i < numFcts; i++) {
	    if (!visited[fineDOFs[i]]) {
887
888
889
890
891
892
	      visited[fineDOFs[i]] = true;

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

	      // for all coarse vertex DOFs
893
	      for (int j = 0; j < dim + 1; j++) {
894
895
		partBegin = vertexPartitions_[coarseDOFs[j]].begin();
		partEnd = vertexPartitions_[coarseDOFs[j]].end();
896
		for (partIt = partBegin; partIt != partEnd; ++partIt) {
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
		  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);
    }

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

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

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

920
    for (wIt = wBegin; wIt != wEnd; ++wIt) {
921
922
923
924
      DegreeOfFreedom dof = wIt->first;
      (*globalSolution)[dof] = 0.0;
    }
    
925
    for (wIt = wBegin; wIt != wEnd; ++wIt) {
926
927
928
929
930
      DegreeOfFreedom dof = wIt->first;
      std::map<int, double>::iterator partIt, partBegin, partEnd;
      partBegin = wIt->second.begin();
      partEnd = wIt->second.end();
    
931
      for (partIt = partBegin; partIt != partEnd; ++partIt) {
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
	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;
949
    ElInfo *elInfo = stack.traverseFirst(mesh, 
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
					 -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;
985
    ElInfo *elInfo = stack.traverseFirst(mesh, 
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
					 -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;
1019
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
    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);
    }

1032
    refinementManager->refineMesh(mesh);
1033
1034
1035
1036
1037
1038
  }


  void ParallelProblem::createOverlap(int level, int overlap, bool openOverlap,
				      std::map<Element*, int> &overlapDistance)
  {
1039
    int i, dim = mesh->getDim();
1040
1041
1042
1043
1044
1045

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

    TraverseStack stack;
1046
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
    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)
  {
1109
    int dim = mesh->getDim();
1110
1111
1112
1113
1114
1115
1116
1117
1118

    TraverseStack stack;
    ElInfo *elInfo;

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

    // first: partition elements ...
    int index, partition;
1119
    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
1120
    while (elInfo) {
1121
1122
1123
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
1124
1125
      if (partitionData) {
	if (partitionData->getLevel() == 0) {
1126
1127
1128
	  index = element->getIndex();
	  partition = partitionVec_[index];
	}
1129
1130
1131
	if (partitionData->getLevel() == level) {
	  for (int i = 0; i < dim + 1; i++) {
	    vertexPartitions_[element->getDOF(i, 0)].insert(partition);
1132
1133
1134
1135
1136
1137
	  }
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }

1138
    if (overlap > 1 || openOverlap == false) {
1139
      // exchange mesh structure codes
1140
      MeshStructure *structures = NEW MeshStructure[mpiSize];
1141
1142
1143
1144
      exchangeMeshStructureCodes(structures);

      // merge codes
      int rank;
1145
1146
1147
      for(rank = 0; rank < mpiSize; rank++) {
	if(rank != mpiRank) {
	  structures[mpiRank].merge(&structures[rank]);
1148
1149
1150
	}
      }
    
1151
      MeshStructure &compositeStructure = structures[mpiRank];
1152
1153
1154
1155
1156
1157
      compositeStructure.reset();

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

1158
1159
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
      while (elInfo) {
1160
1161
1162
1163
	Element *element = elInfo->getElement();
	PartitionElementData *partitionData = 
	  dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
	
1164
	if (partitionData && partitionData->getLevel() == level) {
1165
1166
	  int compositeIndex = compositeStructure.getCurrentElement();
	  indexElement[compositeIndex] = element;
1167
	  if (partitionData->getPartitionStatus() == OVERLAP) {
1168
	    int distance = overlapDistance[element];
1169
	    if (distance < overlap || !openOverlap) {
1170
1171
1172
1173
1174
	      innerOverlapElements.push_back(compositeIndex);
	    } 
	  }
	}
	
1175
	if (element->isLeaf()) {
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
	  compositeStructure.skipBranch();
	} else {
	  compositeStructure.nextElement();
	}
	elInfo = stack.traverseNext(elInfo);
      }
    
      // === exchange 'inner' overlap elements ===

      // exchange number of overlap elements
1186
      int *numOverlapElements = GET_MEMORY(int, mpiSize);
1187
1188
1189
1190
1191
      int tmp = static_cast<int>(innerOverlapElements.size());
      MPI::COMM_WORLD.Allgather(&tmp, 1, MPI_INT,
			      numOverlapElements, 1, MPI_INT);
      
      // exchange overlap elements
1192
      int *offset = GET_MEMORY(int, mpiSize);
1193
      int sum = 0;
1194
      for(rank = 0; rank < mpiSize; rank++) {
1195
1196
1197
1198
1199
	offset[rank] = sum;
	sum += numOverlapElements[rank];
      }

      int *recvBuffer = GET_MEMORY(int, sum);
1200
      int *sendBuffer = GET_MEMORY(int, numOverlapElements[mpiRank]);
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
   
      int *ptr;
      std::vector<int>::iterator elemIt, elemEnd = innerOverlapElements.end();
  
      for(ptr = sendBuffer, elemIt = innerOverlapElements.begin();
	  elemIt != elemEnd; 
	  ++elemIt, ++ptr) 
      {
	*ptr = *elemIt;
      }
  
1212
      MPI::COMM_WORLD.Allgatherv(sendBuffer, numOverlapElements[mpiRank], 
1213
1214
1215
1216
1217
1218
1219
1220
				 MPI_INT, 
				 recvBuffer, numOverlapElements, offset,
				 MPI_INT);



      // fill vertexPartitions for 'inner' overlap elements
      int el;
1221
      for (rank = 0; rank < mpiSize; rank++) {
1222
1223
	int numElements = numOverlapElements[rank];
	//       MSG("rank %d overlap elements with %d: %d\n", 
1224
	// 	  mpiRank, rank, numOverlapElements[rank]);
1225
	int *elements = recvBuffer + offset[rank];
1226
	for (el = 0; el < numElements; el++) {
1227
	  Element *element = indexElement[elements[el]];
1228
	  for (int i = 0; i < dim + 1; i++) {
1229
1230
1231
1232
1233
1234
1235
1236
	    vertexPartitions_[element->getDOF(i, 0)].insert(rank/* + 1*/);
	  }
	}
      }

      // free memory
      DELETE [] structures;
      FREE_MEMORY(recvBuffer, int, sum);
1237
1238
      FREE_MEMORY(sendBuffer, int, numOverlapElements[mpiRank]);
      FREE_MEMORY(numOverlapElements, int, mpiSize);
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
    }
  }

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

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

  // =========================================================================
  // ===== class ParallelProblemScal =========================================
  // =========================================================================
  
  ParallelProblemScal::ParallelProblemScal(const std::string& name,
1255
					   ProblemScal *prob,
1256
1257
1258
					   ProblemInstatScal *problemInstat,
					   std::vector<DOFVector<double>*> vectors)
    : ParallelProblem(name,
1259
		      prob,
1260
1261
		      problemInstat,
		      vectors,
1262
1263
1264
1265
		      prob->getMesh(),
		      prob->getRefinementManager(),
		      prob->getCoarseningManager()),
      problem(prob),
1266
      problemInstat_(problemInstat),
1267
1268
1269
1270
1271
      oldEstimator(NULL),
      oldMarker(NULL),
      rankSolution(mpiSize),
      usersEstimator(NULL),
      usersMarker(NULL)
1272
  {
1273
1274
    for (int i = 0; i < mpiSize; i++) {
      rankSolution[i] = NEW DOFVector<double>(problem->getFESpace(), "rank solution");
1275
1276
    }

1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
    // Create vectors for debuging information
    if (debugMode) {
      dbgRankSolution.resize(mpiSize);
      for (int i = 0; i < mpiSize; i++) {
	dbgRankSolution[i] = NEW DOFVector<double>(problem->getFESpace(), "debug rank solution");
      }
      
      if (mpiRank == 0) {
	dbgSolution = NEW DOFVector<double>(problem->getFESpace(), "debug solution");
      }
    }
1288

1289
1290
    
    if (problemInstat_) {
1291
1292
1293
1294
1295
1296
      dofVectors_.push_back(problemInstat_->getOldSolution());
    }
  }

  ParallelProblemScal::~ParallelProblemScal()
  {
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
    for (int i = 0; i < mpiSize; i++) {
      DELETE rankSolution[i];
    }

    if (debugMode) {
      for (int i = 0; i < mpiSize; i++) {
	DELETE dbgRankSolution[i];
      }

      if (mpiRank == 0) {
	DELETE dbgSolution;
      }
1309
1310
1311
1312
1313
1314
    }
  }


  void ParallelProblemScal::initParallelization(AdaptInfo *adaptInfo)
  {
1315
    if (mpiSize > 1) {
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
      clock_t partitioningStart = clock();

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

      globalRefineOutOfPartition(adaptInfo);

      refineOverlap(adaptInfo);
      createOverlap(adaptInfo);

      globalRefinements();

1329
1330
      oldEstimator = problem->getEstimator();
      oldMarker = problem->getMarker();
1331
1332

      ConditionalEstimator *condEstimator = NULL;
1333
1334
      if (usersEstimator) {
	problem->setEstimator(usersEstimator);
1335
      } else {
1336
1337
	condEstimator = NEW ConditionalEstimator(oldEstimator);
	problem->setEstimator(condEstimator);
1338
1339
      }

1340
1341
      if (usersMarker) {
	problem->setMarker(usersMarker);
1342
1343
1344
      } else {
	ConditionalMarker *newMarker = NEW ConditionalMarker("parallel marker",
							     -1,
1345
							     oldMarker,
1346
1347
							     globalCoarseGridLevel_,
							     localCoarseGridLevel_);
1348
	problem->setMarker(newMarker);
1349
1350
      }

1351
      if (mpiRank == 0) {
1352
1353
1354
1355
1356
1357
1358
1359
	clock_t partitioningEnd = clock();
	partitioningTime = TIME_USED(partitioningStart, 
				     partitioningEnd);
	computationStart = partitioningEnd;
      }

      // modify file writers
      char number[10];
1360
1361
1362
      sprintf(number, "%d", mpiRank);
      ::std::vector<FileWriterInterface*> fileWriters = problem->getFileWriterList();
      ::std::vector<FileWriterInterface*>::iterator fwIt, fwBegin, fwEnd;
1363
1364
      fwBegin = fileWriters.begin();
      fwEnd = fileWriters.end();
1365
      for (fwIt = fwBegin; fwIt != fwEnd; ++fwIt) {
1366
1367
1368
1369
1370
1371
1372
1373
1374
	(*fwIt)->setFilename((*fwIt)->getFilename() + "_proc" + 
			     ::std::string(number) + "_");
	(*fwIt)->setTraverseProperties(-1, 0, elementInPartition);
      }
    }
  }

  void ParallelProblemScal::exitParallelization(AdaptInfo *adaptInfo)
  {
1375
    if (mpiSize > 1) {
1376
1377
      ParallelProblem::exitParallelization(adaptInfo);

1378
1379
      if (!timeIF_) 
	problem->writeFiles(adaptInfo, true);
1380
1381
1382

      partitioner_->deletePartitionData();

1383
1384
1385
1386
      if (!usersEstimator) 
	DELETE problem->getEstimator();
      if (!usersMarker) 
	DELETE problem->getMarker();
1387

1388
1389
      problem->setEstimator(oldEstimator);
      problem->setMarker(oldMarker);    
1390
1391
1392
1393
1394
1395
1396
    }
  }

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

  ProblemStatBase *ParallelProblemScal::getProblem(int number) 
  { 
1397
    return problem; 
1398
1399
1400
1401
  }

  ProblemStatBase *ParallelProblemScal::getProblem(const std::string& name)
  { 
1402
    return problem;
1403
1404
1405
1406
1407
1408
1409
1410
  }


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

  void ParallelProblemScal::exchangeRankSolutions(AdaptInfo *adaptInfo)
  {
1418
    rankSolution[mpiRank]->copy(*(problem->getSolution()));
1419
    ParallelProblem::exchangeRankSolutions(adaptInfo,
1420
1421
1422
1423
1424
1425
1426
1427
1428
					   mesh,
					   rankSolution);

    if (debugMode) {
      dbgRankSolution[mpiRank]->copy(*(problem->getSolution()));
      ParallelProblem::exchangeRankSolutions(adaptInfo,
					     dbgMesh,
					     dbgRankSolution);
    }
1429
1430
1431
1432
1433
  }

  void ParallelProblemScal::buildGlobalSolution(AdaptInfo *adaptInfo)
  {
    ParallelProblem::buildGlobalSolution(adaptInfo,
1434
1435
					 rankSolution,
					 problem->getSolution());
1436

1437
1438
1439
1440
1441
    if (debugMode && mpiRank == 0) {
      ParallelProblem::buildGlobalSolution(adaptInfo,
					   dbgRankSolution,
					   dbgSolution);
    }
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
  }

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


  // =========================================================================
  // ===== class ParallelProblemVec ==========================================
  // =========================================================================
  
  ParallelProblemVec::ParallelProblemVec(const std::string& name,
1458
					 ProblemVec *prob,
1459
1460
1461
					 ProblemInstatVec *problemInstat,
					 std::vector<DOFVector<double>*> vectors)
    : ParallelProblem(name,
1462
		      prob,
1463
1464
		      problemInstat,
		      vectors,
1465
1466
1467
1468
		      prob->getMesh(0),
		      prob->getRefinementManager(0),
		      prob->getCoarseningManager(0)),
      problem(prob),
1469
1470
      problemInstat_(problemInstat)
  {
1471
    nComponents = problem->getNumComponents();
1472

1473
1474
1475
    std::vector<FiniteElemSpace*> feSpaces(nComponents);
    for (int i = 0; i < nComponents; i++) {
      feSpaces[i] = problem->getFESpace(i);
1476
    }
1477
1478
1479
1480
1481
1482

    rankSolution.resize(mpiSize);
    for (int i = 0; i < mpiSize; i++) {
      rankSolution[i] = NEW SystemVector("rank solution", feSpaces, nComponents);
      for (int j = 0; j < nComponents; j++) {
	rankSolution[i]->setDOFVector(j, NEW DOFVector<double>(feSpaces[j], 
1483
1484
1485
								"rank solution"));
      }
    }
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499

    if (debugMode) {
      dbgRankSolution.resize(mpiSize);
      for (int i = 0; i < mpiSize; i++) {
	dbgRankSolution[i] = NEW SystemVector(*(rankSolution[i]));
      }

      if (mpiRank == 0) {
	dbgSolution = NEW SystemVector(*(problem->getSolution()));
      }
    }

    if (problemInstat_) {
      for (int i = 0; i < nComponents; i++) {
1500
1501
1502
1503
1504
1505
1506
	dofVectors_.push_back(problemInstat_->getOldSolution()->getDOFVector(i));
      }
    }
  }

  ParallelProblemVec::~ParallelProblemVec()
  {
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
    for (int i = 0; i < mpiSize; i++) {
      for (int j = 0; j < nComponents; j++) {
	DELETE rankSolution[i]->getDOFVector(j);
      }
      DELETE rankSolution[i];
    }

    if (debugMode) {
      for (int i = 0; i < mpiSize; i++) {
	for (int j = 0; j < nComponents; j++) {
	  DELETE dbgRankSolution[i]->getDOFVector(j);
	}
	DELETE dbgRankSolution[i];
      }

      if (mpiRank == 0) {
	for (int i = 0; i < nComponents; i++) {
	  DELETE dbgSolution->getDOFVector(i);
	}
	DELETE dbgSolution;
1527
1528
1529
1530
1531
1532
1533
1534
      }
    }
  }

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

1535
    if (mpiSize > 1) {
1536
1537
1538
1539
1540
1541
1542
      partitioner_->createPartitionData();
      setElemWeights(adaptInfo);
      partitionMesh(adaptInfo);

      refineOverlap(adaptInfo);
      createOverlap(adaptInfo);

1543
1544
      oldEstimator = problem->getEstimator();
      oldMarker = problem->getMarker();
1545
1546

      std::vector<ConditionalEstimator*> condEstimator;
Thomas Witkowski's avatar
*