ParallelDomainProblem.cc 30.5 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
94
95
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);
    std::cout << "GR = " << globalRefinement << std::endl;

    refinementManager->globalRefine(mesh, globalRefinement);
96

97
    updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
98

99
    // === Create petsc matrix. ===
100

101
102
    int ierr;
    ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
103
    ierr = MatSetSizes(petscMatrix, nRankDOFs, nRankDOFs, nOverallDOFs, nOverallDOFs);
104
105
106
    ierr = MatSetType(petscMatrix, MATAIJ);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
107
    ierr = VecSetSizes(petscRhsVec, nRankDOFs, nOverallDOFs);
108
    ierr = VecSetType(petscRhsVec, VECMPI);
109
110

    ierr = VecCreate(PETSC_COMM_WORLD, &petscSolVec);
111
    ierr = VecSetSizes(petscSolVec, nRankDOFs, nOverallDOFs);
112
    ierr = VecSetType(petscSolVec, VECMPI);
113
114
115
  }

  void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
116
  {}
117

118

119
120
121
  void ParallelDomainProblemBase::fillPetscMatrix(DOFMatrix *mat, 
						  DOFVector<double> *vec)
  {
122
123
124
125
126
127
128
129
130
131
132
133
134
135
    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) {
136
137
	  int r = mapLocalGlobalDOFs[row(*icursor)];
	  int c = mapLocalGlobalDOFs[col(*icursor)];
138
	  double v = value(*icursor);
139

140
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
141
	}
142

143
144
145
146

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

147
    
148
149
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
150
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
151
152
153
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
154
    }    
155
156
  }

157

158
159
  void ParallelDomainProblemBase::solvePetscMatrix(DOFVector<double> *vec)
  {
160
    clock_t t = clock();
Thomas Witkowski's avatar
Thomas Witkowski committed
161
    
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
    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)
180
      *dofIt = vecPointer[counter++];
181
182
183

    VecRestoreArray(petscSolVec, &vecPointer);

184
185
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
186
187
188

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
189
190
    
    int i = 0;
191
192
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   sendIt = sendDofs.begin();
193
194
	 sendIt != sendDofs.end();
	 ++sendIt, i++) {
195
196
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
197

198
199
      for (int j = 0; j < nSendDOFs; j++)
	sendBuffers[i][j] = (*vec)[(sendIt->second)[j][0]];
200

Thomas Witkowski's avatar
Thomas Witkowski committed
201
202
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
203
204
205
    }

    i = 0;
206
207
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvDofs.begin();
208
209
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
210
211
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
212

Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
215
216
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
217
218
219

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

220
221
    
    i = 0;
222
223
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvDofs.begin();
224
225
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
226
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
227
	(*vec)[(recvIt->second)[j][0]] = recvBuffers[i][j];
228
229
230
231

      delete [] recvBuffers[i];
    }
    
232
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
233
234
235
      delete [] sendBuffers[i];

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

238

239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
  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();
	}
	TEST_EXIT(elNum != -1)("invalid element number\n");
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

273

274
275
276
277
278
279
280
281
282
283
284
285
286
287
  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);
  }

288

289
290
  void ParallelDomainProblemBase::createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
							     std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
291
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
292
293
294
295
    FUNCNAME("ParallelDomainProblemBase::createInteriorBoundaryInfo()");

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

296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
    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) {
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
	    // 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;

340
	    // === And add the part of the interior boundary. ===
341
342

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
345
346
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

347
348
349
350
351
352
	    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;
353
354
355
356
357
358
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
361
362
363
364
365
366
367


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
400

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

Thomas Witkowski's avatar
Thomas Witkowski committed
401
402

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

	  // 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];
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
  }


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


456
457
458
459
460
  void ParallelDomainProblemBase::createLocalGlobalNumbering(
                                 std::vector<const DegreeOfFreedom*>& rankDOFs,
				 std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
				 int& nRankDOFs, 
				 int& nOverallDOFs)
461
  {
462
463
464
    FUNCNAME("ParallelDomainProblemBase::createLocalGlobalNumbering()");

    // === Get rank information about DOFs. ===
465
466
467
468
469
470
471
472
473

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

474
475

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
478
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
479
480
481
    rstart -= nRankDOFs;

   
482
483
484
485
486
487
    // === 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;
488

489
490
491
492
    // 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;
493
494
495
496
497
498
499
500
501

    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
502
503
	DegreeOfFreedom newDofIndex = 0;
	for (int i = 0; i < nRankDOFs; i++) {
504
505
506
507
508
509
510
511
512
513
514
	  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)
515
	    sendNewDofs[*itRanks][it->first] = newDofIndex;
516
517
518
519
	}
      } 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.
520
521
522
523
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
524
525
526
527
      }
    }


528
    // === Send and receive the dof indices at boundary. ===
529
530
531

    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());
532
533
534

    sendDofs.clear();
    recvDofs.clear();
535
    
Thomas Witkowski's avatar
Thomas Witkowski committed
536
537
538
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

539
    int i = 0;
540
541
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
542
543
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
544
545
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
546
      int c = 0;
547
548
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
549
	   dofIt != sendIt->second.end();
550
551
552
	   ++dofIt) {
	sendBuffers[i][c++] = (dofIt->first)[0];
	sendBuffers[i][c++] = dofIt->second;
553

554
	sendDofs[sendIt->first].push_back(dofIt->first);
555
556
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
557
558
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
559
560
561
    }

    i = 0;
562
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
563
564
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
565
566
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
567

Thomas Witkowski's avatar
Thomas Witkowski committed
568
569
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
570
571
572
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
573
    MPI::Request::Waitall(requestCounter, request);
574
575

    
576
    // === Delete send buffers. ===
577
578

    i = 0;
579
580
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
581
582
583
584
585
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) 
      delete [] sendBuffers[i];


586
587
588
589
590
    // === Change dof indices for rank partition. ===

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

592
    for (int i = 0; i < nRankDOFs; i++) {
593
594
595
596
597
598
599
      const_cast<DegreeOfFreedom*>(rankDOFs[i])[0] = i; 
      mapLocalGlobalDOFs[i] = rstart + i;
      mapGlobalLocalDOFs[rstart + i] = i;
      isRankDOF[i] = true;
    }


600
    // === Change dof indices at boundary from other ranks. ===
601
602

    i = 0;
603
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
604
605
606
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

607
      for (int j = 0; j < recvIt->second; j++) {
608

609
610
611
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
	DegreeOfFreedom newLocalDof = mapLocalGlobalDOFs.size();
612

613
	bool found = false;
614
615
616
617

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

619
	  if ((dofIt->first)[0] == oldDof) {
620
	    recvDofs[recvIt->first].push_back(dofIt->first);
621
	    const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newLocalDof;
622
623
	    mapLocalGlobalDOFs[newLocalDof] = newGlobalDof;
	    mapGlobalLocalDOFs[newGlobalDof] = newLocalDof;
624
	    isRankDOF[newLocalDof] = false;
625
	    found = true;
626
627
628
	    break;
	  }
	}
629
630

	TEST_EXIT(found)("Should not happen!\n");
631
632
633
634
635
636
637
      }

      delete [] recvBuffers[i];
    }
  }


638
639
640
  void ParallelDomainProblemBase::updateLocalGlobalNumbering(int& nRankDOFs, 
							     int& nOverallDOFs)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
641
642
    FUNCNAME("ParallelDomainProblemBase::updateLocalGlobalNumbering()");

643
    std::set<const DegreeOfFreedom*> rankDOFs;
644
    std::map<const DegreeOfFreedom*, int> newBoundaryDOFs;
645
    
646
    // === Get all DOFs in ranks partition. ===
647
648
649
650
651
652
653
654
655
656
657
658
659

    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
660
661
    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
    // === rankDOFs to boundaryDOFs                                                ===
662

663
664
    std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDofs = sendDofs;
    std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDofs = recvDofs;
665

666
    for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
Thomas Witkowski's avatar
Thomas Witkowski committed
667
668
	   myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end();
669
	 ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
      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");
	}

693
694
695
696
697
698
699
700
	TEST_EXIT(boundaryDOFs.find(dof1) != boundaryDOFs.end())
	  ("Should never happen!\n");
	TEST_EXIT(boundaryDOFs.find(dof2) != boundaryDOFs.end())
	  ("Should never happen!\n");

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

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

711

712
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
745
746
747
    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");
	}

	TEST_EXIT(boundaryDOFs.find(dof1) != boundaryDOFs.end())
	  ("Should never happen!\n");
	TEST_EXIT(boundaryDOFs.find(dof2) != boundaryDOFs.end())
	  ("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
748

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

    }

761
762
763
764
765
    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
766
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
767
768
769
770
771
    rstart -= nRankDOFs;


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

Thomas Witkowski's avatar
Thomas Witkowski committed
772
    mpiComm.Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793


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

794
795
796
    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

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

800
801
802
803
804
805
806
807
808
809
810
811
    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;
812

Thomas Witkowski's avatar
Thomas Witkowski committed
813
814
      request[requestCounter++] = 
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
815
816
817
818
819
820
821
822
823
    }

    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
824
825
      request[requestCounter++] = 
	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
826
827
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
828
829
830

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

831
832
833
834
835
836
837
838
839
840
841
842
843

    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) {
844
      
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
      int j = 0;
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = recvIt->second.begin();
	   dofIt != recvIt->second.end();
	   ++dofIt) {
	const_cast<DegreeOfFreedom*>(*dofIt)[0] = newDofIndex;
// 	if (mpiRank == 1) 
// 	  std::cout << "SET TO " << newDofIndex << std::endl;

	mapLocalGlobalDOFs[newDofIndex] = recvBuffers[i][j];
	mapGlobalLocalDOFs[recvBuffers[i][j]] = newDofIndex;
	isRankDOF[newDofIndex] = false;
	newDofIndex++;
	j++;
      }

      delete [] recvBuffers[i++];
861
862
    }

863
864
865
866
    // === Update list of dofs that must be communicated for solution exchange. ===

    sendDofs = sendNewDofs;
    recvDofs = recvNewDofs;
867
868
869
870
871
872
873
874
  }


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

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


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

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

922
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
923

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

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

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

962
963
964
965
966
      }
    }
  }


967
968
969
  ParallelDomainProblemScal::ParallelDomainProblemScal(const std::string& name,
						       ProblemScal *problem,
						       ProblemInstatScal *problemInstat)
970
971
972
973
974
    : ParallelDomainProblemBase(name, 
				problem, 
				problemInstat, 
				problem->getFESpace(),
				problem->getRefinementManager()),
975
      probScal(problem)
976
977
978
  {
  }

979
980
981
982
  void ParallelDomainProblemScal::initParallelization(AdaptInfo *adaptInfo)
  {
    ParallelDomainProblemBase::initParallelization(adaptInfo);

Thomas Witkowski's avatar
Thomas Witkowski committed
983
984
985
986
987
    DOFMatrix* m = probScal->getSystemMatrix();

    TEST_EXIT(m)("No DOF Matrix!\n");

    m->setIsRankDOF(isRankDOF);
988
989
  }

990
991
992
993
994
995
996
997
  Flag ParallelDomainProblemScal::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    //    return iterationIF->oneIteration(adaptInfo, toDo);

    Flag flag = dynamic_cast<StandardProblemIteration*>(iterationIF)->buildAndAdapt(adaptInfo, toDo);

    fillPetscMatrix(probScal->getSystemMatrix(), probScal->getRHS());

998
999
    solvePetscMatrix(probScal->getSolution());

1000
//     if (toDo.isSet(SOLVE))
For faster browsing, not all history is shown. View entire blame