ParallelDomainBase.cc 47.3 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
#include <algorithm>

3
#include "ParallelDomainBase.h"
4
5
6
7
8
9
10
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
11
12
#include "DOFMatrix.h"
#include "DOFVector.h"
13
#include "SystemVector.h"
14
#include "VtkWriter.h"
15
#include "ElementDofIterator.h"
16
17
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
18
19

#include "petscksp.h"
20
21
22

namespace AMDiS {

23
24
  PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
  {
25
    if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
26
      std::cout << "  Iteration " << iter << ": " << rnorm << std::endl;
27
28
29
30

    return 0;
  }

31
32
33
34
35
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

36
37
38
39
40
  ParallelDomainBase::ParallelDomainBase(const std::string& name,
					 ProblemIterationInterface *iIF,
					 ProblemTimeInterface *tIF,
					 FiniteElemSpace *fe,
					 RefinementManager *refineManager)
41
42
    : iterationIF(iIF),
      timeIF(tIF),
43
44
      feSpace(fe),
      mesh(fe->getMesh()),
45
      refinementManager(refineManager),
46
      initialPartitionMesh(true),
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
47
      nRankDOFs(0),
48
49
      rstart(0),
      nComponents(1)
50
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
51
52
53
54
55
56
57
    FUNCNAME("ParallelDomainBase::ParalleDomainBase()");

    TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
      ("Only meshes with one DOFAdmin are supported!\n");
    TEST_EXIT(mesh->getDOFAdmin(0).getNumberOfPreDOFs(0) == 0)
      ("Wrong pre dof number for DOFAdmin!\n");

58
59
60
61
62
63
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
  }

64
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
65
66
67
68
  {
    if (mpiSize <= 1)
      return;

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

74
75
76
77
78
79
80
    // create an initial partitioning of the mesh
    partitioner->createPartitionData();
    // set the element weights, which are 1 at the very first begin
    setElemWeights(adaptInfo);
    // and now partition the mesh
    partitionMesh(adaptInfo);   

Thomas Witkowski's avatar
Thomas Witkowski committed
81
82
83
84
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
    DbgCreateElementMap(elMap);
#endif
85

86
    // === Create new global and local DOF numbering. ===
87

88
89
90
    // Set of all DOFs of the rank.
    std::vector<const DegreeOfFreedom*> rankDOFs;
    // Number of DOFs in ranks partition that are owned by the rank.
Thomas Witkowski's avatar
Thomas Witkowski committed
91
    nRankDOFs = 0;
92
    // Number of all DOFs in the macro mesh.
93
    int nOverallDOFs = 0;
94

95
    createLocalGlobalNumbering(rankDOFs, nRankDOFs, nOverallDOFs);
96

Thomas Witkowski's avatar
Thomas Witkowski committed
97
98
    // === Create interior boundary information ===

99
    createInteriorBoundaryInfo(rankDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
100

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
101
102
103
104
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

Thomas Witkowski's avatar
Thomas Witkowski committed
105
#if (DEBUG != 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
106
    DbgTestElementMap(elMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
107
108
    DbgTestInteriorBoundary();
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
109

110

111
    // === Reset all DOFAdmins of the mesh. ===
112

113
    updateDofAdmins();
114

115

116
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
117

Thomas Witkowski's avatar
Thomas Witkowski committed
118
119
120
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

Thomas Witkowski's avatar
Thomas Witkowski committed
121
122
    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);
123

124
125
126
127
128
#if (DEBUG != 0)
      elMap.clear();
      DbgCreateElementMap(elMap);
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
129
      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
130
131
132
133
134
135

      updateDofAdmins();

#if (DEBUG != 0)
      DbgTestElementMap(elMap);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
136
137
    }

138

Thomas Witkowski's avatar
Thomas Witkowski committed
139
#if (DEBUG != 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
140
    DbgTestCommonDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
141
#endif
142

143

144
    // === Create petsc matrix. ===
145

146
147
    int nRankRows = nRankDOFs * nComponents;
    int nOverallRows = nOverallDOFs * nComponents;
148

149
150
151
    MatCreate(PETSC_COMM_WORLD, &petscMatrix);
    MatSetSizes(petscMatrix, nRankRows, nRankRows, nOverallRows, nOverallRows);
    MatSetType(petscMatrix, MATAIJ);
152

153
154
155
156
157
158
159
    VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
    VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
    VecSetType(petscRhsVec, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &petscSolVec);
    VecSetSizes(petscSolVec, nRankRows, nOverallRows);
    VecSetType(petscSolVec, VECMPI);
160
161
  }

162
163
164

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
  {
165
166
167
    MatDestroy(petscMatrix);
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
168
  }
169

170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
  
  void ParallelDomainBase::updateDofAdmins()
  {
    int nAdmins = mesh->getNumberOfDOFAdmin();
    for (int i = 0; i < nAdmins; i++) {
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
      for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
 	admin.setDOFFree(j, false);

      admin.setUsedSize(mapLocalGlobalDOFs.size());
      admin.setUsedCount(mapLocalGlobalDOFs.size());
      admin.setFirstHole(mapLocalGlobalDOFs.size());
    }
  }

188

189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
  void ParallelDomainBase::testForMacroMesh()
  {
    FUNCNAME("ParallelDomainBase::testForMacroMesh()");

    int nMacroElements = 0;

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

      elInfo = stack.traverseNext(elInfo);
    }

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

210
211
  void ParallelDomainBase::setDofMatrix(DOFMatrix* mat, int dispMult, 
					int dispAddRow, int dispAddCol)
212
  {
213
214
215
216
    FUNCNAME("ParallelDomainBase::setDofMatrix()");

    TEST_EXIT(mat)("No DOFMatrix!\n");

217
218
219
220
221
222
223
224
225
226
227
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    traits::row<Matrix>::type row(mat->getBaseMatrix());
    traits::col<Matrix>::type col(mat->getBaseMatrix());
    traits::const_value<Matrix>::type value(mat->getBaseMatrix());

    typedef traits::range_generator<major, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

228
229
230
231
    for (cursor_type cursor = begin<major>(mat->getBaseMatrix()), 
	   cend = end<major>(mat->getBaseMatrix()); cursor != cend; ++cursor)
      for (icursor_type icursor = begin<nz>(cursor), 
	     icend = end<nz>(cursor); icursor != icend; ++icursor)
232
	if (value(*icursor) != 0.0) {
233
234
 	  int r = mapLocalGlobalDOFs[row(*icursor)] * dispMult + dispAddRow;
 	  int c = mapLocalGlobalDOFs[col(*icursor)] * dispMult + dispAddCol;  
235
	  double v = value(*icursor);
236

237
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
238
	}
239
  }
240

241

242
243
244
  void ParallelDomainBase::setDofVector(DOFVector<double>* vec, 
					int dispMult, int dispAdd)
  {
245
246
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
247
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()] * dispMult + dispAdd;
248
249
250
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
251
    }    
252
253
  }

254

255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
  void ParallelDomainBase::fillPetscMatrix(DOFMatrix *mat, DOFVector<double> *vec)
  {
    setDofMatrix(mat);

    MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);

    setDofVector(vec);   
  }

  
  void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
  {
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
	if ((*mat)[i][j]) 
	  setDofMatrix((*mat)[i][j], nComponents, i, j);

    MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);

    for (int i = 0; i < nComponents; i++)
      setDofVector(vec->getDOFVector(i), nComponents, i);
  }


285
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
286
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

289
290
291
292
293
294
    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
    KSPGetPC(ksp, &pc);
295
296
    //    PCSetType(pc, PCNONE);
    PCSetType(pc, PCJACOBI);
297
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
298
299
    KSPSetType(ksp, KSPBCGS);
    //KSPSetType(ksp, KSPCG);
300
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
301
302
    KSPSolve(ksp, petscRhsVec, petscSolVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
305
306
307
308
#if (DEBUG != 0)
    int size = 0;
    VecGetLocalSize(petscSolVec, &size);
    TEST_EXIT(size == nRankDOFs)("Vector and rank DOFs does not fit together!\n");
#endif

309
310
311
    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

Thomas Witkowski's avatar
Thomas Witkowski committed
312
    for (int i = 0; i < nRankDOFs; i++)
313
      (*vec)[mapLocalToDofIndex[i]] = vecPointer[i];
314
315
316

    VecRestoreArray(petscSolVec, &vecPointer);

317
318
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
319
320
321

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
322
323
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
324
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
325
	 sendIt != sendDofs.end(); ++sendIt, i++) {
326
327
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
328

329
      for (int j = 0; j < nSendDOFs; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
330
	sendBuffers[i][j] = (*vec)[*((sendIt->second)[j])];
331

Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
334
335
336
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
337
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
338
	 recvIt != recvDofs.end(); ++recvIt, i++) {
339
340
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
341

Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
344
345
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
346
347
348

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

349
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
350
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
351
	 recvIt != recvDofs.end(); ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
352
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
353
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
354
355
356
357

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

362

363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
  void ParallelDomainBase::solvePetscMatrix(SystemVector *vec)
  {
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
    KSPGetPC(ksp, &pc);
    //    PCSetType(pc, PCNONE);
    PCSetType(pc, PCJACOBI);
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
    //KSPSetType(ksp, KSPCG);
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
    KSPSolve(ksp, petscRhsVec, petscSolVec);

#if (DEBUG != 0)
    int size = 0;
    VecGetLocalSize(petscSolVec, &size);
    TEST_EXIT(size == nRankDOFs)("Vector and rank DOFs does not fit together!\n");
#endif

    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

    for (int i = 0; i < nComponents; i++) {
      DOFVector<double> *dofvec = vec->getDOFVector(i);
      for (int j = 0; j < nRankDOFs; j++)
	(*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];
    }

    VecRestoreArray(petscSolVec, &vecPointer);

    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());

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

      int counter = 0;
      for (int j = 0; j < nComponents; j++) {
	DOFVector<double> *dofvec = vec->getDOFVector(j);
	for (int k = 0; k < nSendDOFs; k++)
	  sendBuffers[i][counter++] = (*dofvec)[*((sendIt->second)[k])];
      }

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

    i = 0;
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt, i++) {
      int nRecvDOFs = recvIt->second.size() * nComponents;
      recvBuffers[i] = new double[nRecvDOFs];

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


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

    i = 0;
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt, i++) {
      int nRecvDOFs = recvIt->second.size();

      int counter = 0;
      for (int j = 0; j < nComponents; j++) {
	DOFVector<double> *dofvec = vec->getDOFVector(j);
	for (int k = 0; k < nRecvDOFs; k++)
	  (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++];
      }

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


454
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
  {
    double localWeightSum = 0.0;
    int elNum = -1;

    elemWeights.clear();

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

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

      if (partitionData && partitionData->getPartitionStatus() == IN) {
	if (partitionData->getLevel() == 0) {
	  elNum = element->getIndex();
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
475
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
476
477
478
479
480
481
482
483
484
485
486
487
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

488

489
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
490
491
492
493
494
495
496
497
498
499
500
501
502
  {
    if (initialPartitionMesh) {
      initialPartitionMesh = false;
      partitioner->fillCoarsePartitionVec(&oldPartitionVec);
      partitioner->partition(&elemWeights, INITIAL);
    } else {
      oldPartitionVec = partitionVec;
      partitioner->partition(&elemWeights, ADAPTIVE_REPART, 100.0 /*0.000001*/);
    }    

    partitioner->fillCoarsePartitionVec(&partitionVec);
  }

503

504
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs)
505
  {
506
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
507
508
509

    // === First, create all the information about the interior boundaries. ===

510
    TraverseStack stack;
511
512
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
513
514
515
516
517
518
519
520
521
522
523
524
    while (elInfo) {
      Element *element = elInfo->getElement();

      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
      if (partitionData->getPartitionStatus() == IN) {
	for (int i = 0; i < 3; i++) {
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->getElementData(PARTITION_ED));
525

526
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
527
528
	    // We have found an element that is at an interior boundary. 

529
530
	    // === Find out, if the boundary part of the element corresponds to the  ===
	    // === rank or to the rank "on the other side" of the interoir boundary. ===
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555

	    const DegreeOfFreedom* boundDOF1 = NULL;
	    const DegreeOfFreedom* boundDOF2 = NULL;
	    
	    switch (i) {
	    case 0:
	      boundDOF1 = element->getDOF(1);
	      boundDOF2 = element->getDOF(2);
	      break;
	    case 1:
	      boundDOF1 = element->getDOF(0);
	      boundDOF2 = element->getDOF(2);
	      break;
	    case 2:
	      boundDOF1 = element->getDOF(0);
	      boundDOF2 = element->getDOF(1);
	      break;
	    default:
	      ERROR_EXIT("Should never happen!\n");
	    }

	    bool isRankDOF1 = (find(rankDOFs.begin(), rankDOFs.end(), boundDOF1) != rankDOFs.end());
	    bool isRankDOF2 = (find(rankDOFs.begin(), rankDOFs.end(), boundDOF2) != rankDOFs.end());
	    bool ranksBoundary = isRankDOF1 || isRankDOF2;

556
	    // === And add the part of the interior boundary. ===
557
558

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
559
560
561
562
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

563
564
565
566
567
568
	    bound.rankObject.el = element;
	    bound.rankObject.subObjAtBoundary = EDGE;
	    bound.rankObject.ithObjAtBoundary = i;
	    bound.neighbourObject.el = elInfo->getNeighbour(i);
	    bound.neighbourObject.subObjAtBoundary = EDGE;
	    bound.neighbourObject.ithObjAtBoundary = -1;
569
570
571
572
573
574
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
575
576
577


    // === Once we have this information, we must care about their order. Eventually ===
Thomas Witkowski's avatar
Thomas Witkowski committed
578
579
    // === all the boundaries have to be in the same order on both ranks that share  ===
    // === the bounday.                                                              ===
Thomas Witkowski's avatar
Thomas Witkowski committed
580

Thomas Witkowski's avatar
Thomas Witkowski committed
581
    std::vector<int*> sendBuffers, recvBuffers;
Thomas Witkowski's avatar
Thomas Witkowski committed
582

Thomas Witkowski's avatar
Thomas Witkowski committed
583
584
585
586
    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
587
588
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
589
590
591
      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
      for (int i = 0; i < nSendInt; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
592
593
594
	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
595
596
      request[requestCounter++] =
	mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
597
598
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
599
600
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
601
602
      int nRecvInt = rankIt->second.size();
      int *buffer = new int[nRecvInt];
Thomas Witkowski's avatar
Thomas Witkowski committed
603
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
604

Thomas Witkowski's avatar
Thomas Witkowski committed
605
606
      request[requestCounter++] = 
	mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
607
608
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
609
610
611

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

Thomas Witkowski's avatar
Thomas Witkowski committed
612
613

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
614
615
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
616
617
618
619
620
621
622
623
624
625
626
627
628
629

      // === We have received from rank "rankIt->first" the ordered list of element ===
      // === indices. We now have to sort the corresponding list in this rank to    ===
      // === get the same order.                                                    ===
      
      for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
	// If the expected object is not at place, search for it.
	if ((rankIt->second)[j].neighbourObject.el->getIndex() != recvBuffers[i][j]) {
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
	    if ((rankIt->second)[k].neighbourObject.el->getIndex() == recvBuffers[i][j])
	      break;

	  // The element must always be found, because the list is just in another order.
Thomas Witkowski's avatar
Thomas Witkowski committed
630
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
631
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
634
635
636
637
638
639
640
641
642
643
644

	  // Swap the current with the found element.
	  AtomicBoundary tmpBound = (rankIt->second)[k];
	  (rankIt->second)[k] = (rankIt->second)[j];
	  (rankIt->second)[j] = tmpBound;	
	}
      }

      delete [] recvBuffers[i++];
    }

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


648
  void ParallelDomainBase::removeMacroElements()
649
650
651
652
653
654
655
656
657
658
659
660
  {
    std::vector<MacroElement*> macrosToRemove;
    for (std::deque<MacroElement*>::iterator it = mesh->firstMacroElement();
	 it != mesh->endOfMacroElements();
	 ++it) {
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>
	((*it)->getElement()->getElementData(PARTITION_ED));
      if (partitionData->getPartitionStatus() != IN)
	macrosToRemove.push_back(*it);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
661
    mesh->removeMacroElements(macrosToRemove, feSpace);
662
663
664
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
665
666
667
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
						      int& nRankDOFs, 
						      int& nOverallDOFs)
668
  {
669
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
670
671

    // === Get rank information about DOFs. ===
672
673
674

    // Stores to each DOF pointer the set of ranks the DOF is part of.
    std::map<const DegreeOfFreedom*, std::set<int> > partitionDOFs;
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
675
    DofContainer rankAllDofs;
676
    DofToRank boundaryDofs;
677

678
    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDofs);
679
680
681
682

    nRankDOFs = rankDOFs.size();
    nOverallDOFs = partitionDOFs.size();

683

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

686
    rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
687
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
688
689
    rstart -= nRankDOFs;

690
691
692
693
    
    // === Create for all dofs in rank new indices. The new index must start at ===
    // === index 0, must be continues and have the same order as the indices    ===
    // === had before.                                                          ===
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
694

Thomas Witkowski's avatar
Thomas Witkowski committed
695

696
697
698
    // Do not change the indices now, but create a new indexing a store it here.
    DofIndexMap rankDofsNewLocalIndex;
    isRankDof.clear();
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
699
700
701
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
702
703
704
      rankDofsNewLocalIndex[*dofIt] = i;
      // First, we set all dofs in ranks partition to be owend by the rank. Later,
      // the dofs in ranks partition that are owned by other rank are set to false.
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
705
      isRankDof[i] = true;
706
      i++;
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
707
708
    }

709
710
711

    // === Create for all rank owned dofs a new global indexing.                ===

712
713
    // Stores for dofs in rank a new global index.
    DofIndexMap rankDofsNewGlobalIndex;
714
715
716
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
717
718
719
    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
720
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
721
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
722
723
724
      i++;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
725
 
726
727
728
729
730
731
    // === Create information which dof indices must be send and which must  ===
    // === be received.                                                      ===

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

733
734
735
736
    // Maps to each rank the number of new DOF indices this rank will receive from.
    // All these DOFs are at an interior boundary on this rank, but are owned by
    // another rank.
    std::map<int, int> recvNewDofs;
737

738
    for (DofToRank::iterator it = boundaryDofs.begin(); it != boundaryDofs.end(); ++it) {
739
740
741
742
743
744
745
746

      if (it->second == mpiRank) {
	// If the boundary dof is a rank dof, it must be send to other ranks.

	// Search for all ranks that have this dof too.
	for (std::set<int>::iterator itRanks = partitionDOFs[it->first].begin();
	     itRanks != partitionDOFs[it->first].end();
	     ++itRanks) {
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
747
	  if (*itRanks != mpiRank) {
748
	    TEST_EXIT_DBG(rankDofsNewGlobalIndex.count(it->first) == 1)
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
749
750
	      ("DOF Key not found!\n");

751
	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
752
	  }
753
754
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
755
756
	// If the boundary dof is not a rank dof, its new dof index (and later
	// also the dof values) must be received from another rank.
757
758
759
760
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
761
762
763
      }
    }

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
764

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

Thomas Witkowski's avatar
Thomas Witkowski committed
767
    std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size());
768
769
770

    sendDofs.clear();
    recvDofs.clear();
771
    
Thomas Witkowski's avatar
Thomas Witkowski committed
772
773
774
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
775
    i = 0;
776
777
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
778
779
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
780
781
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
782

783
      int c = 0;
784
785
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
786
	   dofIt != sendIt->second.end();
787
	   ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
788
	sendBuffers[i][c++] = *(dofIt->first);
789
	sendBuffers[i][c++] = dofIt->second;
790

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

Thomas Witkowski's avatar
Thomas Witkowski committed
794
795
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
796
797
798
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
805
806
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
807
808
809
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
810
    MPI::Request::Waitall(requestCounter, request);
811
812

    
813
    // === Delete send buffers. ===
814

Thomas Witkowski's avatar
Thomas Witkowski committed
815
816
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
817
818


819
    // === Change dof indices at boundary from other ranks. ===
820

Thomas Witkowski's avatar
Thomas Witkowski committed
821
822
    // Within this small data structure we track which dof index was already changed.
    // This is used to avoid the following situation: Assume, there are two dof indices
823
    // a and b in boundaryDofs. Then we have to change index a to b and b to c. When
Thomas Witkowski's avatar
Thomas Witkowski committed
824
825
826
827
828
    // the second rule applies, we have to avoid that not the first b, resulted from
    // changing a to b, is set to c, but the second one. Therefore, after the first
    // rule was applied, the dof pointer is set to false in this data structure and 
    // is not allowed to be changed anymore.
    std::map<const DegreeOfFreedom*, bool> dofChanged;
829
    for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
830
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
831
832
      dofChanged[dofIt->first] = false;

833
    i = 0;
834
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
835
836
837
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

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

840
841
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
842

843
	bool found = false;
844

Thomas Witkowski's avatar
Thomas Witkowski committed
845
846
	// Iterate over all boundary dofs to find the dof, which index we have to change.

847
848
	for (DofToRank::iterator dofIt = boundaryDofs.begin(); 
	     dofIt != boundaryDofs.end(); ++dofIt) {
849

Thomas Witkowski's avatar
Thomas Witkowski committed
850
851
852
	  if (*(dofIt->first) == oldDof && !dofChanged[dofIt->first]) {
	    dofChanged[dofIt->first] = true;

853
	    recvDofs[recvIt->first].push_back(dofIt->first);
854
	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
855
	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
856

857
	    found = true;
858
859
860
	    break;
	  }
	}
861

Thomas Witkowski's avatar
Thomas Witkowski committed
862
	TEST_EXIT_DBG(found)("Should not happen!\n");
863
864
865
866
      }

      delete [] recvBuffers[i];
    }
Thomas Witkowski's avatar
n  
Thomas Witkowski committed
867

868
    // === Create now the local to global index and local to dof index mappings.  ===
869

870
871
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
872
873
874
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
875
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs)
876
  {
877
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
878

879
    typedef std::set<const DegreeOfFreedom*> DofSet;
Thomas Witkowski's avatar
Thomas Witkowski committed
880

881
    // === Get all DOFs in ranks partition. ===
882

Thomas Witkowski's avatar
Thomas Witkowski committed
883
    ElementDofIterator elDofIt(feSpace);
884
    DofSet rankDOFSet;
885
886
887
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
888
      Element *element = elInfo->getElement();     
Thomas Witkowski's avatar
Thomas Witkowski committed
889
890
      elDofIt.reset(element);
      do {
891
	rankDOFSet.insert(elDofIt.getDofPtr());
Thomas Witkowski's avatar
Thomas Witkowski committed
892
      } while(elDofIt.next());
893
894
895
896

      elInfo = stack.traverseNext(elInfo);
    }

897
898
899
900
901
902
903
    DofContainer rankAllDofs;
    for (DofSet::iterator dofIt = rankDOFSet.begin(); dofIt != rankDOFSet.end(); ++dofIt)
      rankAllDofs.push_back(*dofIt);
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
    DofContainer rankDOFs = rankAllDofs;


Thomas Witkowski's avatar
Thomas Witkowski committed
904
    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
905
    // === rankDOFs to boundaryDOFs.                                               ===
906

907
908
    sendDofs.clear();
    recvDofs.clear();
909

Thomas Witkowski's avatar
Thomas Witkowski committed
910
911
912
    for (RankToBoundMap::iterator it =  myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end(); ++it) {

Thomas Witkowski's avatar
Thomas Witkowski committed
913
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
914
915
	   boundIt != it->second.end(); ++boundIt) {

916
917
	DofContainer dofs;
	DofContainer &dofsToSend = sendDofs[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
918
919
920

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
921
922
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
923
924
	  break;
	case 1:
925
926
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
927
928
	  break;
	case 2:
929
930
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
Thomas Witkowski's avatar
Thomas Witkowski committed
931
932
933
934
935
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

936
937
938
939
	for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt) {
	  if (find(dofsToSend.begin(), dofsToSend.end(), *dofIt) == dofsToSend.end())
	    dofsToSend.push_back(*dofIt);
	}
940

941
942
943
944
945
946
947
	dofs.clear();
	addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
			 dofs);	
  	addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
  		       dofs);
	for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
	  dofsToSend.push_back(dofs[i]);
948
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
949
	
950
951
952
      }
    }    

Thomas Witkowski's avatar
Thomas Witkowski committed
953
954
955
    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end(); ++it) {

956
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
957
958
	   boundIt != it->second.end(); ++boundIt) {

959
960
	DofContainer dofs;
	DofContainer &dofsToRecv = recvDofs[it->first];
961
962
963

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
964
965
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
966
967
	  break;
	case 1:
968
969
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
970
971
	  break;
	case 2:
972
973
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
974
975
976
977
978
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
	for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt) {
	  DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), *dofIt);
	  if (eraseIt != rankDOFs.end())
	    rankDOFs.erase(eraseIt);
	  if (find(dofsToRecv.begin(), dofsToRecv.end(), *dofIt) == dofsToRecv.end())
	    dofsToRecv.push_back(*dofIt);
	}

	dofs.clear();
	addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 
		       dofs);
	addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
			 dofs);

	for (int i = static_cast<int>(dofs.size()) - 1; i >= 0; i--) {
	  TEST_EXIT_DBG(find(rankDOFs.begin(), rankDOFs.end(), dofs[i]) != rankDOFs.end())
995
996
	    ("Should never happen!\n");

997
	  DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), dofs[i]);
998
999
1000
	  if (eraseIt != rankDOFs.end())
	    rankDOFs.erase(eraseIt);

1001
	  dofsToRecv.push_back(dofs[i]);
1002
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
1003
1004
1005
      }
    }

1006
1007
1008
1009
1010
    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1011
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
1012
1013
1014
1015
1016
    rstart -= nRankDOFs;


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

Thomas Witkowski's avatar
Thomas Witkowski committed
1017
    mpiComm.Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);