ParallelDomainProblem.cc 31 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
14
15
#include "VtkWriter.h"

#include "petscksp.h"
16
17
18
19

namespace AMDiS {

  ParallelDomainProblemBase::ParallelDomainProblemBase(const std::string& name,
20
21
						       ProblemIterationInterface *iIF,
						       ProblemTimeInterface *tIF,
22
23
						       FiniteElemSpace *fe,
						       RefinementManager *refineManager)
24
25
    : iterationIF(iIF),
      timeIF(tIF),
26
27
      feSpace(fe),
      mesh(fe->getMesh()),
28
      refinementManager(refineManager),
29
      initialPartitionMesh(true),
30
      nRankDOFs(0)
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
  {
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
  }

  void ParallelDomainProblemBase::initParallelization(AdaptInfo *adaptInfo)
  {
    if (mpiSize <= 1)
      return;

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


51
    // === Create new global and local DOF numbering. ===
52

53
54
55
    // Set of all DOFs of the rank.
    std::vector<const DegreeOfFreedom*> rankDOFs;
    // Number of DOFs in ranks partition that are owned by the rank.
56
    int nRankDOFs = 0;
57
    // Number of all DOFs in the macro mesh.
58
    int nOverallDOFs = 0;
59
60

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

Thomas Witkowski's avatar
Thomas Witkowski committed
62
63
64

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

65
    createInteriorBoundaryInfo(rankDOFs, boundaryDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
66

Thomas Witkowski's avatar
Thomas Witkowski committed
67

68
69
    // === Remove all macro elements that are not part of the rank partition. ===

70
    removeMacroElements();
71
      
72

73
    // === Reset all DOFAdmins of the mesh. ===
74
75
76

    int nAdmins = mesh->getNumberOfDOFAdmin();
    for (int i = 0; i < nAdmins; i++) {
77
78
79
80
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));

      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
81
      for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
82
83
 	admin.setDOFFree(j, false);

84
      admin.setUsedSize(mapLocalGlobalDOFs.size());
85
86
      admin.setUsedCount(mapLocalGlobalDOFs.size());
      admin.setFirstHole(mapLocalGlobalDOFs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
87
88
    }

89
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
90

Thomas Witkowski's avatar
Thomas Witkowski committed
91
92
93
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

Thomas Witkowski's avatar
Thomas Witkowski committed
94
95
    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);
96

Thomas Witkowski's avatar
Thomas Witkowski committed
97
98
99
100
101
102
103
      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
    }


#if (DEBUG != 0)
    testInteriorBoundary();
#endif
104

105
    // === Create petsc matrix. ===
106

107
108
    int ierr;
    ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
109
    ierr = MatSetSizes(petscMatrix, nRankDOFs, nRankDOFs, nOverallDOFs, nOverallDOFs);
110
111
112
    ierr = MatSetType(petscMatrix, MATAIJ);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
113
    ierr = VecSetSizes(petscRhsVec, nRankDOFs, nOverallDOFs);
114
    ierr = VecSetType(petscRhsVec, VECMPI);
115
116

    ierr = VecCreate(PETSC_COMM_WORLD, &petscSolVec);
117
    ierr = VecSetSizes(petscSolVec, nRankDOFs, nOverallDOFs);
118
    ierr = VecSetType(petscSolVec, VECMPI);
119
120
121
  }

  void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
122
  {}
123

124

125
126
127
  void ParallelDomainProblemBase::fillPetscMatrix(DOFMatrix *mat, 
						  DOFVector<double> *vec)
  {
128
129
130
131
132
133
134
135
136
137
138
139
140
141
    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;

    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)
	if (value(*icursor) != 0.0) {
142
143
	  int r = mapLocalGlobalDOFs[row(*icursor)];
	  int c = mapLocalGlobalDOFs[col(*icursor)];
144
	  double v = value(*icursor);
145

146
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
147
	}
148

149
150
151
152

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

153
    
154
155
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
156
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
157
158
159
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
160
    }    
161
162
  }

163

164
165
  void ParallelDomainProblemBase::solvePetscMatrix(DOFVector<double> *vec)
  {
166
    clock_t t = clock();
Thomas Witkowski's avatar
Thomas Witkowski committed
167
    
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
    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);
    KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, 0);
    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)
186
      *dofIt = vecPointer[counter++];
187
188
189

    VecRestoreArray(petscSolVec, &vecPointer);

190
191
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
192
193
194

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
195
196
    
    int i = 0;
197
198
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   sendIt = sendDofs.begin();
199
200
	 sendIt != sendDofs.end();
	 ++sendIt, i++) {
201
202
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
203

204
205
      for (int j = 0; j < nSendDOFs; j++)
	sendBuffers[i][j] = (*vec)[(sendIt->second)[j][0]];
206

Thomas Witkowski's avatar
Thomas Witkowski committed
207
208
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
209
210
211
    }

    i = 0;
212
213
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvDofs.begin();
214
215
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
216
217
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
218

Thomas Witkowski's avatar
Thomas Witkowski committed
219
220
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
221
222
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
223
224
225

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

226
227
    
    i = 0;
228
229
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvDofs.begin();
230
231
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
232
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
233
	(*vec)[(recvIt->second)[j][0]] = recvBuffers[i][j];
234
235
236
237

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

    std::cout << "SOLUTION = " << TIME_USED(t,clock()) << std::endl;
242
243
  }

244

245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
  double ParallelDomainProblemBase::setElemWeights(AdaptInfo *adaptInfo) 
  {
    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
266
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
267
268
269
270
271
272
273
274
275
276
277
278
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

279

280
281
282
283
284
285
286
287
288
289
290
291
292
293
  void ParallelDomainProblemBase::partitionMesh(AdaptInfo *adaptInfo)
  {
    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);
  }

294

295
296
  void ParallelDomainProblemBase::createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
							     std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
297
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
300
301
    FUNCNAME("ParallelDomainProblemBase::createInteriorBoundaryInfo()");

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

302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
    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));
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
	    // We have found an element that is at an interior boundary. 

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

	    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;

346
	    // === And add the part of the interior boundary. ===
347
348

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
351
352
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

353
354
355
356
357
358
	    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;
359
360
361
362
363
364
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
367
368
369
370
371
372
373


    // === 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
374
375
376
377
378
379
    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
380
381
	 rankIt != myIntBoundary.boundary.end();
	 ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
384
      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
      for (int i = 0; i < nSendInt; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
387
	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
388
389
      request[requestCounter++] =
	mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
390
391
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = 
	   otherIntBoundary.boundary.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
394
395
	 rankIt != otherIntBoundary.boundary.end();
	 ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
396
397
      int nRecvInt = rankIt->second.size();
      int *buffer = new int[nRecvInt];
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
      recvBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
      request[requestCounter++] = 
	mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
402
403
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
404
405
406

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

Thomas Witkowski's avatar
Thomas Witkowski committed
407
408

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
409
410
    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = 
	   otherIntBoundary.boundary.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
	 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
427
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
428
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
429
430
431
432
433
434
435
436
437
438
439
440
441

	  // 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];
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
  }


  void ParallelDomainProblemBase::removeMacroElements()
  {
    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);
  }


462
463
464
465
466
  void ParallelDomainProblemBase::createLocalGlobalNumbering(
                                 std::vector<const DegreeOfFreedom*>& rankDOFs,
				 std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
				 int& nRankDOFs, 
				 int& nOverallDOFs)
467
  {
468
469
470
    FUNCNAME("ParallelDomainProblemBase::createLocalGlobalNumbering()");

    // === Get rank information about DOFs. ===
471
472
473
474
475
476
477
478
479

    // 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();

480
481

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
484
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
485
486
487
    rstart -= nRankDOFs;

   
488
489
490
491
492
493
    // === 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;
494

495
496
497
498
    // 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;
499
500
501
502
503
504
505
506
507

    for (std::map<const DegreeOfFreedom*, int>::iterator it = boundaryDOFs.begin();
	 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
508
509
	DegreeOfFreedom newDofIndex = 0;
	for (int i = 0; i < nRankDOFs; i++) {
510
511
512
513
514
515
516
517
518
519
520
	  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)
521
	    sendNewDofs[*itRanks][it->first] = newDofIndex;
522
523
524
525
	}
      } 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.
526
527
528
529
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
530
531
532
533
      }
    }


534
    // === Send and receive the dof indices at boundary. ===
535
536
537

    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());
538
539
540

    sendDofs.clear();
    recvDofs.clear();
541
    
Thomas Witkowski's avatar
Thomas Witkowski committed
542
543
544
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

545
    int i = 0;
546
547
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
548
549
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
550
551
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
552
      int c = 0;
553
554
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
555
	   dofIt != sendIt->second.end();
556
557
558
	   ++dofIt) {
	sendBuffers[i][c++] = (dofIt->first)[0];
	sendBuffers[i][c++] = dofIt->second;
559

560
	sendDofs[sendIt->first].push_back(dofIt->first);
561
562
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
563
564
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
565
566
567
    }

    i = 0;
568
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
569
570
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
571
572
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
573

Thomas Witkowski's avatar
Thomas Witkowski committed
574
575
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
576
577
578
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
579
    MPI::Request::Waitall(requestCounter, request);
580
581

    
582
    // === Delete send buffers. ===
583
584

    i = 0;
585
586
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
587
588
589
590
591
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) 
      delete [] sendBuffers[i];


592
593
594
595
596
    // === Change dof indices for rank partition. ===

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

598
    for (int i = 0; i < nRankDOFs; i++) {
599
600
601
602
603
604
605
      const_cast<DegreeOfFreedom*>(rankDOFs[i])[0] = i; 
      mapLocalGlobalDOFs[i] = rstart + i;
      mapGlobalLocalDOFs[rstart + i] = i;
      isRankDOF[i] = true;
    }


606
    // === Change dof indices at boundary from other ranks. ===
607
608

    i = 0;
609
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
610
611
612
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

613
      for (int j = 0; j < recvIt->second; j++) {
614

615
616
617
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
	DegreeOfFreedom newLocalDof = mapLocalGlobalDOFs.size();
618

619
	bool found = false;
620
621
622
623

	for (std::map<const DegreeOfFreedom*, int>::iterator dofIt = boundaryDOFs.begin();
	     dofIt != boundaryDOFs.end();
	     ++dofIt) {
624

625
	  if ((dofIt->first)[0] == oldDof) {
626
	    recvDofs[recvIt->first].push_back(dofIt->first);
627
	    const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newLocalDof;
628
629
	    mapLocalGlobalDOFs[newLocalDof] = newGlobalDof;
	    mapGlobalLocalDOFs[newGlobalDof] = newLocalDof;
630
	    isRankDOF[newLocalDof] = false;
631
	    found = true;
632
633
634
	    break;
	  }
	}
635

Thomas Witkowski's avatar
Thomas Witkowski committed
636
	TEST_EXIT_DBG(found)("Should not happen!\n");
637
638
639
640
641
642
643
      }

      delete [] recvBuffers[i];
    }
  }


644
645
646
  void ParallelDomainProblemBase::updateLocalGlobalNumbering(int& nRankDOFs, 
							     int& nOverallDOFs)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
647
648
    FUNCNAME("ParallelDomainProblemBase::updateLocalGlobalNumbering()");

649
    std::set<const DegreeOfFreedom*> rankDOFs;
650
    std::map<const DegreeOfFreedom*, int> newBoundaryDOFs;
651
    
652
    // === Get all DOFs in ranks partition. ===
653
654
655
656
657
658
659
660
661
662
663
664
665

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *element = elInfo->getElement();
      
      for (int i = 0; i < 3; i++) 
	rankDOFs.insert(element->getDOF(i));

      elInfo = stack.traverseNext(elInfo);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
666
667
    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
    // === rankDOFs to boundaryDOFs                                                ===
668

669
670
    std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDofs = sendDofs;
    std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDofs = recvDofs;
671

672
    for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
Thomas Witkowski's avatar
Thomas Witkowski committed
673
674
	   myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end();
675
	 ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
      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
699
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
700
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
701
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
702
703
704
705
706
	  ("Should never happen!\n");

	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
	newBoundaryDOFs[dof2] = boundaryDOFs[dof2];

Thomas Witkowski's avatar
Thomas Witkowski committed
707
708
709
	std::vector<const DegreeOfFreedom*> boundDOFs;
	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		   boundDOFs);
710
711
	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
712
	  sendNewDofs[it->first].push_back(boundDOFs[i]);
713
	}
714
715
716
      }
    }    

717

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
    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
745
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
746
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
747
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
748
749
750
751
752
753
	  ("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
754

755
756
757
	std::vector<const DegreeOfFreedom*> boundDOFs;
	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		   boundDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
758
759
	for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) {
	//	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
760
761
	  rankDOFs.erase(boundDOFs[i]);
	  newBoundaryDOFs[boundDOFs[i]] = it->first;
762
	  recvNewDofs[it->first].push_back(boundDOFs[i]);
763
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
764
765
766
767
      }

    }

768
769
770
771
772
    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
773
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
774
775
776
777
778
    rstart -= nRankDOFs;


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

Thomas Witkowski's avatar
Thomas Witkowski committed
779
    mpiComm.Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800


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

801
802
803
    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

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

807
808
809
810
811
812
813
814
815
816
817
818
    i = 0;
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   sendIt = sendNewDofs.begin();
	 sendIt != sendNewDofs.end(); 
	 ++sendIt, i++) {
      int nSendDofs = sendIt->second.size();
      sendBuffers[i] = new int[nSendDofs];
      int c = 0;
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = sendIt->second.begin();
	   dofIt != sendIt->second.end();
	   ++dofIt)
	sendBuffers[i][c++] = (*dofIt)[0] + rstart;
819

Thomas Witkowski's avatar
Thomas Witkowski committed
820
821
      request[requestCounter++] = 
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
822
823
824
825
826
827
828
829
830
    }

    i = 0;
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
      int nRecvDofs = recvIt->second.size();
      recvBuffers[i] = new int[nRecvDofs];
	
Thomas Witkowski's avatar
Thomas Witkowski committed
831
832
      request[requestCounter++] = 
	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
833
834
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
835
836
837

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

838
839
840
841
842
843
844
845
846
847
848
849
850

    i = 0;
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   sendIt = sendNewDofs.begin();
	 sendIt != sendNewDofs.end(); 
	 ++sendIt)
      delete [] sendBuffers[i++];

    i = 0;
    int newDofIndex = nRankDOFs;
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end(); ++recvIt) {
851
      
852
853
854
855
856
857
858
859
860
861
862
863
864
      int j = 0;
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = recvIt->second.begin();
	   dofIt != recvIt->second.end();
	   ++dofIt) {
	const_cast<DegreeOfFreedom*>(*dofIt)[0] = newDofIndex;
	mapLocalGlobalDOFs[newDofIndex] = recvBuffers[i][j];
	mapGlobalLocalDOFs[recvBuffers[i][j]] = newDofIndex;
	isRankDOF[newDofIndex] = false;
	newDofIndex++;
	j++;
      }

      delete [] recvBuffers[i++];
865
866
    }

867
868
869
870
    // === Update list of dofs that must be communicated for solution exchange. ===

    sendDofs = sendNewDofs;
    recvDofs = recvNewDofs;
871
872
873
874
875
876
877
878
  }


  void ParallelDomainProblemBase::addAllDOFs(Element *el, int ithEdge, 
					     std::vector<const DegreeOfFreedom*>& dofs)
  {
    FUNCNAME("ParallelDomainProblemBase::addAllDOFs()");

Thomas Witkowski's avatar
Thomas Witkowski committed
879
880
881
    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
882
	addAllDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
883
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
884
	addAllDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
885
886
887
888
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
889
	addAllDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
890
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
891
	addAllDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
892
893
894
895
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
896
	addAllDOFs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
897
	dofs.push_back(el->getFirstChild()->getDOF(2));
898
	addAllDOFs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
899
900
901
902
903
904
905
906
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
  }


907
  void ParallelDomainProblemBase::createDOFMemberInfo(
908
		       std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
909
910
911
		       std::vector<const DegreeOfFreedom*>& rankDOFs,
		       std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
  {
912
    // === Determine to each dof the set of partitions the dof belongs to. ===
913
914
915
916
917
918
919
920
921
922
923
924
925

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

      // Determine to each dof the partition(s) it corresponds to.
      for (int i = 0; i < 3; i++) 
	partitionDOFs[element->getDOF(i)].insert(partitionVec[element->getIndex()]);
          
      elInfo = stack.traverseNext(elInfo);
    }

926
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
927

928
    // iterate over all DOFs
929
930
931
    for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDOFs.begin();
	 it != partitionDOFs.end();
	 ++it) {
932
933

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

938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
	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;
	}
965

966
967
968
969
970
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
  void ParallelDomainProblemBase::testInteriorBoundary()
  {
    // Maps to each neighbour rank an array of WorldVectors. This array contain the 
    // coordinates of all dofs this rank shares on the interior boundary with the 
    // neighbour rank.
    std::map<int, std::vector<WorldVector<double> > > sendCoords;
    std::map<int, int> nRecvCoords;
    
    for (std::map<const DegreeOfFreedom*, int>::iterator it = boundaryDOFs.begin();
	 it != boundaryDOFs.end();
	 ++it) {
      
    }
  }


987
988
989
  ParallelDomainProblemScal::ParallelDomainProblemScal(const std::string& name,
						       ProblemScal *problem,
						       ProblemInstatScal *problemInstat)
990
991
992
993
994
    : ParallelDomainProblemBase(name, 
				problem, 
				problemInstat, 
				problem->getFESpace(),
				problem->getRefinementManager()),
995
      probScal(problem)
996
997
998
  {
  }

999
1000
  void ParallelDomainProblemScal::initParallelization(AdaptInfo *adaptInfo)
  {