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
  PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
24
  {    
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
  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);
384
385
    TEST_EXIT(size == nRankDOFs * nComponents)
      ("Vector and rank DOFs does not fit together!\n");
386
387
388
389
390
391
392
393
394
395
396
397
398
#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);

399

400
401
402
403
404
405
406
407
408
    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++) {
409
410
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs * nComponents];
411
412
413
414
415
416
417
418

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

419
420
      request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs * nComponents,
						MPI_DOUBLE, sendIt->first, 0);
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
    }

    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);
444
445
 	for (int k = 0; k < nRecvDOFs; k++)
 	  (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++];
446
447
448
449
450
451
452
453
454
455
      }

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


456
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
  {
    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
477
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
478
479
480
481
482
483
484
485
486
487
488
489
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

490

491
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
492
493
494
495
496
497
498
499
500
501
502
503
504
  {
    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);
  }

505

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

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

512
    TraverseStack stack;
513
514
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
515
516
517
518
519
520
521
522
523
524
525
526
    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));
527

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

531
532
	    // === 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. ===
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557

	    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;

558
	    // === And add the part of the interior boundary. ===
559
560

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

565
566
567
568
569
570
	    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;
571
572
573
574
575
576
 	  }
	}
      }

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


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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
611
612
613

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

Thomas Witkowski's avatar
Thomas Witkowski committed
614
615

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

      // === 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
632
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
633
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
634
635
636
637
638
639
640
641
642
643
644
645
646

	  // 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];
647
648
649
  }


650
  void ParallelDomainBase::removeMacroElements()
651
652
653
654
655
656
657
658
659
660
661
662
  {
    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
663
    mesh->removeMacroElements(macrosToRemove, feSpace);
664
665
666
  }


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

    // === Get rank information about DOFs. ===
674
675
676

    // 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
677
    DofContainer rankAllDofs;
678
    DofToRank boundaryDofs;
679

680
    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDofs);
681
682
683
684

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

685

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

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

692
693
694
695
    
    // === 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
696

Thomas Witkowski's avatar
Thomas Witkowski committed
697

698
699
700
    // 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
701
702
703
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
704
705
706
      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
707
      isRankDof[i] = true;
708
      i++;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
709
710
    }

711
712
713

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
727
 
728
729
730
731
732
733
    // === 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;
734

735
736
737
738
    // 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;
739

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

      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
749
	  if (*itRanks != mpiRank) {
750
	    TEST_EXIT_DBG(rankDofsNewGlobalIndex.count(it->first) == 1)
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
751
752
	      ("DOF Key not found!\n");

753
	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
754
	  }
755
756
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
757
758
	// 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.
759
760
761
762
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
763
764
765
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
766

767
    // === Send and receive the dof indices at boundary. ===
768

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

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

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

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

793
	sendDofs[sendIt->first].push_back(dofIt->first);
794
795
      }

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

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

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


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

    
815
    // === Delete send buffers. ===
816

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


821
    // === Change dof indices at boundary from other ranks. ===
822

Thomas Witkowski's avatar
Thomas Witkowski committed
823
824
    // 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
825
    // 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
826
827
828
829
830
    // 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;
831
    for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
832
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
833
834
      dofChanged[dofIt->first] = false;

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

840
      for (int j = 0; j < recvIt->second; j++) {
841

842
843
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
844

845
	bool found = false;
846

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

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

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

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

859
	    found = true;
860
861
862
	    break;
	  }
	}
863

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

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

870
    // === Create now the local to global index and local to dof index mappings.  ===
871

872
873
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
874
875
876
  }


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

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

883
    // === Get all DOFs in ranks partition. ===
884

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

      elInfo = stack.traverseNext(elInfo);
    }

899
900
901
902
903
904
905
    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
906
    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
907
    // === rankDOFs to boundaryDOFs.                                               ===
908

909
910
    sendDofs.clear();
    recvDofs.clear();
911

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

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

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

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

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

943
944
945
946
947
948
949
	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]);
950
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
951
	
952
953
954
      }
    }    

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

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

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

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

981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
	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())
997
998
	    ("Should never happen!\n");

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

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

1008
1009
1010
1011
1012
    nRankDOFs = rankDOFs.size();

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

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


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

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


1022
1023
1024
    // Do not change the indices now, but create a new indexing a store it here.
    DofIndexMap rankDofsNewLocalIndex;
    isRankDof.clear();
1025
    int i = 0;
1026
1027
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
1028

1029
1030
1031
      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
1032
      isRankDof[i] = true;
1033
      i++;
1034
1035
    }

1036
    // Stores for all rank owned dofs a new global index.
1037
    DofIndexMap rankDofsNewGlobalIndex;
1038
1039
1040
1041
1042
1043
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
1044
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
1045
1046
1047
1048
1049
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
      i++;
    }


1050
1051
    // === Send new DOF indices. ===

1052
1053
    std::vector<int*> sendBuffers(sendDofs.size());
    std::vector<int*> recvBuffers(recvDofs.size());
1054

1055
    MPI::Request request[sendDofs.size() + recvDofs.size()];
Thomas Witkowski's avatar
Thomas Witkowski committed
1056
1057
    int requestCounter = 0;

1058
    i = 0;
1059
1060
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt, i++) {