ParallelDomainProblem.cc 29.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

91
    refinementManager->globalRefine(mesh, 12);
92

93
    updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
94

95
    // === Create petsc matrix. ===
96

97
98
    int ierr;
    ierr = MatCreate(PETSC_COMM_WORLD, &petscMatrix);
99
    ierr = MatSetSizes(petscMatrix, nRankDOFs, nRankDOFs, nOverallDOFs, nOverallDOFs);
100
101
102
    ierr = MatSetType(petscMatrix, MATAIJ);

    ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
103
    ierr = VecSetSizes(petscRhsVec, nRankDOFs, nOverallDOFs);
104
    ierr = VecSetType(petscRhsVec, VECMPI);
105
106

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

  void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
112
  {}
113

114

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

136
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
137
	}
138

139
140
141
142

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

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

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
150
    }    
151
152
  }

153

154
155
  void ParallelDomainProblemBase::solvePetscMatrix(DOFVector<double> *vec)
  {
156
157
    clock_t t = clock();

158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
    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)
176
      *dofIt = vecPointer[counter++];
177
178
179

    VecRestoreArray(petscSolVec, &vecPointer);

180
181
182
183
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
    
    int i = 0;
184
185
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   sendIt = sendDofs.begin();
186
187
	 sendIt != sendDofs.end();
	 ++sendIt, i++) {
188
189
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
190

191
192
      for (int j = 0; j < nSendDOFs; j++)
	sendBuffers[i][j] = (*vec)[(sendIt->second)[j][0]];
193

194
      mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
195
196
197
    }

    i = 0;
198
199
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvDofs.begin();
200
201
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
202
203
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
204

205
      mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
206
207
208
209
210
211
    }

    
    mpiComm.Barrier();
    
    i = 0;
212
213
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   recvIt = recvDofs.begin();
214
215
	 recvIt != recvDofs.end();
	 ++recvIt, i++) {
216
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
217
	(*vec)[(recvIt->second)[j][0]] = recvBuffers[i][j];
218
219
220
221

      delete [] recvBuffers[i];
    }
    
222
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
223
224
225
      delete [] sendBuffers[i];

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

228

229
230
231
232
233
234
235
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
  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;
  }

263

264
265
266
267
268
269
270
271
272
273
274
275
276
277
  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);
  }

278

279
280
  void ParallelDomainProblemBase::createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
							     std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
281
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
282
283
284
285
    FUNCNAME("ParallelDomainProblemBase::createInteriorBoundaryInfo()");

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

286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
    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) {
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
	    // 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;

330
	    // === And add the part of the interior boundary. ===
331
332

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
335
336
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

337
338
339
340
341
342
	    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;
343
344
345
346
347
348
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
351
352
353
354
355
356
357
358
359
360
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
386
387
388
389
390
391
392
393
394
395
396
397


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

    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end();
	 ++rankIt) {
      int* buffer = new int[rankIt->second.size()];
      for (int i = 0; i < static_cast<int>(rankIt->second.size()); i++)
	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
      sendBuffers.push_back(buffer);
      
      mpiComm.Isend(buffer, rankIt->second.size(), MPI_INT, rankIt->first, 0);
    }

    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end();
	 ++rankIt) {
      int *buffer = new int[rankIt->second.size()];
      recvBuffers.push_back(buffer);
      
      mpiComm.Irecv(buffer, rankIt->second.size(), MPI_INT, rankIt->first, 0);      
    }

    mpiComm.Barrier();

    int i = 0;
    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = otherIntBoundary.boundary.begin();
	 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.
398
399
	  TEST_EXIT(k < static_cast<int>(rankIt->second.size()))
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
402
403
404
405
406
407
408
409
410
411
412

	  // 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];
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
  }


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


433
434
435
436
437
  void ParallelDomainProblemBase::createLocalGlobalNumbering(
                                 std::vector<const DegreeOfFreedom*>& rankDOFs,
				 std::map<const DegreeOfFreedom*, int>& boundaryDOFs,
				 int& nRankDOFs, 
				 int& nOverallDOFs)
438
  {
439
440
441
    FUNCNAME("ParallelDomainProblemBase::createLocalGlobalNumbering()");

    // === Get rank information about DOFs. ===
442
443
444
445
446
447
448
449
450

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

451
452

    // === Get starting position for global rank dof ordering. ====
453
454
455
456
457
458

    int rstart = 0;
    MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
    rstart -= nRankDOFs;

   
459
460
461
462
463
464
    // === 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;
465

466
467
468
469
    // 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;
470
471
472
473
474
475
476
477
478

    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
479
480
	DegreeOfFreedom newDofIndex = 0;
	for (int i = 0; i < nRankDOFs; i++) {
481
482
483
484
485
486
487
488
489
490
491
	  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)
492
	    sendNewDofs[*itRanks][it->first] = newDofIndex;
493
494
495
496
	}
      } 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.
497
498
499
500
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
501
502
503
504
      }
    }


505
    // === Send and receive the dof indices at boundary. ===
506
507
508

    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());
509
510
511

    sendDofs.clear();
    recvDofs.clear();
512
513
    
    int i = 0;
514
515
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
516
517
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
518
519
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
520
      int c = 0;
521
522
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
523
	   dofIt != sendIt->second.end();
524
525
526
	   ++dofIt) {
	sendBuffers[i][c++] = (dofIt->first)[0];
	sendBuffers[i][c++] = dofIt->second;
527

528
	sendDofs[sendIt->first].push_back(dofIt->first);
529
530
      }

531
      mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
532
533
534
    }

    i = 0;
535
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
536
537
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
538
539
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
540

541
      mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
542
543
544
545
546
547
    }


    mpiComm.Barrier();

    
548
    // === Delete send buffers. ===
549
550

    i = 0;
551
552
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
553
554
555
556
557
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) 
      delete [] sendBuffers[i];


558
559
560
561
562
    // === Change dof indices for rank partition. ===

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

564
    for (int i = 0; i < nRankDOFs; i++) {
565
566
567
568
569
570
571
      const_cast<DegreeOfFreedom*>(rankDOFs[i])[0] = i; 
      mapLocalGlobalDOFs[i] = rstart + i;
      mapGlobalLocalDOFs[rstart + i] = i;
      isRankDOF[i] = true;
    }


572
    // === Change dof indices at boundary from other ranks. ===
573
574

    i = 0;
575
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
576
577
578
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

579
      for (int j = 0; j < recvIt->second; j++) {
580

581
582
583
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
	DegreeOfFreedom newLocalDof = mapLocalGlobalDOFs.size();
584

585
	bool found = false;
586
587
588
589

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

591
	  if ((dofIt->first)[0] == oldDof) {
592
	    recvDofs[recvIt->first].push_back(dofIt->first);
593
	    const_cast<DegreeOfFreedom*>(dofIt->first)[0] = newLocalDof;
594
595
	    mapLocalGlobalDOFs[newLocalDof] = newGlobalDof;
	    mapGlobalLocalDOFs[newGlobalDof] = newLocalDof;
596
	    isRankDOF[newLocalDof] = false;
597
	    found = true;
598
599
600
	    break;
	  }
	}
601
602

	TEST_EXIT(found)("Should not happen!\n");
603
604
605
606
607
608
609
      }

      delete [] recvBuffers[i];
    }
  }


610
611
612
  void ParallelDomainProblemBase::updateLocalGlobalNumbering(int& nRankDOFs, 
							     int& nOverallDOFs)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
613
614
    FUNCNAME("ParallelDomainProblemBase::updateLocalGlobalNumbering()");

615
    std::set<const DegreeOfFreedom*> rankDOFs;
616
    std::map<const DegreeOfFreedom*, int> newBoundaryDOFs;
617
    
618
    // === Get all DOFs in ranks partition. ===
619
620
621
622
623
624
625
626
627
628
629
630
631

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

635
636
    std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDofs = sendDofs;
    std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDofs = recvDofs;
637

638
    for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
Thomas Witkowski's avatar
Thomas Witkowski committed
639
640
	   myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end();
641
	 ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
      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");
	}

665
666
667
668
669
670
671
672
	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
673
674
675
	std::vector<const DegreeOfFreedom*> boundDOFs;
	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		   boundDOFs);
676
677
	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
678
	  sendNewDofs[it->first].push_back(boundDOFs[i]);
679
	}
680
681
682
      }
    }    

683

684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
    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
720

721
722
723
724
725
726
	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;
727
	  recvNewDofs[it->first].push_back(boundDOFs[i]);
728
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
729
730
731
732
      }

    }

733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
    nRankDOFs = rankDOFs.size();

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

    int rstart = 0;
    MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
    rstart -= nRankDOFs;


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

    MPI_Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);


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

766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
    std::vector<int*> sendBuffers(sendNewDofs.size());
    std::vector<int*> recvBuffers(recvNewDofs.size());

    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;
781

782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
      mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
    }

    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];
	
      mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
    }

    mpiComm.Barrier();

    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) {
809
      
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
      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++];
826
827
    }

828
829
830
831
    // === Update list of dofs that must be communicated for solution exchange. ===

    sendDofs = sendNewDofs;
    recvDofs = recvNewDofs;
832
833
834
835
836
837
838
839
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
840
841
842
    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
843
	addAllDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
844
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
845
	addAllDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
846
847
848
849
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
850
	addAllDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
851
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
852
	addAllDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
853
854
855
856
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
857
	addAllDOFs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
858
	dofs.push_back(el->getFirstChild()->getDOF(2));
859
	addAllDOFs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
860
861
862
863
864
865
866
867
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
  }


868
  void ParallelDomainProblemBase::createDOFMemberInfo(
869
		       std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
870
871
872
		       std::vector<const DegreeOfFreedom*>& rankDOFs,
		       std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
  {
873
    // === Determine to each dof the set of partitions the dof belongs to. ===
874
875
876
877
878
879
880
881
882
883
884
885
886

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

887
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
888

889
    // iterate over all DOFs
890
891
892
    for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDOFs.begin();
	 it != partitionDOFs.end();
	 ++it) {
893
894

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

899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
	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;
	}
926

927
928
929
930
931
      }
    }
  }


932
933
934
  ParallelDomainProblemScal::ParallelDomainProblemScal(const std::string& name,
						       ProblemScal *problem,
						       ProblemInstatScal *problemInstat)
935
936
937
938
939
    : ParallelDomainProblemBase(name, 
				problem, 
				problemInstat, 
				problem->getFESpace(),
				problem->getRefinementManager()),
940
      probScal(problem)
941
942
943
  {
  }

944
945
946
947
948
949
950
  void ParallelDomainProblemScal::initParallelization(AdaptInfo *adaptInfo)
  {
    ParallelDomainProblemBase::initParallelization(adaptInfo);

    probScal->getSystemMatrix()->setIsRankDOF(isRankDOF);
  }

951
952
953
954
955
956
957
958
  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());

959
960
    solvePetscMatrix(probScal->getSolution());

961
962
963
964
965
966
967
968
969
970
//     if (toDo.isSet(SOLVE))
//       iterationIF->getProblem()->solve(adaptInfo, false);

//     if (toDo.isSet(SOLVE_RHS))
//       iterationIF->getProblem()->solve(adaptInfo, true);

//     if (toDo.isSet(ESTIMATE)) 
//       iterationIF->getProblem()->estimate(adaptInfo);


971
    return flag;
972
  }
973
974

}