ParallelDomainProblem.cc 43.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#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"
11
12
#include "DOFMatrix.h"
#include "DOFVector.h"
13
#include "VtkWriter.h"
14
#include "ElementDofIterator.h"
15
16

#include "petscksp.h"
17
18
19

namespace AMDiS {

20
21
22
23
24
25
26
27
28
29
30
31
32
  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)
33
34
    : iterationIF(iIF),
      timeIF(tIF),
35
36
      feSpace(fe),
      mesh(fe->getMesh()),
37
      refinementManager(refineManager),
38
      initialPartitionMesh(true),
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
39
40
      nRankDOFs(0),
      rstart(0)
41
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
42
43
44
45
46
47
48
    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");

49
50
51
52
53
54
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
  }

55
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
56
57
58
59
  {
    if (mpiSize <= 1)
      return;

60
61
62
63
64
    // 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();

65
66
67
68
69
70
71
    // 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
72
73
74
75
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
    DbgCreateElementMap(elMap);
#endif
76

77
    // === Create new global and local DOF numbering. ===
78

79
80
81
    // Set of all DOFs of the rank.
    std::vector<const DegreeOfFreedom*> rankDOFs;
    // Number of DOFs in ranks partition that are owned by the rank.
82
    int nRankDOFs = 0;
83
    // Number of all DOFs in the macro mesh.
84
    int nOverallDOFs = 0;
85
86

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

Thomas Witkowski's avatar
Thomas Witkowski committed
88
89
    // === Create interior boundary information ===

90
    createInteriorBoundaryInfo(rankDOFs, boundaryDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
91

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
92
93
94
95
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

Thomas Witkowski's avatar
Thomas Witkowski committed
96
#if (DEBUG != 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
97
    DbgTestElementMap(elMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
98
99
    DbgTestInteriorBoundary();
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
100

101
    // === Reset all DOFAdmins of the mesh. ===
102
103
104

    int nAdmins = mesh->getNumberOfDOFAdmin();
    for (int i = 0; i < nAdmins; i++) {
105
106
107
108
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));

      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
109
      for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
110
111
 	admin.setDOFFree(j, false);

112
      admin.setUsedSize(mapLocalGlobalDOFs.size());
113
114
      admin.setUsedCount(mapLocalGlobalDOFs.size());
      admin.setFirstHole(mapLocalGlobalDOFs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
117

118
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
119

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
126
127
128
      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
129
#if (DEBUG != 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
130
    DbgTestCommonDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
131
#endif
132

133
    // === Create petsc matrix. ===
134

135
136
    int ierr;
    ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
137
    ierr = MatSetSizes(petscMatrix, nRankDOFs, nRankDOFs, nOverallDOFs, nOverallDOFs);
138
139
140
    ierr = MatSetType(petscMatrix, MATAIJ);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
141
    ierr = VecSetSizes(petscRhsVec, nRankDOFs, nOverallDOFs);
142
    ierr = VecSetType(petscRhsVec, VECMPI);
143
144

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

149
150
151
152

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

154

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

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

200
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
201
	}
202

203
204
205
206

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

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

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
214
    }    
215
216
  }

217

218
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
219
220
221
222
223
224
225
226
227
228
  {
    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
    KSPGetPC(ksp, &pc);
    PCSetType(pc, PCJACOBI);
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
229
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
230
231
232
233
234
235
236
    KSPSolve(ksp, petscRhsVec, petscSolVec);

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

    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    int counter = 0;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
237
238
239
240
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
      *dofIt = vecPointer[mapGlobalLocalDOFs[rstart + counter]];
      counter++;
    }
241
242
243

    VecRestoreArray(petscSolVec, &vecPointer);

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
244
#if 0
245
246
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
247
248
249

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
250
251
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
252
253
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); 
254
	 ++sendIt, i++) {
255
256
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
257

258
259
      for (int j = 0; j < nSendDOFs; j++)
	sendBuffers[i][j] = (*vec)[(sendIt->second)[j][0]];
260

Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
263
264
265
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); 
268
	 ++recvIt, i++) {
269
270
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
274
275
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
278

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

279
280
    
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
282
283
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
284
285
286
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++) {
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
      }
287
288
289
290

      delete [] recvBuffers[i];
    }
    
291
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
292
      delete [] sendBuffers[i];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
293
#endif
294
295
  }

296

297
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
  {
    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
318
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
319
320
321
322
323
324
325
326
327
328
329
330
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

331

332
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
333
334
335
336
337
338
339
340
341
342
343
344
345
  {
    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);
  }

346

Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs,
						      DofToRank& boundaryDOFs)
349
  {
350
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353

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

354
    TraverseStack stack;
355
356
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
357
358
359
360
361
362
363
364
365
366
367
368
    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));
369

370
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
371
372
	    // We have found an element that is at an interior boundary. 

373
374
	    // === 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. ===
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399

	    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;

400
	    // === And add the part of the interior boundary. ===
401
402

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
403
404
405
406
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

407
408
409
410
411
412
	    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;
413
414
415
416
417
418
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
421


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
427
428
429
430
    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
433
434
435
      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
      for (int i = 0; i < nSendInt; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438
	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
439
440
      request[requestCounter++] =
	mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
441
442
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
443
444
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
445
446
      int nRecvInt = rankIt->second.size();
      int *buffer = new int[nRecvInt];
Thomas Witkowski's avatar
Thomas Witkowski committed
447
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
448

Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
      request[requestCounter++] = 
	mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
453
454
455

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

Thomas Witkowski's avatar
Thomas Witkowski committed
456
457

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
460
461
462
463
464
465
466
467
468
469
470
471
472
473

      // === 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
474
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
475
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
478
479
480
481
482
483
484
485
486
487
488

	  // 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];
489
490
491
  }


492
  void ParallelDomainBase::removeMacroElements()
493
494
495
496
497
498
499
500
501
502
503
504
  {
    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
505
    mesh->removeMacroElements(macrosToRemove, feSpace);
506
507
508
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
509
510
511
512
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
						      DofToRank& boundaryDOFs,
						      int& nRankDOFs, 
						      int& nOverallDOFs)
513
  {
514
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
515
516

    // === Get rank information about DOFs. ===
517
518
519

    // 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
520
    DofContainer rankAllDofs;
521

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
522
    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDOFs);
523
524
525
526

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

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

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
529
    //    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
530
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
531
532
    rstart -= nRankDOFs;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
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
558
559
560
561
562
563

    typedef std::map<const DegreeOfFreedom*, DegreeOfFreedom> DofIndexMap;
    DofIndexMap rankDofsNewLocalIndex, rankDofsNewGlobalIndex;
    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++;
    }

    if (mpiRank == 0) {
      for (DofContainer::iterator dofIt = rankAllDofs.begin();
	   dofIt != rankAllDofs.end(); ++dofIt) {
	std::cout << "OLD DOF = " << **dofIt << "   NEW L DOF = " << rankDofsNewLocalIndex[*dofIt];
	if (rankDofsNewGlobalIndex.find(*dofIt) != rankDofsNewGlobalIndex.end())
	  std::cout << "    NEW G DOF = " << rankDofsNewGlobalIndex[*dofIt];
	std::cout << std::endl;
      }
    }

    exit(0);


  
564
565
566
567
568
569
    // === 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;
570

571
572
573
574
    // 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;
575

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

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

	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
	  }
592
593
594
595
	}
      } else {
	// 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.
596
597
598
599
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
600
601
602
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
603

604
    // === Send and receive the dof indices at boundary. ===
605

Thomas Witkowski's avatar
Thomas Witkowski committed
606
    std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size());
607
608
609

    sendDofs.clear();
    recvDofs.clear();
610
    
Thomas Witkowski's avatar
Thomas Witkowski committed
611
612
613
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

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

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

630
	sendDofs[sendIt->first].push_back(dofIt->first);
631
632
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
633
634
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
635
636
637
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
646
647
648
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
649
    MPI::Request::Waitall(requestCounter, request);
650
651

    
652
    // === Delete send buffers. ===
653
654

    i = 0;
655
656
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
657
658
659
660
661
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) 
      delete [] sendBuffers[i];


662
663
664
665
    // === Change dof indices for rank partition. ===

    mapLocalGlobalDOFs.clear();
    mapGlobalLocalDOFs.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
666
    isRankDof.clear();
667

668
    // === Change dof indices at boundary from other ranks. ===
669

Thomas Witkowski's avatar
Thomas Witkowski committed
670
671
672
673
674
675
676
677
    // 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
678
679
    for (DofToRank::iterator dofIt = boundaryDOFs.begin(); dofIt != boundaryDOFs.end();
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
680
681
      dofChanged[dofIt->first] = false;

682
    i = 0;
683
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
684
685
686
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

687
      for (int j = 0; j < recvIt->second; j++) {
688

689
690
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
691

692
	bool found = false;
693

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

Thomas Witkowski's avatar
Thomas Witkowski committed
696
697
	for (DofToRank::iterator dofIt = boundaryDOFs.begin(); 
	     dofIt != boundaryDOFs.end(); ++dofIt) {
698

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

702
	    recvDofs[recvIt->first].push_back(dofIt->first);
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
703
704
	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
705

706
	    found = true;
707
708
709
	    break;
	  }
	}
710

Thomas Witkowski's avatar
Thomas Witkowski committed
711
	TEST_EXIT_DBG(found)("Should not happen!\n");
712
713
714
715
      }

      delete [] recvBuffers[i];
    }
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
716
717
718
719
720
721
722
723
724
725

    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;
    }
726
727
728
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
729
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs)
730
  {
731
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
732

733
    std::set<const DegreeOfFreedom*> rankDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
734
    DofToRank newBoundaryDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
735

Thomas Witkowski's avatar
Thomas Witkowski committed
736

737
    // === Get all DOFs in ranks partition. ===
738

Thomas Witkowski's avatar
Thomas Witkowski committed
739
    ElementDofIterator elDofIt(feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
740

741
742
743
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
744
      Element *element = elInfo->getElement();     
Thomas Witkowski's avatar
Thomas Witkowski committed
745
746
747
748
      elDofIt.reset(element);
      do {
	rankDOFs.insert(elDofIt.getDofPtr());
      } while(elDofIt.next());
749
750
751
752

      elInfo = stack.traverseNext(elInfo);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
756
757
    RankToDofContainer sendNewDofs;
    RankToDofContainer recvNewDofs;
758

Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
761
    for (RankToBoundMap::iterator it =  myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end(); ++it) {

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

Thomas Witkowski's avatar
Thomas Witkowski committed
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
	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
785
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
786
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
787
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
788
789
790
791
	  ("Should never happen!\n");

	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
Thomas Witkowski's avatar
Thomas Witkowski committed
792
793
794
795
796
797
798
799

  	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
800
	DofContainer boundDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
801
	
Thomas Witkowski's avatar
Thomas Witkowski committed
802
803
	addAllVertexDOFs(boundIt->rankObject.el, 
			 boundIt->rankObject.ithObjAtBoundary,
Thomas Witkowski's avatar
Thomas Witkowski committed
804
805
806
807
			 boundDOFs);	
  	addAllEdgeDOFs(boundIt->rankObject.el, 
  		       boundIt->rankObject.ithObjAtBoundary,
  		       boundDOFs);
808
809
810

	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
811
	  sendNewDofs[it->first].push_back(boundDOFs[i]);
812
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
813
	
814
815
816
      }
    }    

Thomas Witkowski's avatar
Thomas Witkowski committed
817
818
819
    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end(); ++it) {

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

823
824
825
826
827
828
829
830
831
832
833
834
835
	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
836
837
	  dof1 = boundIt->rankObject.el->getDOF(1);
	  dof2 = boundIt->rankObject.el->getDOF(0);
838
839
840
841
842
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
843
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
844
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
845
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
846
847
848
849
850
851
	  ("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
852

Thomas Witkowski's avatar
Thomas Witkowski committed
853
854
855
856
857
858
859
  	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
860
	DofContainer boundDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
861
	
Thomas Witkowski's avatar
Thomas Witkowski committed
862
	addAllEdgeDOFs(boundIt->rankObject.el, 
Thomas Witkowski's avatar
Thomas Witkowski committed
863
864
 		       boundIt->rankObject.ithObjAtBoundary,
 		       boundDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
865
866
867
868
	addAllVertexDOFs(boundIt->rankObject.el, 
			 boundIt->rankObject.ithObjAtBoundary,
			 boundDOFs);

Thomas Witkowski's avatar
Thomas Witkowski committed
869
	for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) {
870
871
	  rankDOFs.erase(boundDOFs[i]);
	  newBoundaryDOFs[boundDOFs[i]] = it->first;
872
	  recvNewDofs[it->first].push_back(boundDOFs[i]);
873
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
874
875
876
      }
    }

877
878
879
880
881
    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
882
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
883
884
885
886
887
    rstart -= nRankDOFs;


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

Thomas Witkowski's avatar
Thomas Witkowski committed
888
    mpiComm.Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
889
890
891
892
893
894


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

    mapLocalGlobalDOFs.clear();
    mapGlobalLocalDOFs.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
895
    isRankDof.clear();
896
897
898

    int i = 0;
    for (std::set<const DegreeOfFreedom*>::iterator dofIt = rankDOFs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
899
900
	 dofIt != rankDOFs.end(); ++dofIt, i++) {
      *(const_cast<DegreeOfFreedom*>(*dofIt)) = i;
901
902
      mapLocalGlobalDOFs[i] = rstart + i;
      mapGlobalLocalDOFs[rstart + i] = i;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
903
      isRankDof[i] = true;
904
905
906
907
    }

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

908
909
910
    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

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

914
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
915
    for (RankToDofContainer::iterator sendIt = sendNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
916
	 sendIt != sendNewDofs.end(); ++sendIt, i++) {
917
918
919
      int nSendDofs = sendIt->second.size();
      sendBuffers[i] = new int[nSendDofs];
      int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
920
      for (DofContainer::iterator dofIt = sendIt->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
921
	   dofIt != sendIt->second.end(); ++dofIt)
922
	sendBuffers[i][c++] = (*dofIt)[0] + rstart;
923

Thomas Witkowski's avatar
Thomas Witkowski committed
924
925
      request[requestCounter++] = 
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
926
927
928
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
929
    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
930
	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
931
932
933
      int nRecvDofs = recvIt->second.size();
      recvBuffers[i] = new int[nRecvDofs];
	
Thomas Witkowski's avatar
Thomas Witkowski committed
934
935
      request[requestCounter++] = 
	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
936
937
    }

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

940
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
941
    for (RankToDofContainer::iterator sendIt = sendNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
942
	 sendIt != sendNewDofs.end(); ++sendIt)
943
944
945
946
      delete [] sendBuffers[i++];

    i = 0;
    int newDofIndex = nRankDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
947
    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
948
	 recvIt != recvNewDofs.end(); ++recvIt) {      
949
      int j = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
950
      for (DofContainer::iterator dofIt = recvIt->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
951
	   dofIt != recvIt->second.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
952
	(*const_cast<DegreeOfFreedom*>(*dofIt)) = newDofIndex;
953
954
	mapLocalGlobalDOFs[newDofIndex] = recvBuffers[i][j];
	mapGlobalLocalDOFs[recvBuffers[i][j]] = newDofIndex;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
955
	isRankDof[newDofIndex] = false;
956
957
958
959
960
	newDofIndex++;
	j++;
      }

      delete [] recvBuffers[i++];
961
962
    }

963
964
965
966
    // === Update list of dofs that must be communicated for solution exchange. ===

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

968
969
970
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
971
972
  void ParallelDomainBase::addAllVertexDOFs(Element *el, int ithEdge, 
					    DofContainer& dofs)
973
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
974
    FUNCNAME("ParallelDomainBase::addAllVertexDOFs()");
975

Thomas Witkowski's avatar
Thomas Witkowski committed
976
977
978
    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
979
	addAllVertexDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
980
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
981
	addAllVertexDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
982
983
984
985
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
986
	addAllVertexDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
987
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
988
	addAllVertexDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
989
990
991
992
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
993
	addAllVertexDOFs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
994
	dofs.push_back(el->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
995
	addAllVertexDOFs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
996
997
998
999
1000
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
For faster browsing, not all history is shown. View entire blame