ParallelProblem.cc 55.5 KB
Newer Older
1
2
3
#include <queue>
#include <time.h>
#include "boost/lexical_cast.hpp"
4
#include "ParallelProblem.h"
5
6
7
8
9
10
11
12
13
14
#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"
15
#include "MacroElement.h"
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#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"
30
#include "VtkWriter.h"
31
32
33
34
#include "mpi.h"

namespace AMDiS {

35
36
  using boost::lexical_cast;

37
38
  bool elementInPartition(ElInfo *elInfo)
  {
39
40
41
42
43
44
45
    PartitionElementData *elementData = dynamic_cast<PartitionElementData*>
      (elInfo->getElement()->getElementData(PARTITION_ED));
    if (elementData && elementData->getPartitionStatus() == IN) {
      return true;
    } else {
      return false;
    }
46
47
48
49
50
51
52
  }

  class MyDualTraverse : public DualTraverse
  {
  public:
    MyDualTraverse(int coarseLevel)
      : coarseLevel_(coarseLevel)
Thomas Witkowski's avatar
Thomas Witkowski committed
53
    {}
54
55
56
57
58

    bool skipEl1(ElInfo *elInfo)
    {
      PartitionElementData *elementData = dynamic_cast<PartitionElementData*>
	(elInfo->getElement()->getElementData(PARTITION_ED));
59
60
      if (elementData) {
	if (elInfo->getElement()->isLeaf() && 
61
62
63
64
65
66
	   elementData->getLevel() < coarseLevel_)
	  return false;
	if(elementData->getLevel() == coarseLevel_)
	  return false;
      }
      return true;
Thomas Witkowski's avatar
Thomas Witkowski committed
67
    }
Naumann, Andreas's avatar
Naumann, Andreas committed
68
69
70
71
    
    /*bool skipEl2(ElInfo *elInfo)
    {
	    return skipEl1(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
72
    }*/
73
74
75
76
77
78
79
80
  private:
    int coarseLevel_;
  };

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

Thomas Witkowski's avatar
Thomas Witkowski committed
81
  ParallelProblemBase::ParallelProblemBase(std::string name,
82
					   ProblemIterationInterface *iterationIF,
83
84
					   ProblemTimeInterface *timeIF)
    : iterationIF_(iterationIF),
85
      timeIF_(timeIF),
86
87
      debugMode(0),
      debugServerProcess(false)
88
  {
89
90
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
91
92
93

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

Naumann, Andreas's avatar
Naumann, Andreas committed
94
    groupIsSet=false;
95
96
97
98
    if (debugMode) {
      MPI::Group group = MPI::COMM_WORLD.Get_group();
      
      int rankSize = mpiSize - 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
99
100
      int *ranks = new int[rankSize];
      for (int i = 0; i < rankSize; i++)
101
	ranks[i] = i + 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
102

Naumann, Andreas's avatar
Naumann, Andreas committed
103
      std::cout<<"set amdis-group"<<std::endl;
104
      amdisGroup = group.Incl(rankSize, ranks);
Naumann, Andreas's avatar
Naumann, Andreas committed
105
      groupIsSet=true;
106
107
108
109
110
111
112
113
114
115
      mpiComm = MPI::COMM_WORLD.Create(amdisGroup);
      
      if (mpiComm != MPI::COMM_NULL) {
	mpiRank = mpiComm.Get_rank();
	mpiSize = mpiComm.Get_size();
	debugServerProcess = false;
      } else {
	debugServerProcess = true;
      }
      
Thomas Witkowski's avatar
Thomas Witkowski committed
116
      delete [] ranks;
117
118
119
    } else {
      mpiComm = MPI::COMM_WORLD;
    }
120
121
  }

122
123
124
125
  void ParallelProblemBase::exitParallelization(AdaptInfo *adaptInfo)
  {
    if (!timeIF_) 
      closeTimestep(adaptInfo);
Naumann, Andreas's avatar
Naumann, Andreas committed
126
127
128
129
130
131
132
133
#ifdef DEBUG
    std::cout<<"free amdisGroup"<<std::endl;
#endif
    if(groupIsSet)
    	amdisGroup.Free();
#ifdef DEBUG
    std::cout<<"freeD amdisGroup"<<std::endl;
#endif    
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
  }

  void ParallelProblemBase::closeTimestep(AdaptInfo *adaptInfo)
  {
    if (mpiSize > 1 && doBuildGlobalSolution(adaptInfo)) {
      if (debugMode && mpiRank == 0) {
	// Send adaptInfo inforamtion to debug server
	double sendTime = adaptInfo->getTime();
	double sendTimestep = adaptInfo->getTimestep();
	MPI::COMM_WORLD.Send(&sendTime, 1, MPI_DOUBLE, 0, 100);
	MPI::COMM_WORLD.Send(&sendTimestep, 1, MPI_DOUBLE, 0, 100);
      }
      synchronizeMeshes(adaptInfo);	
      exchangeRankSolutions(adaptInfo);
      buildGlobalSolution(adaptInfo);
    }
    
    if (timeIF_) 
      timeIF_->closeTimestep(adaptInfo);
  }
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185

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

	synchronizeMeshes(adaptInfo);
	partitionMesh(adaptInfo);
	refineOverlap(adaptInfo);
	createOverlap(adaptInfo);
	synchronizeMeshes(adaptInfo);
	exchangeDOFVectors(adaptInfo);
	coarsenOutOfPartition(adaptInfo);
	
	clock_t partitioningEnd = clock();
	partitioningTime = TIME_USED(partitioningStart, 
				     partitioningEnd);
	computationStart = partitioningEnd;
      }
      
      flag |= iterationIF_->oneIteration(adaptInfo, toDo & ~(MARK | ADAPT));
    } else {
      flag = iterationIF_->oneIteration(adaptInfo, toDo);
    }
    
    // synchronize adaption flag
Thomas Witkowski's avatar
Thomas Witkowski committed
186
    unsigned long *flagBuffer = new unsigned long[mpiSize];
187
188
189
190
191
    
    unsigned long localFlag = flag.getFlags();
    
    mpiComm.Allgather(&localFlag, 1, MPI_UNSIGNED_LONG,
		      flagBuffer, 1, MPI_UNSIGNED_LONG);
Thomas Witkowski's avatar
Thomas Witkowski committed
192
    for (int i = 0; i < mpiSize; i++)
193
      flag.setFlag(flagBuffer[i]);
Thomas Witkowski's avatar
Thomas Witkowski committed
194
195
    
    delete [] flagBuffer;
196
197
    
    return flag;
Thomas Witkowski's avatar
Thomas Witkowski committed
198
  }
199

200
201
202
203
  // =========================================================================
  // ===== class ParallelProblem =============================================
  // =========================================================================

Thomas Witkowski's avatar
Thomas Witkowski committed
204
  ParallelProblem::ParallelProblem(std::string name,
205
206
207
				   ProblemIterationInterface *iterationIF,
				   ProblemTimeInterface *timeIF,
				   std::vector<DOFVector<double>*> vectors,
208
				   Mesh *mesh_,
209
210
				   RefinementManager *rm,
				   CoarseningManager *cm)
211
    : ParallelProblemBase(name, iterationIF, timeIF),
212
      name_(name),
213
      mesh(mesh_),
214
215
      refinementManager(rm),
      coarseningManager(cm),
216
217
218
219
220
221
222
223
224
225
226
      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),
227
      repartTimeFactor_(10.0)
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
  {
    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");

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268

    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() 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    delete partitioner;
270
271
272
273
274
275
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
    double *weightSum = new double[mpiSize];
    int *partArray = new int[mpiSize];
278
    int part = 0;
279

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

282
    if (mpiRank == 0) {
283
284

      double average = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
285
      for (int i = 0; i < mpiSize; i++)
286
287
	average += weightSum[i];

288
      average /= mpiSize;
289

290
291
      for (int i = 0; i < mpiSize; i++) {
	if ((weightSum[i] / average) > upperPartThreshold_) {
292
293
294
	  part = 1;
	  break;
	}
295
	if ((weightSum[i] / average) < lowerPartThreshold_) {
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
	  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_);
	}
      }

325
      for (int i = 0; i < mpiSize; i++) {
326
327
328
329
	partArray[i] = part;
      }      
    }

330
331
    mpiComm.Scatter(partArray, 1, MPI_INT,
		    &part, 1, MPI_INT, 0);
332
    
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
    delete [] weightSum;
    delete [] partArray;
335
336
337
338

    return (part == 1);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
  bool ParallelProblem::doBuildGlobalSolution(AdaptInfo *adaptInfo) 
  {
341
342
343
344
345
346
    return true;
  }

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

348
    if (initial) {
349
      initial = false;
350
      partitioner->fillCoarsePartitionVec(&oldPartitionVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
351
      partitioner->partition(&elemWeights, INITIAL);
352
    } else {
353
      oldPartitionVec = partitionVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
354
      partitioner->partition(&elemWeights, ADAPTIVE_REPART, 100.0 /*0.000001*/);
355
356
    }    

357
    partitioner->fillCoarsePartitionVec(&partitionVec);
358
359
360
361
  }

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

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

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

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

379
380
	if (partitionData->getPartitionStatus() == IN) {
	  for (int i = 0; i < dim + 1; i++) {
381
	    DegreeOfFreedom dof = dofs[i][0];
382
383
384
385
	    if (inOut[dof] == 2) 
	      inOut[dof] = 3;
	    if (inOut[dof] == 0) 
	      inOut[dof] = 1;
386
387
	  }
	} else {
388
	  for (int i = 0; i < dim + 1; i++) {
389
	    DegreeOfFreedom dof = dofs[i][0];
390
391
392
393
	    if (inOut[dof] == 1) 
	      inOut[dof] = 3;
	    if (inOut[dof] == 0) 
	      inOut[dof] = 2;
394
395
396
397
398
399
400
401
402
	  }
	}

	elInfo = stack.traverseNext(elInfo);
      }

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

	int level = partitionData->getLevel();

411
412
	if (level < localCoarseGridLevel_) {
	  if (partitionData->getPartitionStatus() != IN) {
413
	    const DegreeOfFreedom **dofs = element->getDOF();
414
	    for (int i = 0; i < dim + 1; i++) {
415
	      DegreeOfFreedom dof = dofs[i][0];
416
	      if (inOut[dof] == 3) {
417
418
		element->setMark(1);
		marked = true;
419
420
		if ((level + 1) < localCoarseGridLevel_) 
		  finished = false;
421
422
423
424
425
426
		break;
	      }
	    }
	  } else {
	    element->setMark(1);
	    marked = true;
427
428
	    if ((level + 1) < localCoarseGridLevel_) 
	      finished = false;
429
430
431
432
433
	  }
	}

	elInfo = stack.traverseNext(elInfo);
      }
434
435
      if (marked) 
	refinementManager->refineMesh(mesh);
436
437
438
439
440
441
    }
  }

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

452
    refinementManager->refineMesh(mesh);
453
454
455
456
457
  }

  void ParallelProblem::coarsenOutOfPartition(AdaptInfo *adaptInfo)
  {
    Flag meshCoarsened = 1;    
458
    while (meshCoarsened.getFlags() != 0) {
459
      TraverseStack stack;
460
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
461
      while (elInfo) {
462
463
464
465
	Element *element = elInfo->getElement();
	PartitionElementData *partitionData = 
	  dynamic_cast<PartitionElementData*>
	  (element->getElementData(PARTITION_ED));
466
	if (partitionData->getPartitionStatus() == OUT) {
467
468
469
470
471
	  int mark = min(0, -partitionData->getLevel() + globalCoarseGridLevel_);
	  element->setMark(mark);
	}
	elInfo = stack.traverseNext(elInfo);
      }
472
      meshCoarsened = coarseningManager->coarsenMesh(mesh);
473
    }
474
    mpiComm.Barrier();
475
476
477
478
  }

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

    // broadcast code sizes
Thomas Witkowski's avatar
Thomas Witkowski committed
484
    int *codeSize = new int[mpiSize];
485
    int tmp = static_cast<int>(myCode.size());
486
487
488
489
490
    mpiComm.Allgather(&tmp, 1, MPI_INT, codeSize, 1, MPI_INT);
    if (debugMode) {
      // send code sizes also to debug server
      MPI::COMM_WORLD.Gather(&tmp, 1, MPI_INT, NULL, 1, MPI_INT, 0);
    }
491
492

    // broadcast number of elements
Thomas Witkowski's avatar
Thomas Witkowski committed
493
    int *elements = new int[mpiSize];
494
    tmp = structures[mpiRank].getNumElements();
495
496
497
498
499
    mpiComm.Allgather(&tmp, 1, MPI_INT, elements, 1, MPI_INT);
    if (debugMode) {
      // send number of elements also to debug server
      MPI::COMM_WORLD.Gather(&tmp, 1, MPI_INT, NULL, 1, MPI_INT, 0);
    }
500
501

    // broadcast codes
Thomas Witkowski's avatar
Thomas Witkowski committed
502
    int *codeOffset = new int[mpiSize];
503
    int codeSizeSum = 0;
504
    for (int rank = 0; rank < mpiSize; rank++) {
505
506
507
508
      codeOffset[rank] = codeSizeSum;
      codeSizeSum += codeSize[rank];
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
509
510
    unsigned long int *code = new unsigned long int[codeSizeSum];
    unsigned long int *localCode = new unsigned long int[codeSize[mpiRank]];
511
512
513
    unsigned long int *ptr;
    std::vector<unsigned long int>::const_iterator it, end = myCode.end();
  
Thomas Witkowski's avatar
Thomas Witkowski committed
514
    for (ptr = localCode, it = myCode.begin(); it != end; ++it, ++ptr)
515
      *ptr = *it;
516
517
518
519
520
521
522
523
524
525

    mpiComm.Allgatherv(localCode, codeSize[mpiRank], 
		       MPI_UNSIGNED_LONG, 
		       code, codeSize, codeOffset,
		       MPI_UNSIGNED_LONG);
    if (debugMode) {
      // send codes also to debug server
      MPI::COMM_WORLD.Send(localCode, codeSize[mpiRank],
			   MPI_UNSIGNED_LONG, 0, 100);
    }
526
    
527
528
    for (int rank = 0; rank < mpiSize; rank++) {
      if (rank != mpiRank) {
529
530
531
532
	std::vector<unsigned long int> remoteCode;
	unsigned long int *ptr;
	unsigned long int *begin = code + codeOffset[rank]; 
	unsigned long int *end = begin + codeSize[rank];
533
	for (ptr = begin; ptr != end; ++ptr) {
534
535
536
537
538
539
	  remoteCode.push_back(*ptr);
	}
	structures[rank].init(remoteCode, elements[rank]);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
540
541
542
543
544
    delete [] elements;
    delete [] code;
    delete [] localCode;
    delete [] codeOffset;
    delete [] codeSize;    
545
546
547
548
549
550
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
551
    MeshStructure *structures = new MeshStructure[mpiSize];
552

553
    // build composite mesh structure
554
555
556
    exchangeMeshStructureCodes(structures);

    // merge codes
Thomas Witkowski's avatar
Thomas Witkowski committed
557
558
    for (int rank = 0; rank < mpiSize; rank++)
      if (rank != mpiRank)
559
560
561
562
563
564
565
	structures[mpiRank].merge(&structures[rank]);
  
    // build finest mesh on the rank partition
    structures[mpiRank].fitMeshToStructure(mesh,
					   refinementManager,
					   true);

Thomas Witkowski's avatar
Thomas Witkowski committed
566
    delete [] structures;
567
568
569
570
571
572
573
574
  }


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

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

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

592
    const FiniteElemSpace *feSpace = rankSolutions[0]->getFeSpace();
593
    int dim = workMesh->getDim();
594
595
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int numFcts = basFcts->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
596
597
    DegreeOfFreedom *coarseDOFs = new DegreeOfFreedom[numFcts];
    DegreeOfFreedom *fineDOFs = new DegreeOfFreedom[numFcts];
598
    DOFAdmin *admin = feSpace->getAdmin();
599

600
601
    std::vector<std::vector<DegreeOfFreedom> > sendOrder(mpiSize);
    std::vector<std::vector<DegreeOfFreedom> > recvOrder(mpiSize);
602
603
604
605
606
607

    elementPartitions_.clear();

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

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

615
616
      if (partitionData) {
	if (partitionData->getLevel() == 0) {
617
	  elementPartition = partitionVec[element->getIndex()];
618
619
620
621
	}

	PartitionStatus status = partitionData->getPartitionStatus();

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

	    // collect other partitions element belongs to
627
	    for (int i = 0; i < dim + 1; i++) {
628
629
	      std::set<int>::iterator setBegin = vertexPartitions[coarseDOFs[i]].begin();
	      std::set<int>::iterator setEnd = vertexPartitions[coarseDOFs[i]].end();
630
631
	      for (std::set<int>::iterator setIt = setBegin; setIt != setEnd; ++setIt) {
		elementPartitions_[element].insert(*setIt);
632
633
634
635
636
637
638
	      }	
	    }

	    coarseElement = element;
	  }


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

642
643
	    for (int i = 0; i < numFcts; i++) {
	      if (status == OVERLAP) {
644
645
646
		// send dofs
		sendOrder[elementPartition].push_back(fineDOFs[i]);
	      } 
647
	      if (status == IN) {
648
		// recv dofs
649
650
651
652
653
		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) {
654
655
656
657
658
659
660
661
662
663
664
665
666
 		    recvOrder[*setIt].push_back(fineDOFs[i]);
		  }
		}
	      }
	    }
	  }
	}
      }
      
      elInfo = stack.traverseNext(elInfo);
    }

    // create send and recv buffers and fill send buffers
667
    DOFVector<double> *solution = rankSolutions[mpiRank];
668
669
670
671
672
673

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

674
675
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
676
677
678
679
680
	int sendSize = static_cast<int>(sendOrder[partition].size());
	int recvSize = static_cast<int>(recvOrder[partition].size());

	sendBufferSize[partition] = sendSize;	
	recvBufferSize[partition] = recvSize;
681
	if (sendSize > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
682
	  sendBuffer[partition] = new double[sendSize];
683
684
685
686
687
	  std::vector<DegreeOfFreedom>::iterator dofIt;
	  dofIt = sendOrder[partition].begin();
	  double *bufferIt, *bufferBegin, *bufferEnd;
	  bufferBegin = sendBuffer[partition];
	  bufferEnd = bufferBegin + sendSize;
Thomas Witkowski's avatar
Thomas Witkowski committed
688
	  for (bufferIt = bufferBegin; bufferIt < bufferEnd; ++bufferIt, ++dofIt)
689
	    *bufferIt = (*solution)[*dofIt];
690
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
691
692
	if (recvSize > 0)
	  recvBuffer[partition] = new double[recvSize];
693
694
695
696
      }
    }

    // non-blocking sends
697
698
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
699
700
701
702
703
704
	if (sendBufferSize[partition] > 0) {
	  mpiComm.Isend(sendBuffer[partition],
			sendBufferSize[partition],
			MPI_DOUBLE,
			partition,
			0);
705
706
707
	}
      }
    }    
708
   
709
    // blocking recieves
710
711
712
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (recvBufferSize[partition] > 0) {
713
714
715
716
717
	  mpiComm.Recv(recvBuffer[partition],
		       recvBufferSize[partition],
		       MPI_DOUBLE,
		       partition,
		       0);
718
719
720
721
722
	}
      }
    }    

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

    // copy values into rank solutions
726
727
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
728
	std::vector<DegreeOfFreedom>::iterator dofIt = recvOrder[partition].begin();
729
730
	for (int i = 0; i < recvBufferSize[partition]; i++) {
	  (*(rankSolutions[partition]))[*dofIt] = recvBuffer[partition][i];
731
732
733
734
	  ++dofIt;
	}
      }
    }    
735

736
    // free send and recv buffers
737
738
    for (int partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
741
742
	if (sendBufferSize[partition] > 0)
	  delete [] sendBuffer[partition];
	if (recvBufferSize[partition] > 0)
	  delete [] recvBuffer[partition];
743
744
      }
    }    
745

Thomas Witkowski's avatar
Thomas Witkowski committed
746
747
    delete [] coarseDOFs;
    delete [] fineDOFs;
748
749
750
751
752
  }

  void ParallelProblem::exchangeDOFVector(AdaptInfo *adaptInfo,
					  DOFVector<double> *values)
  {
753
754
    partitioner->fillLeafPartitionVec(&oldPartitionVec, &oldPartitionVec);
    partitioner->fillLeafPartitionVec(&partitionVec, &partitionVec);
755
756
757
758

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

762
    const FiniteElemSpace *feSpace = values->getFeSpace();
763
764
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int numFcts = basFcts->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
765
    DegreeOfFreedom *dofs = new DegreeOfFreedom[numFcts];
766
767
768
769
770
771
772
773
774
775
    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();
776
777
      int oldPartition = oldPartitionVec[index];
      int newPartition = partitionVec[index];
778

Thomas Witkowski's avatar
Thomas Witkowski committed
779
      if (oldPartition != newPartition) {
780
781
782
	// get dof indices
	basFcts->getLocalIndices(element, admin, dofs);

Thomas Witkowski's avatar
Thomas Witkowski committed
783
784
785
	if (oldPartition == mpiRank) {
	  // send element values to new partition
	  for (int i = 0; i < numFcts; i++)
786
787
	    sendOrder[newPartition].push_back(dofs[i]);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
788
789
790
	if (newPartition == mpiRank) {
	  // recv element values from old partition
	  for (int i = 0; i < numFcts; i++)
791
792
793
794
795
796
797
	    recvOrder[oldPartition].push_back(dofs[i]);
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
798
    delete [] dofs;
799
800
801
802
803
804
805
806

    // === 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;
807
808
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
809
810
811
812
813
	int sendSize = static_cast<int>(sendOrder[partition].size());
	int recvSize = static_cast<int>(recvOrder[partition].size());
      
	sendBufferSize[partition] = sendSize;	
	recvBufferSize[partition] = recvSize;
814
	if (sendSize > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
815
	  sendBuffer[partition] = new double[sendSize];
816
817
818
819
820
	  std::vector<DegreeOfFreedom>::iterator dofIt;
	  dofIt = sendOrder[partition].begin();
	  double *bufferIt, *bufferBegin, *bufferEnd;
	  bufferBegin = sendBuffer[partition];
	  bufferEnd = bufferBegin + sendSize;
Thomas Witkowski's avatar
Thomas Witkowski committed
821
	  for (bufferIt = bufferBegin; bufferIt < bufferEnd; ++bufferIt, ++dofIt)
822
823
	    *bufferIt = (*values)[*dofIt];
	}
824
	if (recvSize > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
825
	  recvBuffer[partition] = new double[recvSize];
826
827
828
829
      }
    }

    // === non-blocking sends ===
830
831
832
833
834
835
836
837
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (sendBufferSize[partition] > 0) {
	  mpiComm.Isend(sendBuffer[partition],
			  sendBufferSize[partition],
			  MPI_DOUBLE,
			  partition,
			  0);
838
839
840
841
842
	}
      }
    }    
    
    // === blocking receives ===
843
844
845
846
847
848
849
850
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (recvBufferSize[partition] > 0) {
	  mpiComm.Recv(recvBuffer[partition],
			 recvBufferSize[partition],
			 MPI_DOUBLE,
			 partition,
			 0);
851
852
853
854
855
	}
      }
    }    

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

    // === copy received values into DOFVector ===
859
860
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
861
	std::vector<DegreeOfFreedom>::iterator dofIt = recvOrder[partition].begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
862
	for (int i = 0; i < recvBufferSize[partition]; i++) {
863
864
865
866
867
868
869
	  (*values)[*dofIt] = recvBuffer[partition][i];
	  ++dofIt;
	}
      }
    }    

    // === free send and receive buffers ===
870
871
872
    for (partition = 0; partition < mpiSize; partition++) {
      if (partition != mpiRank) {
	if (sendBufferSize[partition] > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
873
	  delete [] sendBuffer[partition];
874
	if (recvBufferSize[partition] > 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
875
	  delete [] recvBuffer[partition];
876
877
878
879
880
881
882
883
884
885
      }
    }
  }

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

886
    const FiniteElemSpace *feSpace = globalSolution->getFeSpace();
887
    int dim = mesh->getDim();
888
889
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int numFcts = basFcts->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
890
891
892
    DegreeOfFreedom *coarseDOFs = new DegreeOfFreedom[numFcts];
    DegreeOfFreedom *fineDOFs = new DegreeOfFreedom[numFcts];
    DOFAdmin *admin = feSpace->getAdmin();
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
    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;

912
    MyDualTraverse dualTraverse(localCoarseGridLevel_);
913
914
915
    ElInfo *elInfo1, *elInfo2;
    ElInfo *large, *small;

916
917
918
919
920
921
922
923
924
    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);
925
926
927
928
929
930
931

    while (cont) {
      Element *element1 = elInfo1->getElement();
      Element *element2 = elInfo2->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>
	(element1->getElementData(PARTITION_ED));
Naumann, Andreas's avatar
Naumann, Andreas committed
932
933
934
935
936
937
#ifdef DEBUG
	if(element1->getElementData()==NULL)
		std::cout<<"even no elementData, id:"<<element1->getIndex()<<std::endl;
	if(partitionData==NULL) std::cout<<"elem_id:"<<element1->getIndex()<<std::endl;
#endif
	TEST_EXIT_DBG(partitionData!=NULL)("PARTITION_ED data not found\n");
938
      if (partitionData->getPartitionStatus() == IN) {
939
940

	// get coarse dofs
941
	if (element1 != lastCoarseElement) {
942
943
944
945
	  basFcts->getLocalIndices(element1, admin, coarseDOFs);
	  lastCoarseElement = element1;
	}
      
946
	if (elementPartitions_[element1].size() > 1) {
947
948
949
950
	  // get fine dofs
	  basFcts->getLocalIndices(element2, admin, fineDOFs);
      
	  // for all fine DOFs
951
952
	  for (int i = 0; i < numFcts; i++) {
	    if (!visited[fineDOFs[i]]) {
953
954
	      visited[fineDOFs[i]] = true;

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

	      // for all coarse vertex DOFs
959
	      for (int j = 0; j < dim + 1; j++) {
960
961
		partBegin = vertexPartitions[coarseDOFs[j]].begin();
		partEnd = vertexPartitions[coarseDOFs[j]].end();
962
		for (partIt = partBegin; partIt != partEnd; ++partIt) {
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
		  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);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
978
979
    delete [] coarseDOFs;
    delete [] fineDOFs;
980
981
982
983
984
985

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

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

986
    for (wIt = wBegin; wIt != wEnd; ++wIt) {
987
988
989
990
      DegreeOfFreedom dof = wIt->first;
      (*globalSolution)[dof] = 0.0;
    }
    
991
    for (wIt = wBegin; wIt != wEnd; ++wIt) {
992
993
994
995
996
      DegreeOfFreedom dof = wIt->first;
      std::map<int, double>::iterator partIt, partBegin, partEnd;
      partBegin = wIt->second.begin();
      partEnd = wIt->second.end();
    
997
      for (partIt = partBegin; partIt != partEnd; ++partIt) {
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
	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;

1012
1013
    if (!add) 
      errVec.clear();
1014
1015

    TraverseStack stack;
1016
    ElInfo *elInfo = stack.traverseFirst(mesh, 
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
					 -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;

Thomas Witkowski's avatar
Thomas Witkowski committed
1049
    elemWeights.clear();
1050
1051

    TraverseStack stack;
1052
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
1053
					 Mesh::CALL_EVERY_EL_PREORDER);
Thomas Witkowski's avatar
Thomas Witkowski committed
1054
    while (elInfo) {
1055
1056
1057
1058
1059
1060
      Element *element = elInfo->getElement();

      // get partition data
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));

Thomas Witkowski's avatar
Thomas Witkowski committed
1061
1062
      if (partitionData && partitionData->getPartitionStatus() == IN) {
	if (partitionData->getLevel() == 0) {
1063
1064
1065
	  elNum = element->getIndex();
	}
	TEST_EXIT(elNum != -1)("invalid element number\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
1066
1067
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }


  void ParallelProblem::globalRefinements() 
  {
    if (globalRefinements_ <= 0) 
      return;

    TraverseStack stack;
1085
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
    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);
    }

1098
    refinementManager->refineMesh(mesh);
1099
1100
1101
1102
1103
1104
  }


  void ParallelProblem::createOverlap(int level, int overlap, bool openOverlap,
				      std::map<Element*, int> &overlapDistance)
  {
1105
    int dim = mesh->getDim();
1106
1107
1108
1109
1110
1111

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

    TraverseStack stack;
1112
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
1113
    while (elInfo) {
1114
1115
1116
1117
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
      
1118
1119
1120
      if (partitionData) {
	if (partitionData->getLevel() == level) {
	  for (int i = 0; i < dim + 1; i++) {	    
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
	    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;

1138
    while (!overlapQueue.empty()) {
1139
1140
1141
1142
1143
1144
1145
1146
      // get first element in queue
      Element *element = overlapQueue.front();
      overlapQueue.pop();

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

1147
1148
      if (distance >= overlap) 
	continue;
1149
1150
      
      // get all adjacent elements
1151
      for (int i = 0; i < dim + 1; i++) {
1152
1153
	itBegin = vertexElements[element->getDOF(i, 0)].begin();
	itEnd = vertexElements[element->getDOF(i, 0)].end();
1154
1155
	for (it = itBegin; it != itEnd; ++it) {
	  if (overlapDistance[*it] == -1) {
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
	    // 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)
  {
1176
    int dim = mesh->getDim();
1177
1178

    // clear partition dof vector
1179
    vertexPartitions.clear();
1180
1181

    // first: partition elements ...
1182
1183
1184
    int partition;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
1185
    while (elInfo) {
1186
1187
1188
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
1189
1190
      if (partitionData) {
	if (partitionData->getLevel() == 0) {
1191
	  partition = partitionVec[element->getIndex()];
1192
	}
1193
1194
	if (partitionData->getLevel() == level) {
	  for (int i = 0; i < dim + 1; i++) {
1195
	    vertexPartitions[element->getDOF(i, 0)].insert(partition);
1196
1197
1198
1199
1200
1201
	  }
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }

1202
    if (overlap > 1 || openOverlap == false) {
1203
      // exchange mesh structure codes
Thomas Witkowski's avatar
Thomas Witkowski committed
1204
      MeshStructure *structures = new MeshStructure[mpiSize];
1205
1206
1207
      exchangeMeshStructureCodes(structures);

      // merge codes
1208
1209
      for (int rank = 0; rank < mpiSize; rank++) {
	if (rank != mpiRank) {
1210
	  structures[mpiRank].merge(&structures[rank]);
1211
1212
1213
	}
      }
    
1214
      MeshStructure &compositeStructure = structures[mpiRank];
1215
1216
1217
1218
1219
1220
      compositeStructure.reset();

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

1221
1222
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
      while (elInfo) {
1223
1224
1225
1226
	Element *element = elInfo->getElement();
	PartitionElementData *partitionData = 
	  dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
	
1227
	if (partitionData && partitionData->getLevel() == level) {
1228
1229
	  int compositeI