ParallelDomainProblem.cc 45.9 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
3
#include <boost/lambda/lambda.hpp>
#include <algorithm>

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

#include "petscksp.h"
20
21

namespace AMDiS {
Thomas Witkowski's avatar
Thomas Witkowski committed
22
23
  
  using namespace boost::lambda;
24

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

    return 0;
  }

  ParallelDomainBase::ParallelDomainBase(const std::string& name,
					 ProblemIterationInterface *iIF,
					 ProblemTimeInterface *tIF,
					 FiniteElemSpace *fe,
					 RefinementManager *refineManager)
38
39
    : iterationIF(iIF),
      timeIF(tIF),
40
41
      feSpace(fe),
      mesh(fe->getMesh()),
42
      refinementManager(refineManager),
43
      initialPartitionMesh(true),
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
44
45
      nRankDOFs(0),
      rstart(0)
46
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
47
48
49
50
51
52
53
    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");

54
55
56
57
58
59
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
  }

60
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
61
62
63
64
  {
    if (mpiSize <= 1)
      return;

65
66
67
68
69
    // 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();

70
71
72
73
74
75
76
    // 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
77
78
79
80
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
    DbgCreateElementMap(elMap);
#endif
81

82
    // === Create new global and local DOF numbering. ===
83

84
85
86
    // 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
87
    nRankDOFs = 0;
88
    // Number of all DOFs in the macro mesh.
89
    int nOverallDOFs = 0;
90
91

    createLocalGlobalNumbering(rankDOFs, boundaryDOFs, nRankDOFs, nOverallDOFs);
92

Thomas Witkowski's avatar
Thomas Witkowski committed
93
94
    // === Create interior boundary information ===

95
    createInteriorBoundaryInfo(rankDOFs, boundaryDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
96

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
97
98
99
100
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

Thomas Witkowski's avatar
Thomas Witkowski committed
101
#if (DEBUG != 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
102
    DbgTestElementMap(elMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
103
104
    DbgTestInteriorBoundary();
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
105

106

107
    // === Reset all DOFAdmins of the mesh. ===
108

109
    updateDofAdmins();
110

111
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
112

Thomas Witkowski's avatar
Thomas Witkowski committed
113
114
115
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

Thomas Witkowski's avatar
Thomas Witkowski committed
116
117
    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);
118

119
120
121
122
123
#if (DEBUG != 0)
      elMap.clear();
      DbgCreateElementMap(elMap);
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
124
      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
125
126
127
128
129
130

      updateDofAdmins();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
133
#if (DEBUG != 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
134
    DbgTestCommonDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
135
#endif
136

137

138
    // === Create petsc matrix. ===
139

140
141
    int ierr;
    ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
142
    ierr = MatSetSizes(petscMatrix, nRankDOFs, nRankDOFs, nOverallDOFs, nOverallDOFs);
143
144
145
    ierr = MatSetType(petscMatrix, MATAIJ);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
146
    ierr = VecSetSizes(petscRhsVec, nRankDOFs, nOverallDOFs);
147
    ierr = VecSetType(petscRhsVec, VECMPI);
148
149

    ierr = VecCreate(PETSC_COMM_WORLD, &petscSolVec);
150
    ierr = VecSetSizes(petscSolVec, nRankDOFs, nOverallDOFs);
151
    ierr = VecSetType(petscSolVec, VECMPI);
152
153
  }

154
155
156
157

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
  {
  }
158

159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
  
  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());
    }
  }

177

178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
  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");
  }


  void ParallelDomainBase::fillPetscMatrix(DOFMatrix *mat, 
					   DOFVector<double> *vec)
202
  {
203
204
205
206
207
208
209
210
211
212
213
    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;

214
215
216
217
    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)
218
	if (value(*icursor) != 0.0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
219
220
 	  int r = mapLocalGlobalDOFs[row(*icursor)];
 	  int c = mapLocalGlobalDOFs[col(*icursor)];  
221
	  double v = value(*icursor);
222

223
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
224
	}
225

226
227
228
229

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

230
    
231
232
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
233
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
234
235
236
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
237
    }    
238
239
  }

240

241
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
242
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
243
244
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

245
246
247
248
249
250
    KSP ksp;
    PC pc;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
259
260
261
262
263
264
#if (DEBUG != 0)
    int size = 0;
    VecGetLocalSize(petscSolVec, &size);
    TEST_EXIT(size == nRankDOFs)("Vector and rank DOFs does not fit together!\n");
#endif

265
266
267
    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    for (int i = 0; i < nRankDOFs; i++)
269
      (*vec)[mapLocalToDofIndex[i]] = vecPointer[i];
270
271
272

    VecRestoreArray(petscSolVec, &vecPointer);

273
274
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
277

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
278
279
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); 
282
	 ++sendIt, i++) {
283
284
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
285

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

Thomas Witkowski's avatar
Thomas Witkowski committed
289
290
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
291
292
293
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
294
295
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); 
296
	 ++recvIt, i++) {
297
298
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
299

Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
302
303
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
304
305
306

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

307
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
308
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
309
310
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
311
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
312
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
313
314
315
316

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

321

322
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
  {
    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
343
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
344
345
346
347
348
349
350
351
352
353
354
355
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

356

357
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
358
359
360
361
362
363
364
365
366
367
368
369
370
  {
    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);
  }

371

Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs,
						      DofToRank& boundaryDOFs)
374
  {
375
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
376
377
378

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

379
    TraverseStack stack;
380
381
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
382
383
384
385
386
387
388
389
390
391
392
393
    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));
394

395
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
396
397
	    // We have found an element that is at an interior boundary. 

398
399
	    // === 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. ===
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

	    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;

425
	    // === And add the part of the interior boundary. ===
426
427

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
428
429
430
431
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

432
433
434
435
436
437
	    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;
438
439
440
441
442
443
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
446


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
454
455
    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
460
      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
      for (int i = 0; i < nSendInt; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
463
	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
464
465
      request[requestCounter++] =
	mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
466
467
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
      int nRecvInt = rankIt->second.size();
      int *buffer = new int[nRecvInt];
Thomas Witkowski's avatar
Thomas Witkowski committed
472
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
473

Thomas Witkowski's avatar
Thomas Witkowski committed
474
475
      request[requestCounter++] = 
	mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
478
479
480

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

Thomas Witkowski's avatar
Thomas Witkowski committed
481
482

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
483
484
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
487
488
489
490
491
492
493
494
495
496
497
498

      // === 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
499
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
500
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
501
502
503
504
505
506
507
508
509
510
511
512
513

	  // 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];
514
515
516
  }


517
  void ParallelDomainBase::removeMacroElements()
518
519
520
521
522
523
524
525
526
527
528
529
  {
    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
530
    mesh->removeMacroElements(macrosToRemove, feSpace);
531
532
533
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
534
535
536
537
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
						      DofToRank& boundaryDOFs,
						      int& nRankDOFs, 
						      int& nOverallDOFs)
538
  {
539
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
540
541

    // === Get rank information about DOFs. ===
542
543
544

    // 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
545
    DofContainer rankAllDofs;
546

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
547
    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDOFs);
548
549
550
551

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

552

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

555
    rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
556
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
557
558
    rstart -= nRankDOFs;

559
560
561
562
    
    // === 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
563

Thomas Witkowski's avatar
Thomas Witkowski committed
564

565
566
567
    // 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
568
569
570
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
571
572
573
      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
574
      isRankDof[i] = true;
575
      i++;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
576
577
    }

578
579
580

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

581
582
    // Stores for dofs in rank a new global index.
    DofIndexMap rankDofsNewGlobalIndex;
583
584
585
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
586
587
588
    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
589
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
590
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
591
592
593
      i++;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
594
 
595
596
597
598
599
600
    // === 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;
601

602
603
604
605
    // 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;
606

Thomas Witkowski's avatar
Thomas Witkowski committed
607
    for (DofToRank::iterator it = boundaryDOFs.begin(); it != boundaryDOFs.end(); ++it) {
608
609
610
611
612
613
614
615

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

620
	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
621
	  }
622
623
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
624
625
	// 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.
626
627
628
629
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
630
631
632
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
633

634
    // === Send and receive the dof indices at boundary. ===
635

Thomas Witkowski's avatar
Thomas Witkowski committed
636
    std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size());
637
638
639

    sendDofs.clear();
    recvDofs.clear();
640
    
Thomas Witkowski's avatar
Thomas Witkowski committed
641
642
643
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
644
    i = 0;
645
646
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
647
648
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
649
650
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
651

652
      int c = 0;
653
654
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
655
	   dofIt != sendIt->second.end();
656
	   ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
657
	sendBuffers[i][c++] = *(dofIt->first);
658
	sendBuffers[i][c++] = dofIt->second;
659

660
	sendDofs[sendIt->first].push_back(dofIt->first);
661
662
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
663
664
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
665
666
667
    }

    i = 0;
668
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
669
670
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
671
672
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
673

Thomas Witkowski's avatar
Thomas Witkowski committed
674
675
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
676
677
678
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
679
    MPI::Request::Waitall(requestCounter, request);
680
681

    
682
    // === Delete send buffers. ===
683

Thomas Witkowski's avatar
Thomas Witkowski committed
684
685
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
686
687


688
689
690
    // === Change dof indices for rank partition. ===

    mapLocalGlobalDOFs.clear();
691

692
    // === Change dof indices at boundary from other ranks. ===
693

Thomas Witkowski's avatar
Thomas Witkowski committed
694
695
696
697
698
699
700
701
    // 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
    // a and b in boundaryDOFs. Then we have to change index a to b and b to c. When
    // 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;
Thomas Witkowski's avatar
Thomas Witkowski committed
702
703
    for (DofToRank::iterator dofIt = boundaryDOFs.begin(); dofIt != boundaryDOFs.end();
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
704
705
      dofChanged[dofIt->first] = false;

706
    i = 0;
707
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
708
709
710
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

711
      for (int j = 0; j < recvIt->second; j++) {
712

713
714
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
715

716
	bool found = false;
717

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

Thomas Witkowski's avatar
Thomas Witkowski committed
720
721
	for (DofToRank::iterator dofIt = boundaryDOFs.begin(); 
	     dofIt != boundaryDOFs.end(); ++dofIt) {
722

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

726
	    recvDofs[recvIt->first].push_back(dofIt->first);
727
	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
728
	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
729

730
	    found = true;
731
732
733
	    break;
	  }
	}
734

Thomas Witkowski's avatar
Thomas Witkowski committed
735
	TEST_EXIT_DBG(found)("Should not happen!\n");
736
737
738
739
      }

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

741
742
743
    
    // === Create now the local to global index, and vice verse, mappings.     ===

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
744
745
746
    for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin();
	 dofIt != rankDofsNewLocalIndex.end(); ++dofIt) {
      DegreeOfFreedom localDof = dofIt->second;
747
      DegreeOfFreedom globalDof = rankDofsNewGlobalIndex[dofIt->first];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
748
749
750

      *const_cast<DegreeOfFreedom*>(dofIt->first) = localDof;
      mapLocalGlobalDOFs[localDof] = globalDof;
Thomas Witkowski's avatar
Thomas Witkowski committed
751
752
    }

753
754
755
756
    mapLocalToDofIndex.clear();
    for (DofIndexMap::iterator dofIt = rankOwnedDofsNewLocalIndex.begin();
	 dofIt != rankOwnedDofsNewLocalIndex.end(); ++dofIt)
      mapLocalToDofIndex[dofIt->second] = *(dofIt->first);
757
758
759
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
760
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs)
761
  {
762
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
763

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

766
    // === Get all DOFs in ranks partition. ===
767

Thomas Witkowski's avatar
Thomas Witkowski committed
768
    ElementDofIterator elDofIt(feSpace);
769
    DofSet rankDOFSet;
770
771
772
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
773
      Element *element = elInfo->getElement();     
Thomas Witkowski's avatar
Thomas Witkowski committed
774
775
      elDofIt.reset(element);
      do {
776
	rankDOFSet.insert(elDofIt.getDofPtr());
Thomas Witkowski's avatar
Thomas Witkowski committed
777
      } while(elDofIt.next());
778
779
780
781

      elInfo = stack.traverseNext(elInfo);
    }

782
783
784
785
786
787
788
    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
789
    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
790
    // === rankDOFs to boundaryDOFs.                                               ===
791

792
    DofToRank newBoundaryDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
793
794
    RankToDofContainer sendNewDofs;
    RankToDofContainer recvNewDofs;
795

Thomas Witkowski's avatar
Thomas Witkowski committed
796
797
798
    for (RankToBoundMap::iterator it =  myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end(); ++it) {

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

802
	const DegreeOfFreedom *dof1, *dof2;
Thomas Witkowski's avatar
Thomas Witkowski committed
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
	  dof1 = boundIt->rankObject.el->getDOF(1);
	  dof2 = boundIt->rankObject.el->getDOF(2);
	  break;
	case 1:
	  dof1 = boundIt->rankObject.el->getDOF(0);
	  dof2 = boundIt->rankObject.el->getDOF(2);
	  break;
	case 2:
	  dof1 = boundIt->rankObject.el->getDOF(0);
	  dof2 = boundIt->rankObject.el->getDOF(1);
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
821
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
822
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
823
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
824
825
826
827
	  ("Should never happen!\n");

	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
Thomas Witkowski's avatar
Thomas Witkowski committed
828

829
830
831
832
833
834
835
	DofContainer &dofsToSend = sendNewDofs[it->first];

  	if (find(dofsToSend.begin(), dofsToSend.end(), dof1) == dofsToSend.end())
 	  dofsToSend.push_back(dof1);
  	if (find(dofsToSend.begin(), dofsToSend.end(), dof2) == dofsToSend.end())
 	  dofsToSend.push_back(dof2);

Thomas Witkowski's avatar
Thomas Witkowski committed
836
837
838
	DofContainer boundDOFs;
	addAllVertexDOFs(boundIt->rankObject.el, 
			 boundIt->rankObject.ithObjAtBoundary,
Thomas Witkowski's avatar
Thomas Witkowski committed
839
840
841
842
			 boundDOFs);	
  	addAllEdgeDOFs(boundIt->rankObject.el, 
  		       boundIt->rankObject.ithObjAtBoundary,
  		       boundDOFs);
843
844
845

	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
846
	  dofsToSend.push_back(boundDOFs[i]);
847
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
848
	
849
850
851
      }
    }    

Thomas Witkowski's avatar
Thomas Witkowski committed
852
853
854
    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end(); ++it) {

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

858
	const DegreeOfFreedom *dof1, *dof2;
859
860
861
862
863
864
865
866
867
868
869

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
	  dof1 = boundIt->rankObject.el->getDOF(1);
	  dof2 = boundIt->rankObject.el->getDOF(2);
	  break;
	case 1:
	  dof1 = boundIt->rankObject.el->getDOF(0);
	  dof2 = boundIt->rankObject.el->getDOF(2);
	  break;
	case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
870
871
	  dof1 = boundIt->rankObject.el->getDOF(1);
	  dof2 = boundIt->rankObject.el->getDOF(0);
872
873
874
875
876
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
877
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
878
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
879
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
880
881
	  ("Should never happen!\n");

882
883
884
885
886
887
888
	DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), dof1);
	if (eraseIt != rankDOFs.end())
	  rankDOFs.erase(eraseIt);
	eraseIt = find(rankDOFs.begin(), rankDOFs.end(), dof2);
	if (eraseIt != rankDOFs.end())
	  rankDOFs.erase(eraseIt);

889
890
	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
Thomas Witkowski's avatar
Thomas Witkowski committed
891

892
893
894
895
896
897
898
	DofContainer &dofsToRecv = recvNewDofs[it->first];
  	if (find(dofsToRecv.begin(), dofsToRecv.end(), dof1) == dofsToRecv.end())
 	  dofsToRecv.push_back(dof1);
  	if (find(dofsToRecv.begin(), dofsToRecv.end(), dof2) == dofsToRecv.end())
 	  dofsToRecv.push_back(dof2);

	DofContainer boundDOFs;	
Thomas Witkowski's avatar
Thomas Witkowski committed
899
	addAllEdgeDOFs(boundIt->rankObject.el, 
Thomas Witkowski's avatar
Thomas Witkowski committed
900
901
 		       boundIt->rankObject.ithObjAtBoundary,
 		       boundDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
902
903
904
905
	addAllVertexDOFs(boundIt->rankObject.el, 
			 boundIt->rankObject.ithObjAtBoundary,
			 boundDOFs);

Thomas Witkowski's avatar
Thomas Witkowski committed
906
	for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) {
907
908
909
910
911
912
913
	  TEST_EXIT_DBG(find(rankDOFs.begin(), rankDOFs.end(), boundDOFs[i]) != rankDOFs.end())
	    ("Should never happen!\n");

	  eraseIt = find(rankDOFs.begin(), rankDOFs.end(), boundDOFs[i]);
	  if (eraseIt != rankDOFs.end())
	    rankDOFs.erase(eraseIt);

914
	  newBoundaryDOFs[boundDOFs[i]] = it->first;
915
	  dofsToRecv.push_back(boundDOFs[i]);
916
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
917
918
919
      }
    }

920
921
922
923
924
    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
925
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
926
927
928
929
930
    rstart -= nRankDOFs;


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

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


934
935
936
    // Do not change the indices now, but create a new indexing a store it here.
    DofIndexMap rankDofsNewLocalIndex;
    isRankDof.clear();
937
    int i = 0;
938
939
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
940

941
942
943
      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
944
      isRankDof[i] = true;
945
      i++;
946
947
    }

948
    // Stores for all rank owned dofs a new global index.
949
    DofIndexMap rankDofsNewGlobalIndex;
950
951
952
953
954
955
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
956
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
957
958
959
960
961
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
      i++;
    }


962
963
    // === Send new DOF indices. ===

964
965
966
    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

Thomas Witkowski's avatar
Thomas Witkowski committed
967
968
969
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

970
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
971
    for (RankToDofContainer::iterator sendIt = sendNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
972
	 sendIt != sendNewDofs.end(); ++sendIt, i++) {
973
974
975
      int nSendDofs = sendIt->second.size();
      sendBuffers[i] = new int[nSendDofs];
      int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
976
      for (DofContainer::iterator dofIt = sendIt->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
977
	   dofIt != sendIt->second.end(); ++dofIt)
978
	sendBuffers[i][c++] = rankDofsNewGlobalIndex[*dofIt];
979

Thomas Witkowski's avatar
Thomas Witkowski committed
980
981
      request[requestCounter++] = 
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
982
983
984
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
985
    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
986
	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
987
988
989
      int nRecvDofs = recvIt->second.size();
      recvBuffers[i] = new int[nRecvDofs];
	
Thomas Witkowski's avatar
Thomas Witkowski committed
990
991
      request[requestCounter++] = 
	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
992
993
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
994
995
    MPI::Request::Waitall(requestCounter, request);

996
997
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
998
999

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1000
    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1001
	 recvIt != recvNewDofs.end(); ++recvIt) {      </