ParallelDomainProblem.cc 28.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, 1);
92

93
//     updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
94

95
//     exit(0);
96

97
    // === Create petsc matrix. ===
98

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

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

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

  void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
114
  {}
115

116

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

138
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
139
	}
140

141
142
143
144

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

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

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
152
    }    
153
154
  }

155

156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
  void ParallelDomainProblemBase::solvePetscMatrix(DOFVector<double> *vec)
  {
    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
      delete [] sendBuffers[i];    
224
225
  }

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

261

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

276

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

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

284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
    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) {
299
300
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
	    // 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;

328
	    // === And add the part of the interior boundary. ===
329
330

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

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

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
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


    // === 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.
396
397
	  TEST_EXIT(k < static_cast<int>(rankIt->second.size()))
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
400
401
402
403
404
405
406
407
408
409
410

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


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


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

    // === Get rank information about DOFs. ===
440
441
442
443
444
445
446
447
448

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

449
450

    // === Get starting position for global rank dof ordering. ====
451
452
453
454
455
456

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

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

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

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


503
    // === Send and receive the dof indices at boundary. ===
504
505
506

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

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

526
	sendDofs[sendIt->first].push_back(dofIt->first);
527
528
      }

529
      mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_INT, sendIt->first, 0);
530
531
532
    }

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

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


    mpiComm.Barrier();

    
546
    // === Delete send buffers. ===
547
548

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


556
557
558
559
560
    // === Change dof indices for rank partition. ===

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

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


570
    // === Change dof indices at boundary from other ranks. ===
571
572

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

577
      for (int j = 0; j < recvIt->second; j++) {
578

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

583
	bool found = false;
584
585
586
587

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

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

	TEST_EXIT(found)("Should not happen!\n");
601
602
603
604
605
      }

      delete [] recvBuffers[i];
    }

606
607
#if 1
    // === Create local information from sendDofs and recvDofs
608

609
610
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   it = sendDofs.begin();
611
612
	 it != sendDofs.end();
	 ++it)
613
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
614
	   dofIt != it->second.end();
615
616
617
618
	   dofIt++) {
	//	std::cout << "AND SET S " << (*dofIt)[0] << " = " << mapGlobalLocalDOFs[(*dofIt)[0]] << std::endl;
	//	const_cast<DegreeOfFreedom*>(*dofIt)[0] = mapGlobalLocalDOFs[(*dofIt)[0]];
      }
619

620
621
    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
	   it = recvDofs.begin();
622
623
	 it != recvDofs.end();
	 ++it)
624
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
625
	   dofIt != it->second.end();
626
627
628
629
630
	   dofIt++) {
	//	std::cout << "AND SET R " << (*dofIt)[0] << " = " << mapGlobalLocalDOFs[(*dofIt)[0]] << std::endl;
	//	const_cast<DegreeOfFreedom*>(*dofIt)[0] = mapGlobalLocalDOFs[(*dofIt)[0]];
      }
#endif
631
632
633
  }


634
635
636
  void ParallelDomainProblemBase::updateLocalGlobalNumbering(int& nRankDOFs, 
							     int& nOverallDOFs)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
637
638
    FUNCNAME("ParallelDomainProblemBase::updateLocalGlobalNumbering()");

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

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

659
660
661
    std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDOFs;
    std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDOFs;

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

689
690
691
692
693
694
695
696
	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
697
698
699
	std::vector<const DegreeOfFreedom*> boundDOFs;
	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		   boundDOFs);
700
701
702
703
	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
	  sendNewDOFs[it->first].push_back(boundDOFs[i]);
	}
704
705
706
      }
    }    

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
    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
744

745
746
747
748
749
750
751
752
	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;
	  recvNewDOFs[it->first].push_back(boundDOFs[i]);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
753
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
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
    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. ===

    for (int rank = 0; rank < mpiSize; rank++) {
      if (rank == mpiRank) 
	continue;

      
    }

  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
805
806
807
    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
808
	addAllDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
809
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
810
	addAllDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
811
812
813
814
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
815
	addAllDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
816
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
817
	addAllDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
818
819
820
821
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
822
	addAllDOFs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
823
	dofs.push_back(el->getFirstChild()->getDOF(2));
824
	addAllDOFs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
825
826
827
828
829
830
831
832
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
  }


833
  void ParallelDomainProblemBase::createDOFMemberInfo(
834
		       std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
835
836
837
		       std::vector<const DegreeOfFreedom*>& rankDOFs,
		       std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
  {
838
    // === Determine to each dof the set of partitions the dof belongs to. ===
839
840
841
842
843
844
845
846
847
848
849
850
851

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

852
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
853

854
    // iterate over all DOFs
855
856
857
    for (std::map<const DegreeOfFreedom*, std::set<int> >::iterator it = partitionDOFs.begin();
	 it != partitionDOFs.end();
	 ++it) {
858
859

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

864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
	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;
	}
891

892
893
894
895
896
      }
    }
  }


897
898
899
  ParallelDomainProblemScal::ParallelDomainProblemScal(const std::string& name,
						       ProblemScal *problem,
						       ProblemInstatScal *problemInstat)
900
901
902
903
904
    : ParallelDomainProblemBase(name, 
				problem, 
				problemInstat, 
				problem->getFESpace(),
				problem->getRefinementManager()),
905
      probScal(problem)
906
907
908
  {
  }

909
910
911
912
913
914
915
  void ParallelDomainProblemScal::initParallelization(AdaptInfo *adaptInfo)
  {
    ParallelDomainProblemBase::initParallelization(adaptInfo);

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

916
917
918
919
920
921
922
923
  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());

924
925
    solvePetscMatrix(probScal->getSolution());

926
927
928
929
930
931
932
933
934
935
//     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);


936
    return flag;
937
  }
938
939

}