ParallelDomainProblem.cc 43.7 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
27
28
29
30
31
32
33
34
35
36
37
  PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
  {
    if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
      std::cout << "  Iteration " << iter << ": " << rnorm << "\n";

    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
    // === Reset all DOFAdmins of the mesh. ===
107
108
109

    int nAdmins = mesh->getNumberOfDOFAdmin();
    for (int i = 0; i < nAdmins; i++) {
110
111
112
113
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));

      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
114
      for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
115
116
 	admin.setDOFFree(j, false);

117
      admin.setUsedSize(mapLocalGlobalDOFs.size());
118
119
      admin.setUsedCount(mapLocalGlobalDOFs.size());
      admin.setFirstHole(mapLocalGlobalDOFs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
    }

122
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
123

Thomas Witkowski's avatar
Thomas Witkowski committed
124
125
126
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

Thomas Witkowski's avatar
Thomas Witkowski committed
127
128
    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);
129

Thomas Witkowski's avatar
Thomas Witkowski committed
130
131
132
      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
    }

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
    // === Create petsc matrix. ===
138

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

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

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

153
154
155
156

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

158

159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
  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)
183
  {
184
185
186
187
188
189
190
191
192
193
194
    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;

195
196
197
198
    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)
199
	if (value(*icursor) != 0.0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
200
201
 	  int r = mapLocalGlobalDOFs[row(*icursor)];
 	  int c = mapLocalGlobalDOFs[col(*icursor)];  
202
	  double v = value(*icursor);
203

204
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
205
	}
206

207
208
209
210

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

211
    
212
213
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
214
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
215
216
217
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
218
    }    
219
220
  }

221

222
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
223
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

226
227
228
229
230
231
    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
    KSPGetPC(ksp, &pc);
Thomas Witkowski's avatar
Thomas Witkowski committed
232
    PCSetType(pc, PCNONE);
233
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
Thomas Witkowski's avatar
Thomas Witkowski committed
234
235
    //    KSPSetType(ksp, KSPBCGS);
    KSPSetType(ksp, KSPCG);
236
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
237
238
    KSPSolve(ksp, petscRhsVec, petscSolVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
239
240
241
242
243
244
#if (DEBUG != 0)
    int size = 0;
    VecGetLocalSize(petscSolVec, &size);
    TEST_EXIT(size == nRankDOFs)("Vector and rank DOFs does not fit together!\n");
#endif

245
246
    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
247
    vec->set(1.0);
248

Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
    for (int i = 0; i < nRankDOFs; i++)
      (*vec)[mapGlobalLocalDOFs[i]] = vecPointer[i];
251
252
253

    VecRestoreArray(petscSolVec, &vecPointer);

Thomas Witkowski's avatar
Thomas Witkowski committed
254
#if 1
255
256
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
257
258
259

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
260
261
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
262
263
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); 
264
	 ++sendIt, i++) {
265
266
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
267

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

Thomas Witkowski's avatar
Thomas Witkowski committed
271
272
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
273
274
275
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); 
278
	 ++recvIt, i++) {
279
280
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
281

Thomas Witkowski's avatar
Thomas Witkowski committed
282
283
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
284
285
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
286
287
288

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

289
290
    
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
291
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
292
293
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
294
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
295
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
296
297
298
299

      delete [] recvBuffers[i];
    }
    
300
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
301
      delete [] sendBuffers[i];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
302
#endif
303
304
  }

305

306
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
  {
    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
327
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
328
329
330
331
332
333
334
335
336
337
338
339
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

340

341
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
342
343
344
345
346
347
348
349
350
351
352
353
354
  {
    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);
  }

355

Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs,
						      DofToRank& boundaryDOFs)
358
  {
359
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
362

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

363
    TraverseStack stack;
364
365
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
366
367
368
369
370
371
372
373
374
375
376
377
    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));
378

379
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
380
381
	    // We have found an element that is at an interior boundary. 

382
383
	    // === 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. ===
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

	    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;

409
	    // === And add the part of the interior boundary. ===
410
411

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
414
415
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

416
417
418
419
420
421
	    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;
422
423
424
425
426
427
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
428
429
430


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438
439
    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
442
443
444
      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
      for (int i = 0; i < nSendInt; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
445
446
447
	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
448
449
      request[requestCounter++] =
	mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
450
451
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
454
455
      int nRecvInt = rankIt->second.size();
      int *buffer = new int[nRecvInt];
Thomas Witkowski's avatar
Thomas Witkowski committed
456
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
457

Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
      request[requestCounter++] = 
	mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
460
461
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
462
463
464

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

Thomas Witkowski's avatar
Thomas Witkowski committed
465
466

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
467
468
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
469
470
471
472
473
474
475
476
477
478
479
480
481
482

      // === 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
483
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
484
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
487
488
489
490
491
492
493
494
495
496
497

	  // 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];
498
499
500
  }


501
  void ParallelDomainBase::removeMacroElements()
502
503
504
505
506
507
508
509
510
511
512
513
  {
    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
514
    mesh->removeMacroElements(macrosToRemove, feSpace);
515
516
517
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
518
519
520
521
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
						      DofToRank& boundaryDOFs,
						      int& nRankDOFs, 
						      int& nOverallDOFs)
522
  {
523
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
524
525

    // === Get rank information about DOFs. ===
526
527
528

    // 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
529
    DofContainer rankAllDofs;
530

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
531
    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDOFs);
532
533
534
535

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

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

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
538
    //    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
539
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
540
541
    rstart -= nRankDOFs;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
542

Thomas Witkowski's avatar
Thomas Witkowski committed
543
544
545
    rankDofsNewLocalIndex.clear();
    rankDofsNewGlobalIndex.clear();

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
546
547
548
549
550
551
552
553
554
555
556
557
558
559
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
      rankDofsNewLocalIndex[*dofIt] = i++;
      isRankDof[i] = true;
    }

    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
      i++;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
560
 
561
562
563
564
565
566
    // === 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;
567

568
569
570
571
    // 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;
572

Thomas Witkowski's avatar
Thomas Witkowski committed
573
    for (DofToRank::iterator it = boundaryDOFs.begin(); it != boundaryDOFs.end(); ++it) {
574
575
576
577
578
579
580
581

      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
582
583
584
585
586
587
	  if (*itRanks != mpiRank) {
	    TEST_EXIT_DBG(rankDofsNewGlobalIndex.count(it->first) == 1)
	      ("DOF Key not found!\n");

	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
	  }
588
589
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
590
591
	// 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.
592
593
594
595
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
596
597
598
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
599

600
    // === Send and receive the dof indices at boundary. ===
601

Thomas Witkowski's avatar
Thomas Witkowski committed
602
    std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size());
603
604
605

    sendDofs.clear();
    recvDofs.clear();
606
    
Thomas Witkowski's avatar
Thomas Witkowski committed
607
608
609
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
610
    i = 0;
611
612
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
613
614
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
615
616
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
617

618
      int c = 0;
619
620
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
621
	   dofIt != sendIt->second.end();
622
	   ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
623
	sendBuffers[i][c++] = *(dofIt->first);
624
	sendBuffers[i][c++] = dofIt->second;
625

626
	sendDofs[sendIt->first].push_back(dofIt->first);
627
628
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
631
632
633
    }

    i = 0;
634
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
635
636
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
637
638
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
639

Thomas Witkowski's avatar
Thomas Witkowski committed
640
641
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
642
643
644
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
645
    MPI::Request::Waitall(requestCounter, request);
646
647

    
648
    // === Delete send buffers. ===
649

Thomas Witkowski's avatar
Thomas Witkowski committed
650
651
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
652
653


654
655
656
657
    // === Change dof indices for rank partition. ===

    mapLocalGlobalDOFs.clear();
    mapGlobalLocalDOFs.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
658
    isRankDof.clear();
659

660
    // === Change dof indices at boundary from other ranks. ===
661

Thomas Witkowski's avatar
Thomas Witkowski committed
662
663
664
665
666
667
668
669
    // 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
670
671
    for (DofToRank::iterator dofIt = boundaryDOFs.begin(); dofIt != boundaryDOFs.end();
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
672
673
      dofChanged[dofIt->first] = false;

674
    i = 0;
675
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
676
677
678
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

679
      for (int j = 0; j < recvIt->second; j++) {
680

681
682
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
683

684
	bool found = false;
685

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

Thomas Witkowski's avatar
Thomas Witkowski committed
688
689
	for (DofToRank::iterator dofIt = boundaryDOFs.begin(); 
	     dofIt != boundaryDOFs.end(); ++dofIt) {
690

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

694
	    recvDofs[recvIt->first].push_back(dofIt->first);
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
695
696
	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
697

698
	    found = true;
699
700
701
	    break;
	  }
	}
702

Thomas Witkowski's avatar
Thomas Witkowski committed
703
	TEST_EXIT_DBG(found)("Should not happen!\n");
704
705
706
707
      }

      delete [] recvBuffers[i];
    }
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
708
709
710
711
712
713
714
715
716
717

    for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin();
	 dofIt != rankDofsNewLocalIndex.end(); ++dofIt) {
      DegreeOfFreedom localDof = dofIt->second;
      DegreeOfFreedom globalDof = rankDofsNewGlobalIndex[dofIt->first];

      *const_cast<DegreeOfFreedom*>(dofIt->first) = localDof;
      mapLocalGlobalDOFs[localDof] = globalDof;
      mapGlobalLocalDOFs[globalDof] = localDof;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744

#if 0
    if (mpiRank == 0) {
      for (DofContainer::iterator dofIt = recvDofs[1].begin();
	   dofIt != recvDofs[1].end(); ++dofIt) {
	std::cout << "RECV " << **dofIt << std::endl;
      }
    }

    exit(0);
#endif

#if 0
    if (mpiRank == 0) {
      for (DofContainer::iterator dofIt = rankAllDofs.begin();
	   dofIt != rankAllDofs.end(); ++dofIt) {
	std::cout << "DOF = " << **dofIt << " GLOBAL index = " 
		  << rankDofsNewGlobalIndex[*dofIt] << " FROM MAP = " 
		  << mapLocalGlobalDOFs[**dofIt]
		  << "  AND BACK "
		  << mapGlobalLocalDOFs[mapLocalGlobalDOFs[**dofIt]] 
		  << std::endl;
      }
    }

    exit(0);
#endif
745
746
747
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
748
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs)
749
  {
750
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
751

752
    std::set<const DegreeOfFreedom*> rankDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
753
    DofToRank newBoundaryDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
754

Thomas Witkowski's avatar
Thomas Witkowski committed
755

756
    // === Get all DOFs in ranks partition. ===
757

Thomas Witkowski's avatar
Thomas Witkowski committed
758
    ElementDofIterator elDofIt(feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
759

760
761
762
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
763
      Element *element = elInfo->getElement();     
Thomas Witkowski's avatar
Thomas Witkowski committed
764
765
766
767
      elDofIt.reset(element);
      do {
	rankDOFs.insert(elDofIt.getDofPtr());
      } while(elDofIt.next());
768
769
770
771

      elInfo = stack.traverseNext(elInfo);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
    RankToDofContainer sendNewDofs;
    RankToDofContainer recvNewDofs;
777

Thomas Witkowski's avatar
Thomas Witkowski committed
778
779
780
    for (RankToBoundMap::iterator it =  myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end(); ++it) {

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

Thomas Witkowski's avatar
Thomas Witkowski committed
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
	const DegreeOfFreedom *dof1 = NULL;
	const DegreeOfFreedom *dof2 = NULL;

	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
804
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
805
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
806
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
807
808
809
810
	  ("Should never happen!\n");

	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
Thomas Witkowski's avatar
Thomas Witkowski committed
811
812
813
814
815
816
817
818

  	if (find(sendNewDofs[it->first].begin(), sendNewDofs[it->first].end(), dof1) ==
  	    sendNewDofs[it->first].end())
 	  sendNewDofs[it->first].push_back(dof1);
  	if (find(sendNewDofs[it->first].begin(), sendNewDofs[it->first].end(), dof2) ==
  	    sendNewDofs[it->first].end())
 	  sendNewDofs[it->first].push_back(dof2);
		
Thomas Witkowski's avatar
Thomas Witkowski committed
819
	DofContainer boundDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
820
	
Thomas Witkowski's avatar
Thomas Witkowski committed
821
822
	addAllVertexDOFs(boundIt->rankObject.el, 
			 boundIt->rankObject.ithObjAtBoundary,
Thomas Witkowski's avatar
Thomas Witkowski committed
823
824
825
826
			 boundDOFs);	
  	addAllEdgeDOFs(boundIt->rankObject.el, 
  		       boundIt->rankObject.ithObjAtBoundary,
  		       boundDOFs);
827
828
829

	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
830
	  sendNewDofs[it->first].push_back(boundDOFs[i]);
831
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
832
	
833
834
835
      }
    }    

Thomas Witkowski's avatar
Thomas Witkowski committed
836
837
838
    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end(); ++it) {

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

842
843
844
845
846
847
848
849
850
851
852
853
854
	const DegreeOfFreedom *dof1 = NULL;
	const DegreeOfFreedom *dof2 = NULL;

	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
855
856
	  dof1 = boundIt->rankObject.el->getDOF(1);
	  dof2 = boundIt->rankObject.el->getDOF(0);
857
858
859
860
861
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
862
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
863
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
864
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
865
866
867
868
869
870
	  ("Should never happen!\n");

	rankDOFs.erase(dof1);
	rankDOFs.erase(dof2);
	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
Thomas Witkowski's avatar
Thomas Witkowski committed
871

Thomas Witkowski's avatar
Thomas Witkowski committed
872
873
874
875
876
877
878
  	if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof1) ==
  	    recvNewDofs[it->first].end())
 	  recvNewDofs[it->first].push_back(dof1);
  	if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof2) ==
  	    recvNewDofs[it->first].end())
 	  recvNewDofs[it->first].push_back(dof2);
	
Thomas Witkowski's avatar
Thomas Witkowski committed
879
	DofContainer boundDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
880
	
Thomas Witkowski's avatar
Thomas Witkowski committed
881
	addAllEdgeDOFs(boundIt->rankObject.el, 
Thomas Witkowski's avatar
Thomas Witkowski committed
882
883
 		       boundIt->rankObject.ithObjAtBoundary,
 		       boundDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
884
885
886
887
	addAllVertexDOFs(boundIt->rankObject.el, 
			 boundIt->rankObject.ithObjAtBoundary,
			 boundDOFs);

Thomas Witkowski's avatar
Thomas Witkowski committed
888
	for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) {
889
890
	  rankDOFs.erase(boundDOFs[i]);
	  newBoundaryDOFs[boundDOFs[i]] = it->first;
891
	  recvNewDofs[it->first].push_back(boundDOFs[i]);
892
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
893
894
895
      }
    }

896
897
898
899
900
    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
901
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
902
903
904
905
906
    rstart -= nRankDOFs;


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

Thomas Witkowski's avatar
Thomas Witkowski committed
907
    mpiComm.Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
908
909
910
911
912
913


    // === Create new local DOF index numbering. ===

    mapLocalGlobalDOFs.clear();
    mapGlobalLocalDOFs.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
914
    isRankDof.clear();
915
916
917

    int i = 0;
    for (std::set<const DegreeOfFreedom*>::iterator dofIt = rankDOFs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
918
919
	 dofIt != rankDOFs.end(); ++dofIt, i++) {
      *(const_cast<DegreeOfFreedom*>(*dofIt)) = i;
920
921
      mapLocalGlobalDOFs[i] = rstart + i;
      mapGlobalLocalDOFs[rstart + i] = i;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
922
      isRankDof[i] = true;
923
924
925
926
    }

    // === Send new DOF indices. ===

927
928
929
    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

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

933
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
934
    for (RankToDofContainer::iterator sendIt = sendNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
935
	 sendIt != sendNewDofs.end(); ++sendIt, i++) {
936
937
938
      int nSendDofs = sendIt->second.size();
      sendBuffers[i] = new int[nSendDofs];
      int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
939
      for (DofContainer::iterator dofIt = sendIt->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
940
	   dofIt != sendIt->second.end(); ++dofIt)
941
	sendBuffers[i][c++] = (*dofIt)[0] + rstart;
942

Thomas Witkowski's avatar
Thomas Witkowski committed
943
944
      request[requestCounter++] = 
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
945
946
947
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
948
    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
949
	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
950
951
952
      int nRecvDofs = recvIt->second.size();
      recvBuffers[i] = new int[nRecvDofs];
	
Thomas Witkowski's avatar
Thomas Witkowski committed
953
954
      request[requestCounter++] = 
	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
955
956
    }

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

959
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
960
    for (RankToDofContainer::iterator sendIt = sendNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
961
	 sendIt != sendNewDofs.end(); ++sendIt)
962
963
964
965
      delete [] sendBuffers[i++];

    i = 0;
    int newDofIndex = nRankDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
966
    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
967
	 recvIt != recvNewDofs.end(); ++recvIt) {      
968
      int j = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
969
      for (DofContainer::iterator dofIt = recvIt->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
970
	   dofIt != recvIt->second.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
971
	(*const_cast<DegreeOfFreedom*>(*dofIt)) = newDofIndex;
972
973
	mapLocalGlobalDOFs[newDofIndex] = recvBuffers[i][j];
	mapGlobalLocalDOFs[recvBuffers[i][j]] = newDofIndex;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
974
	isRankDof[newDofIndex] = false;
975
976
977
978
979
	newDofIndex++;
	j++;
      }

      delete [] recvBuffers[i++];
980
981
    }

982
983
984
985
    // === Update list of dofs that must be communicated for solution exchange. ===

    sendDofs = sendNewDofs;
    recvDofs = recvNewDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
986

987
988
989
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
990
991
  void ParallelDomainBase::addAllVertexDOFs(Element *el, int ithEdge, 
					    DofContainer& dofs)
992
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
993
    FUNCNAME("ParallelDomainBase::addAllVertexDOFs()");
994

Thomas Witkowski's avatar
Thomas Witkowski committed
995
996
997
    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
998
	addAllVertexDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
999
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1000
	addAllVertexDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1001
1002
1003
1004
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1005
	addAllVertexDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1006
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1007
	addAllVertexDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1008
1009
1010
1011
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1012
	addAllVertexDOFs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1013
	dofs.push_back(el->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1014
	addAllVertexDOFs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1015
1016
1017
1018
1019
1020
1021
1022
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
1023
1024
1025
1026
1027
1028
1029
1030
1031
  void ParallelDomainBase::addAllEdgeDOFs(Element *el, int ithEdge,
					  DofContainer& dofs)
  {
    FUNCNAME("ParallelDomainBase::addAllEdgeDOFs()");

    bool addThisEdge = false;

    switch (ithEdge) {
    case 0:
Thomas Witkowski's avatar
Thomas Witkowski committed
1032
      if (el->getSecondChild()) {