ParallelDomainProblem.cc 38.8 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),
39
      nRankDOFs(0)
40
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
41
42
43
44
45
46
47
    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");

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

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

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

64
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);   


72
    // === Create new global and local DOF numbering. ===
73

74
75
76
    // Set of all DOFs of the rank.
    std::vector<const DegreeOfFreedom*> rankDOFs;
    // Number of DOFs in ranks partition that are owned by the rank.
77
    int nRankDOFs = 0;
78
    // Number of all DOFs in the macro mesh.
79
    int nOverallDOFs = 0;
80
81

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

Thomas Witkowski's avatar
Thomas Witkowski committed
83
84
85

    // === Create interior boundary information ===

86
    createInteriorBoundaryInfo(rankDOFs, boundaryDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
87

Thomas Witkowski's avatar
Thomas Witkowski committed
88

89
90
    // === Remove all macro elements that are not part of the rank partition. ===

91
    removeMacroElements();
92
      
93

94
    // === Reset all DOFAdmins of the mesh. ===
95
96
97

    int nAdmins = mesh->getNumberOfDOFAdmin();
    for (int i = 0; i < nAdmins; i++) {
98
99
100
101
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));

      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
102
      for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
103
104
 	admin.setDOFFree(j, false);

105
      admin.setUsedSize(mapLocalGlobalDOFs.size());
106
107
      admin.setUsedCount(mapLocalGlobalDOFs.size());
      admin.setFirstHole(mapLocalGlobalDOFs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
108
109
    }

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
118
119
120
121
      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
122
123
124
#if (DEBUG != 0)
    testInteriorBoundary();
#endif
125

126
    // === Create petsc matrix. ===
127

128
129
    int ierr;
    ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
130
    ierr = MatSetSizes(petscMatrix, nRankDOFs, nRankDOFs, nOverallDOFs, nOverallDOFs);
131
132
133
    ierr = MatSetType(petscMatrix, MATAIJ);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
134
    ierr = VecSetSizes(petscRhsVec, nRankDOFs, nOverallDOFs);
135
    ierr = VecSetType(petscRhsVec, VECMPI);
136
137

    ierr = VecCreate(PETSC_COMM_WORLD, &petscSolVec);
138
    ierr = VecSetSizes(petscSolVec, nRankDOFs, nOverallDOFs);
139
    ierr = VecSetType(petscSolVec, VECMPI);
140
141
  }

142
143
144
145

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

147

148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
  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)
172
  {
173
174
175
176
177
178
179
180
181
182
183
    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;

184
185
186
187
    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)
188
	if (value(*icursor) != 0.0) {
189
	  int r = mapLocalGlobalDOFs[row(*icursor)];
Thomas Witkowski's avatar
Thomas Witkowski committed
190
	  int c = mapLocalGlobalDOFs[col(*icursor)];  
191
	  double v = value(*icursor);
192

193
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
194
	}
195

196
197
198
199

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

200
    
201
202
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
203
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
204
205
206
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
207
    }    
208
209
  }

210

211
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
212
213
214
215
216
217
218
219
220
221
  {
    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);
222
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
223
224
225
226
227
228
229
230
    KSPSolve(ksp, petscRhsVec, petscSolVec);

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

    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    int counter = 0;
    for (dofIt.reset(); !dofIt.end(); ++dofIt)
231
      *dofIt = vecPointer[counter++];
232
233
234

    VecRestoreArray(petscSolVec, &vecPointer);

235
236
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
237
238
239

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
240
241
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
242
243
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); 
244
	 ++sendIt, i++) {
245
246
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
247

248
249
      for (int j = 0; j < nSendDOFs; j++)
	sendBuffers[i][j] = (*vec)[(sendIt->second)[j][0]];
250

Thomas Witkowski's avatar
Thomas Witkowski committed
251
252
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
253
254
255
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
256
257
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); 
258
	 ++recvIt, i++) {
259
260
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
261

Thomas Witkowski's avatar
Thomas Witkowski committed
262
263
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
264
265
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
268

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

269
270
    
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
271
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
272
273
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
276
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++) {
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
      }
277
278
279
280

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

285

286
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
  {
    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
307
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
308
309
310
311
312
313
314
315
316
317
318
319
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

320

321
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
322
323
324
325
326
327
328
329
330
331
332
333
334
  {
    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);
  }

335

Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs,
						      DofToRank& boundaryDOFs)
338
  {
339
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
340
341
342

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

343
    TraverseStack stack;
344
345
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
346
347
348
349
350
351
352
353
354
355
356
357
    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));
358

359
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
360
361
	    // We have found an element that is at an interior boundary. 

362
363
	    // === 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. ===
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388

	    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;

389
	    // === And add the part of the interior boundary. ===
390
391

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
394
395
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

396
397
398
399
400
401
	    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;
402
403
404
405
406
407
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
410
411
412
413
414
415
416


    // === Once we have this information, we must care about their order. Eventually ===
    // === all the boundaries have to be in the same order on both ranks (because    ===
    // === each boundary is shared by exactly two ranks).                            ===

    std::vector<int*> sendBuffers;
    std::vector<int*> recvBuffers;

Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
419
420
421
422
    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = 
	   myIntBoundary.boundary.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
423
424
	 rankIt != myIntBoundary.boundary.end();
	 ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
427
      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
      for (int i = 0; i < nSendInt; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
428
429
430
	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
      request[requestCounter++] =
	mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
433
434
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
435
436
    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = 
	   otherIntBoundary.boundary.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
437
438
	 rankIt != otherIntBoundary.boundary.end();
	 ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
439
440
      int nRecvInt = rankIt->second.size();
      int *buffer = new int[nRecvInt];
Thomas Witkowski's avatar
Thomas Witkowski committed
441
442
      recvBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
443
444
      request[requestCounter++] = 
	mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
445
446
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
449

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

Thomas Witkowski's avatar
Thomas Witkowski committed
450
451

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = 
	   otherIntBoundary.boundary.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
	 rankIt != otherIntBoundary.boundary.end();
	 ++rankIt) {

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

	  // 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];
485
486
487
  }


488
  void ParallelDomainBase::removeMacroElements()
489
490
491
492
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);
    }

    mesh->removeMacroElements(macrosToRemove);
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
505
506
507
508
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
						      DofToRank& boundaryDOFs,
						      int& nRankDOFs, 
						      int& nOverallDOFs)
509
  {
510
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
511

Thomas Witkowski's avatar
Thomas Witkowski committed
512

513
    // === Get rank information about DOFs. ===
514
515
516
517
518
519
520
521
522

    // Stores to each DOF pointer the set of ranks the DOF is part of.
    std::map<const DegreeOfFreedom*, std::set<int> > partitionDOFs;

    createDOFMemberInfo(partitionDOFs, rankDOFs, boundaryDOFs);

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

523
524

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

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

   
531
532
533
534
535
536
    // === 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;
537

538
539
540
541
    // 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;
542

Thomas Witkowski's avatar
Thomas Witkowski committed
543
    for (DofToRank::iterator it = boundaryDOFs.begin();
544
545
546
547
548
549
550
	 it != boundaryDOFs.end();
	 ++it) {

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

	// search for new dof index in ranks partition for this boundary dof
551
552
	DegreeOfFreedom newDofIndex = 0;
	for (int i = 0; i < nRankDOFs; i++) {
553
554
555
556
557
558
559
560
561
562
563
	  if (rankDOFs[i] == it->first) {
	    newDofIndex = rstart + i;
	    break;
	  }
	}

	// 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) {
	  if (*itRanks != mpiRank)
564
	    sendNewDofs[*itRanks][it->first] = newDofIndex;
565
566
567
568
	}
      } 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.
569
570
571
572
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
573
574
575
      }
    }

576
    // === Send and receive the dof indices at boundary. ===
577
578
579

    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());
580
581
582

    sendDofs.clear();
    recvDofs.clear();
583
    
Thomas Witkowski's avatar
Thomas Witkowski committed
584
585
586
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

587
    int i = 0;
588
589
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
590
591
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
592
593
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
594

595
      int c = 0;
596
597
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
598
	   dofIt != sendIt->second.end();
599
600
601
	   ++dofIt) {
	sendBuffers[i][c++] = (dofIt->first)[0];
	sendBuffers[i][c++] = dofIt->second;
602

603
	sendDofs[sendIt->first].push_back(dofIt->first);
604
605
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
606
607
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
608
609
610
    }

    i = 0;
611
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
612
613
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
614
615
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
616

Thomas Witkowski's avatar
Thomas Witkowski committed
617
618
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
619
620
621
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
622
    MPI::Request::Waitall(requestCounter, request);
623
624

    
625
    // === Delete send buffers. ===
626
627

    i = 0;
628
629
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
630
631
632
633
634
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) 
      delete [] sendBuffers[i];


635
636
637
638
639
    // === Change dof indices for rank partition. ===

    mapLocalGlobalDOFs.clear();
    mapGlobalLocalDOFs.clear();
    isRankDOF.clear();
640

641
    for (int i = 0; i < nRankDOFs; i++) {
642
643
644
645
646
647
648
      const_cast<DegreeOfFreedom*>(rankDOFs[i])[0] = i; 
      mapLocalGlobalDOFs[i] = rstart + i;
      mapGlobalLocalDOFs[rstart + i] = i;
      isRankDOF[i] = true;
    }


649
    // === Change dof indices at boundary from other ranks. ===
650

Thomas Witkowski's avatar
Thomas Witkowski committed
651
652
653
654
655
656
657
658
659
660
661
662
663
664
    // 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;
    for (DofToRank::iterator dofIt = boundaryDOFs.begin();
	 dofIt != boundaryDOFs.end();
	 ++dofIt) {
      dofChanged[dofIt->first] = false;
    }

665
    i = 0;
666
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
667
668
669
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

670
      for (int j = 0; j < recvIt->second; j++) {
671

672
673
674
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
	DegreeOfFreedom newLocalDof = mapLocalGlobalDOFs.size();
675

676
	bool found = false;
677

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

	for (DofToRank::iterator dofIt = boundaryDOFs.begin();
681
682
	     dofIt != boundaryDOFs.end();
	     ++dofIt) {
683

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

687
	    recvDofs[recvIt->first].push_back(dofIt->first);
Thomas Witkowski's avatar
Thomas Witkowski committed
688
689
	    *(const_cast<DegreeOfFreedom*>(dofIt->first)) = newLocalDof;

690
691
	    mapLocalGlobalDOFs[newLocalDof] = newGlobalDof;
	    mapGlobalLocalDOFs[newGlobalDof] = newLocalDof;
692
	    isRankDOF[newLocalDof] = false;
693
	    found = true;
694
695
696
	    break;
	  }
	}
697

Thomas Witkowski's avatar
Thomas Witkowski committed
698
	TEST_EXIT_DBG(found)("Should not happen!\n");
699
700
701
702
703
704
705
      }

      delete [] recvBuffers[i];
    }
  }


706
707
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, 
						      int& nOverallDOFs)
708
  {
709
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
710

711
    std::set<const DegreeOfFreedom*> rankDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
712
    DofToRank newBoundaryDOFs;
713
    
Thomas Witkowski's avatar
Thomas Witkowski committed
714

715
    // === Get all DOFs in ranks partition. ===
716

Thomas Witkowski's avatar
Thomas Witkowski committed
717
718
    ElementDofIterator elDofIt(feSpace->getAdmin());

719
720
721
722
723
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *element = elInfo->getElement();
      
Thomas Witkowski's avatar
Thomas Witkowski committed
724
725
726
727
      elDofIt.reset(element);
      do {
	rankDOFs.insert(elDofIt.getDofPtr());
      } while(elDofIt.next());
728
729
730
731

      elInfo = stack.traverseNext(elInfo);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
    RankToDofContainer sendNewDofs;
    RankToDofContainer recvNewDofs;
737

738
    for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
	   myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end();
741
	 ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end();
	   ++boundIt) {
	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
765
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
766
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
767
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
768
769
770
771
	  ("Should never happen!\n");

	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];
Thomas Witkowski's avatar
Thomas Witkowski committed
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
	
	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);
	
	DofContainer boundDOFs;
	addAllVertexDOFs(boundIt->rankObject.el, 
			 boundIt->rankObject.ithObjAtBoundary,
			 boundDOFs);
	addAllEdgeDOFs(boundIt->rankObject.el, 
		       boundIt->rankObject.ithObjAtBoundary,
		       boundDOFs);
787
788
789

	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
790
	  sendNewDofs[it->first].push_back(boundDOFs[i]);
791
	}
792
793
794
      }
    }    

795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
    for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
	   otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end();
	 ++it) {
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end();
	   ++boundIt) {
	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
822
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
823
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
824
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
825
826
827
828
829
830
	  ("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
831

Thomas Witkowski's avatar
Thomas Witkowski committed
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
	if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof2) ==
	    recvNewDofs[it->first].end())
	  recvNewDofs[it->first].push_back(dof2);
	if (find(recvNewDofs[it->first].begin(), recvNewDofs[it->first].end(), dof1) ==
	    recvNewDofs[it->first].end())
	  recvNewDofs[it->first].push_back(dof1);

	DofContainer boundDOFs;
	addAllEdgeDOFs(boundIt->rankObject.el, 
		       boundIt->rankObject.ithObjAtBoundary,
		       boundDOFs);
	addAllVertexDOFs(boundIt->rankObject.el, 
			 boundIt->rankObject.ithObjAtBoundary,
			 boundDOFs);

Thomas Witkowski's avatar
Thomas Witkowski committed
847
848
	for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) {
	//	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
849
850
	  rankDOFs.erase(boundDOFs[i]);
	  newBoundaryDOFs[boundDOFs[i]] = it->first;
851
	  recvNewDofs[it->first].push_back(boundDOFs[i]);
852
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
853
854
855
      }
    }

856
857
858
859
860
    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
861
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
862
863
864
865
866
    rstart -= nRankDOFs;


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

Thomas Witkowski's avatar
Thomas Witkowski committed
867
    mpiComm.Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888


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

    mapLocalGlobalDOFs.clear();
    mapGlobalLocalDOFs.clear();
    isRankDOF.clear();

    int i = 0;
    for (std::set<const DegreeOfFreedom*>::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end();
	 ++dofIt, i++) {
      const_cast<DegreeOfFreedom*>(*dofIt)[0] = i;
      mapLocalGlobalDOFs[i] = rstart + i;
      mapGlobalLocalDOFs[rstart + i] = i;
      isRankDOF[i] = true;
    }


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

889
890
891
    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

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

895
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
896
    for (RankToDofContainer::iterator sendIt = sendNewDofs.begin();
897
898
899
900
901
	 sendIt != sendNewDofs.end(); 
	 ++sendIt, i++) {
      int nSendDofs = sendIt->second.size();
      sendBuffers[i] = new int[nSendDofs];
      int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
902
      for (DofContainer::iterator dofIt = sendIt->second.begin();
903
904
905
	   dofIt != sendIt->second.end();
	   ++dofIt)
	sendBuffers[i][c++] = (*dofIt)[0] + rstart;
906

Thomas Witkowski's avatar
Thomas Witkowski committed
907
908
      request[requestCounter++] = 
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
909
910
911
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
912
913
914
    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end(); 
	 ++recvIt, i++) {
915
916
917
      int nRecvDofs = recvIt->second.size();
      recvBuffers[i] = new int[nRecvDofs];
	
Thomas Witkowski's avatar
Thomas Witkowski committed
918
919
      request[requestCounter++] = 
	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
920
921
    }

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

924
925

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
926
    for (RankToDofContainer::iterator sendIt = sendNewDofs.begin();
927
928
929
930
931
932
	 sendIt != sendNewDofs.end(); 
	 ++sendIt)
      delete [] sendBuffers[i++];

    i = 0;
    int newDofIndex = nRankDOFs;
Thomas Witkowski's avatar
Thomas Witkowski committed
933
934
935
    for (RankToDofContainer::iterator recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end(); 
	 ++recvIt) {      
936
      int j = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
937
      for (DofContainer::iterator dofIt = recvIt->second.begin();
938
939
	   dofIt != recvIt->second.end();
	   ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
940
	(*const_cast<DegreeOfFreedom*>(*dofIt)) = newDofIndex;
941
942
943
944
945
946
947
948
	mapLocalGlobalDOFs[newDofIndex] = recvBuffers[i][j];
	mapGlobalLocalDOFs[recvBuffers[i][j]] = newDofIndex;
	isRankDOF[newDofIndex] = false;
	newDofIndex++;
	j++;
      }

      delete [] recvBuffers[i++];
949
950
    }

951
952
953
954
    // === Update list of dofs that must be communicated for solution exchange. ===

    sendDofs = sendNewDofs;
    recvDofs = recvNewDofs;
955
956
957
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
958
959
  void ParallelDomainBase::addAllVertexDOFs(Element *el, int ithEdge, 
					    DofContainer& dofs)
960
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
961
    FUNCNAME("ParallelDomainBase::addAllVertexDOFs()");
962

Thomas Witkowski's avatar
Thomas Witkowski committed
963
964
965
    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
966
	addAllVertexDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
967
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
968
	addAllVertexDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
969
970
971
972
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
973
	addAllVertexDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
974
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
975
	addAllVertexDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
976
977
978
979
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
980
	addAllVertexDOFs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
981
	dofs.push_back(el->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
982
	addAllVertexDOFs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
983
984
985
986
987
988
989
990
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
  void ParallelDomainBase::addAllEdgeDOFs(Element *el, int ithEdge,
					  DofContainer& dofs)
  {
    FUNCNAME("ParallelDomainBase::addAllEdgeDOFs()");

    bool addThisEdge = false;

    switch (ithEdge) {
    case 0:
      if (el->getSecondChild())
	addAllEdgeDOFs(el->getSecondChild(), 2, dofs);
      else
	addThisEdge = true;

      break;
    case 1:
      if (el->getFirstChild())
	addAllEdgeDOFs(el->getFirstChild(), 2, dofs);
      else
	addThisEdge = true;

      break;
    case 2:
      if (el->getFirstChild()) {
	addAllEdgeDOFs(el->getFirstChild(), 1, dofs);
	addAllEdgeDOFs(el->getFirstChild(), 0, dofs);
      } else {
	addThisEdge = true;
      } 

      break;
    default:
      ERROR_EXIT("Should never happen!\n");
    }

    if (addThisEdge) {
      ElementDofIterator elDofIter(feSpace->getAdmin());
      elDofIter.reset(el);
      do {
	if (elDofIter.getCurrentPos() == 1 &&
	    elDofIter.getCurrentElementPos() == ithEdge) {
	  dofs.push_back(elDofIter.getDofPtr());
	}
      } while(elDofIter.next());      
    }
  }


  void ParallelDomainBase::createDOFMemberInfo(DofToPartitions& partitionDOFs,
					       DofContainer& rankDOFs,
					       DofToRank& boundaryDOFs)
1042
  {
1043
    // === Determine to each dof the set of partitions the dof belongs to. ===
1044

1045
1046
    ElementDofIterator elDofIt(feSpace->getAdmin());

1047
1048
1049
1050
1051
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *element = elInfo->getElement();

1052
1053
1054
1055
1056
1057
      elDofIt.reset(element);
      do {
	// Determine to each dof the partition(s) it corresponds to.
	partitionDOFs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
      } while(elDofIt.next());

1058
1059
1060
      elInfo = stack.traverseNext(elInfo);
    }

1061

1062
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
1063

1064
    // iterate over all DOFs
Thomas Witkowski's avatar
Thomas Witkowski committed
1065
    for (DofToPartitions::iterator it = partitionDOFs.begin();
1066
1067
	 it != partitionDOFs.end();
	 ++it) {
1068
1069

      // iterate over all partition the current DOF is part of.
1070
1071
1072
      for (std::set<int>::iterator itpart1 = it->second.begin();
	   itpart1 != it->second.end();
	   ++itpart1) {
1073

1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
	if (*itpart1 == mpiRank) {
	  if (it->second.size() == 1) {
	    rankDOFs.push_back(it->first);
	  } else {	    
	    // This dof is at the ranks boundary. It is owned by the rank only if
	    // the rank number is the highest of all ranks containing this dof.

	    bool insert = true;
	    int highestRank = mpiRank;
	    for (std::set<int>::iterator itpart2 = it->second.begin();
		 itpart2 != it->second.end();
		 ++itpart2) {
	      if (*itpart2 > mpiRank)
		insert = false;

	      if (*itpart2 > highestRank)
		highestRank = *itpart2;
	    }

	    if (insert)
	      rankDOFs.push_back(it->first);

	    boundaryDOFs[it->first] = highestRank;
	  }

	  break;
	}
1101

1102
1103
1104
1105
1106
      }
    }
  }


1107
  void ParallelDomainBase::testInteriorBoundary()