ParallelDomainProblem.cc 33.7 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
41
42
43
44
45
46
  {
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
  }

47
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
48
49
50
51
  {
    if (mpiSize <= 1)
      return;

52
53
54
55
56
    // 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();

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


65
    // === Create new global and local DOF numbering. ===
66

67
68
69
    // Set of all DOFs of the rank.
    std::vector<const DegreeOfFreedom*> rankDOFs;
    // Number of DOFs in ranks partition that are owned by the rank.
70
    int nRankDOFs = 0;
71
    // Number of all DOFs in the macro mesh.
72
    int nOverallDOFs = 0;
73
74

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

Thomas Witkowski's avatar
Thomas Witkowski committed
76
77
78

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

79
    createInteriorBoundaryInfo(rankDOFs, boundaryDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
80

Thomas Witkowski's avatar
Thomas Witkowski committed
81

82
83
    // === Remove all macro elements that are not part of the rank partition. ===

84
    removeMacroElements();
85
      
86

87
    // === Reset all DOFAdmins of the mesh. ===
88
89
90

    int nAdmins = mesh->getNumberOfDOFAdmin();
    for (int i = 0; i < nAdmins; i++) {
91
92
93
94
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));

      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
95
      for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
96
97
 	admin.setDOFFree(j, false);

98
      admin.setUsedSize(mapLocalGlobalDOFs.size());
99
100
      admin.setUsedCount(mapLocalGlobalDOFs.size());
      admin.setFirstHole(mapLocalGlobalDOFs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
101
102
    }

103
104
    exit(0);

105
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
106

Thomas Witkowski's avatar
Thomas Witkowski committed
107
108
109
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

Thomas Witkowski's avatar
Thomas Witkowski committed
110
111
    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);
112

Thomas Witkowski's avatar
Thomas Witkowski committed
113
114
115
116
      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
    }


117
118
119
// #if (DEBUG != 0)
//     testInteriorBoundary();
// #endif
120

121
    // === Create petsc matrix. ===
122

123
124
    int ierr;
    ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
125
    ierr = MatSetSizes(petscMatrix, nRankDOFs, nRankDOFs, nOverallDOFs, nOverallDOFs);
126
127
128
    ierr = MatSetType(petscMatrix, MATAIJ);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
129
    ierr = VecSetSizes(petscRhsVec, nRankDOFs, nOverallDOFs);
130
    ierr = VecSetType(petscRhsVec, VECMPI);
131
132

    ierr = VecCreate(PETSC_COMM_WORLD, &petscSolVec);
133
    ierr = VecSetSizes(petscSolVec, nRankDOFs, nOverallDOFs);
134
    ierr = VecSetType(petscSolVec, VECMPI);
135
136
  }

137
138
139
140

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

142

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

179
180
181
182
    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)
183
	if (value(*icursor) != 0.0) {
184
185
	  int r = mapLocalGlobalDOFs[row(*icursor)];
	  int c = mapLocalGlobalDOFs[col(*icursor)];
186
	  double v = value(*icursor);
187

188
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
189
	}
190

191
192
193
194

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

195
    
196
197
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
198
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
199
200
201
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
202
    }    
203
204
  }

205

206
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
207
208
209
210
211
212
213
214
215
216
  {
    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);
217
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
218
219
220
221
222
223
224
225
    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)
226
      *dofIt = vecPointer[counter++];
227
228
229

    VecRestoreArray(petscSolVec, &vecPointer);

230
231
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
232
233
234

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
235
236
    
    int i = 0;
237
238
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   sendIt = sendDofs.begin();
239
240
	 sendIt != sendDofs.end();
	 ++sendIt, i++) {
241
242
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
243

244
245
      for (int j = 0; j < nSendDOFs; j++)
	sendBuffers[i][j] = (*vec)[(sendIt->second)[j][0]];
246

Thomas Witkowski's avatar
Thomas Witkowski committed
247
248
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
249
250
251
    }

    i = 0;
252
253
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvDofs.begin();
254
255
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
256
257
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
258

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

Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
265

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

266
267
    
    i = 0;
268
269
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvDofs.begin();
270
271
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
272
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
273
	(*vec)[(recvIt->second)[j][0]] = recvBuffers[i][j];
274
275
276
277

      delete [] recvBuffers[i];
    }
    
278
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
279
      delete [] sendBuffers[i];
280
281
  }

282

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

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

317

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

332

333
334
  void ParallelDomainBase::createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
						      std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
335
  {
336
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
339

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

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

356
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
357
358
	    // We have found an element that is at an interior boundary. 

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

	    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;

386
	    // === And add the part of the interior boundary. ===
387
388

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

393
394
395
396
397
398
	    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;
399
400
401
402
403
404
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
407
408
409
410
411
412
413


    // === 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
414
415
416
417
418
419
    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
420
421
	 rankIt != myIntBoundary.boundary.end();
	 ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
424
      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
      for (int i = 0; i < nSendInt; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
427
	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
428
429
      request[requestCounter++] =
	mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
430
431
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
446

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

Thomas Witkowski's avatar
Thomas Witkowski committed
447
448

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = 
	   otherIntBoundary.boundary.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
	 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
467
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
468
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
469
470
471
472
473
474
475
476
477
478
479
480
481

	  // 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];
482
483
484
  }


485
  void ParallelDomainBase::removeMacroElements()
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
  {
    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);
  }


502
  void ParallelDomainBase::createLocalGlobalNumbering(
503
504
505
506
                                 std::vector<const DegreeOfFreedom*>& rankDOFs,
				 std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
				 int& nRankDOFs, 
				 int& nOverallDOFs)
507
  {
508
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
509
510

    // === Get rank information about DOFs. ===
511
512
513
514
515
516
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;

    createDOFMemberInfo(partitionDOFs, rankDOFs, boundaryDOFs);

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

520
521

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

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

   
528
529
530
531
532
533
    // === 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;
534

535
536
537
538
    // 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;
539
540
541
542
543
544
545
546
547

    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
548
549
	DegreeOfFreedom newDofIndex = 0;
	for (int i = 0; i < nRankDOFs; i++) {
550
551
552
553
554
555
556
557
558
559
560
	  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)
561
	    sendNewDofs[*itRanks][it->first] = newDofIndex;
562
563
564
565
	}
      } 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.
566
567
568
569
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
570
571
572
573
      }
    }


574
    // === Send and receive the dof indices at boundary. ===
575
576
577

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

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

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

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

601
	sendDofs[sendIt->first].push_back(dofIt->first);
602
603
      }

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

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

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


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

    
623
    // === Delete send buffers. ===
624
625

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


633
634
635
636
637
    // === Change dof indices for rank partition. ===

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

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


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

    i = 0;
650
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
651
652
653
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

654
      for (int j = 0; j < recvIt->second; j++) {
655

656
657
658
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
	DegreeOfFreedom newLocalDof = mapLocalGlobalDOFs.size();
659

660
	bool found = false;
661
662
663
664

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

666
	  if ((dofIt->first)[0] == oldDof) {
667
	    recvDofs[recvIt->first].push_back(dofIt->first);
668
	    const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newLocalDof;
669
670
	    mapLocalGlobalDOFs[newLocalDof] = newGlobalDof;
	    mapGlobalLocalDOFs[newGlobalDof] = newLocalDof;
671
	    isRankDOF[newLocalDof] = false;
672
	    found = true;
673
674
675
	    break;
	  }
	}
676

Thomas Witkowski's avatar
Thomas Witkowski committed
677
	TEST_EXIT_DBG(found)("Should not happen!\n");
678
679
680
681
682
683
684
      }

      delete [] recvBuffers[i];
    }
  }


685
686
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, 
						      int& nOverallDOFs)
687
  {
688
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
689

690
    std::set<const DegreeOfFreedom*> rankDOFs;
691
    std::map<const DegreeOfFreedom*, int> newBoundaryDOFs;
692
    
693
    // === Get all DOFs in ranks partition. ===
694
695
696
697
698
699
700
701
702
703
704
705
706

    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
707
    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
708
    // === rankDOFs to boundaryDOFs.                                               ===
709

710
711
    std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDofs = sendDofs;
    std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDofs = recvDofs;
712

713
    for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
	   myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end();
716
	 ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
      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
740
	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
741
	  ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
742
	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
743
744
745
746
747
	  ("Should never happen!\n");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
748
749
750
	std::vector<const DegreeOfFreedom*> boundDOFs;
	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		   boundDOFs);
751
752
	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
753
	  sendNewDofs[it->first].push_back(boundDOFs[i]);
754
	}
755
756
757
      }
    }    

758

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

796
797
798
	std::vector<const DegreeOfFreedom*> boundDOFs;
	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		   boundDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
799
800
	for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) {
	//	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
801
802
	  rankDOFs.erase(boundDOFs[i]);
	  newBoundaryDOFs[boundDOFs[i]] = it->first;
803
	  recvNewDofs[it->first].push_back(boundDOFs[i]);
804
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
805
806
807
808
      }

    }

809
810
811
812
813
    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
814
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
815
816
817
818
819
    rstart -= nRankDOFs;


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

Thomas Witkowski's avatar
Thomas Witkowski committed
820
    mpiComm.Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841


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

842
843
844
    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

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

848
849
850
851
852
853
854
855
856
857
858
859
    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;
860

Thomas Witkowski's avatar
Thomas Witkowski committed
861
862
      request[requestCounter++] = 
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
863
864
865
866
867
868
869
870
871
    }

    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
872
873
      request[requestCounter++] = 
	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
874
875
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
876
877
878

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

879
880
881
882
883
884
885
886
887
888
889
890
891

    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) {
892
      
893
894
895
896
897
898
899
900
901
902
903
904
905
      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++];
906
907
    }

908
909
910
911
    // === Update list of dofs that must be communicated for solution exchange. ===

    sendDofs = sendNewDofs;
    recvDofs = recvNewDofs;
912
913
914
  }


915
916
  void ParallelDomainBase::addAllDOFs(Element *el, int ithEdge, 
				      std::vector<const DegreeOfFreedom*>& dofs)
917
  {
918
    FUNCNAME("ParallelDomainBase::addAllDOFs()");
919

Thomas Witkowski's avatar
Thomas Witkowski committed
920
921
922
    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
923
	addAllDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
924
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
925
	addAllDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
926
927
928
929
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
930
	addAllDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
931
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
932
	addAllDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
933
934
935
936
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
937
	addAllDOFs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
938
	dofs.push_back(el->getFirstChild()->getDOF(2));
939
	addAllDOFs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
940
941
942
943
944
945
946
947
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
  }


948
  void ParallelDomainBase::createDOFMemberInfo(
949
		       std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
950
951
952
		       std::vector<const DegreeOfFreedom*>& rankDOFs,
		       std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
  {
953
    // === Determine to each dof the set of partitions the dof belongs to. ===
954

955
956
    ElementDofIterator elDofIt(feSpace->getAdmin());

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

962
963
964
965
966
967
      elDofIt.reset(element);
      do {
	// Determine to each dof the partition(s) it corresponds to.
	partitionDOFs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
      } while(elDofIt.next());

968
969
970
      elInfo = stack.traverseNext(elInfo);
    }

971

972
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
973

974
    // iterate over all DOFs
975
976
977
    for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDOFs.begin();
	 it != partitionDOFs.end();
	 ++it) {
978
979

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

984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
	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;
For faster browsing, not all history is shown. View entire blame