ParallelDomainBase.cc 58.1 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
#include <algorithm>

3
#include "ParallelDomainBase.h"
4
5
6
7
8
9
10
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
11
12
#include "DOFMatrix.h"
#include "DOFVector.h"
13
#include "SystemVector.h"
14
#include "VtkWriter.h"
15
#include "ElementDofIterator.h"
16
17
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
18
#include "ElementFileWriter.h"
19
20

#include "petscksp.h"
21
22
23

namespace AMDiS {

24
  PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *)
25
  {    
26
    if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
27
      std::cout << "  Petsc-Iteration " << iter << ": " << rnorm << std::endl;
28
29
30
31

    return 0;
  }

32
33
34
35
36
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

37
  ParallelDomainBase::ParallelDomainBase(ProblemIterationInterface *iIF,
38
39
40
					 ProblemTimeInterface *tIF,
					 FiniteElemSpace *fe,
					 RefinementManager *refineManager)
41
42
    : iterationIF(iIF),
      timeIF(tIF),
43
      name(iIF->getName()),
44
45
      feSpace(fe),
      mesh(fe->getMesh()),
46
      refinementManager(refineManager),
47
      initialPartitionMesh(true),
48
      nRankDofs(0),
49
50
      rstart(0),
      nComponents(1)
51
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
54
55
56
57
58
    FUNCNAME("ParallelDomainBase::ParalleDomainBase()");

    TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
      ("Only meshes with one DOFAdmin are supported!\n");
    TEST_EXIT(mesh->getDOFAdmin(0).getNumberOfPreDOFs(0) == 0)
      ("Wrong pre dof number for DOFAdmin!\n");

59
60
61
62
63
64
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
  }

65
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
66
67
68
69
  {
    if (mpiSize <= 1)
      return;

70
71
72
73
74
    // Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
    // only for macro meshes, so it will not work yet if the mesh is already refined
    // in some way.
    testForMacroMesh();

75
76
77
78
79
80
81
    // 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);   

Thomas Witkowski's avatar
Thomas Witkowski committed
82
83
84
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
    DbgCreateElementMap(elMap);
85
86
87

    if (mpiRank == 0)
      writePartitioningMesh("part.vtu");
Thomas Witkowski's avatar
Thomas Witkowski committed
88
#endif
89

90
    // === Create new global and local DOF numbering. ===
91

92
93
94
    // Set of all DOFs of the rank.
    std::vector<const DegreeOfFreedom*> rankDOFs;
    // Number of DOFs in ranks partition that are owned by the rank.
95
    nRankDofs = 0;
96
    // Number of all DOFs in the macro mesh.
97
    int nOverallDOFs = 0;
98

99
    createLocalGlobalNumbering(rankDOFs, nRankDofs, nOverallDOFs);
100

Thomas Witkowski's avatar
Thomas Witkowski committed
101
102
    // === Create interior boundary information ===

103
    createInteriorBoundaryInfo(rankDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
104

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
105
106
107
108
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

Thomas Witkowski's avatar
Thomas Witkowski committed
109
#if (DEBUG != 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
110
    DbgTestElementMap(elMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
111
112
    DbgTestInteriorBoundary();
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
113

114

115
    // === Reset all DOFAdmins of the mesh. ===
116

117
    updateDofAdmins();
118

119

120
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
121

Thomas Witkowski's avatar
Thomas Witkowski committed
122
123
124
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

Thomas Witkowski's avatar
Thomas Witkowski committed
125
126
    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);
127

128
129
130
131
132
#if (DEBUG != 0)
      elMap.clear();
      DbgCreateElementMap(elMap);
#endif

133
      updateLocalGlobalNumbering(nRankDofs, nOverallDOFs);
134
135
136
137
138
139

      updateDofAdmins();

#if (DEBUG != 0)
      DbgTestElementMap(elMap);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
    }

142

Thomas Witkowski's avatar
Thomas Witkowski committed
143
#if (DEBUG != 0)
144
    DbgTestCommonDofs(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
145
#endif
146

147
    nRankRows = nRankDofs * nComponents;
148
    nOverallRows = nOverallDOFs * nComponents;
149
150
  }

151
152

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
153
  {}
154

155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
  
  void ParallelDomainBase::updateDofAdmins()
  {
    int nAdmins = mesh->getNumberOfDOFAdmin();
    for (int i = 0; i < nAdmins; i++) {
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
      for (int j = 0; j < static_cast<int>(mapLocalGlobalDOFs.size()); j++)
 	admin.setDOFFree(j, false);

      admin.setUsedSize(mapLocalGlobalDOFs.size());
      admin.setUsedCount(mapLocalGlobalDOFs.size());
      admin.setFirstHole(mapLocalGlobalDOFs.size());
    }
  }

173

174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
  void ParallelDomainBase::testForMacroMesh()
  {
    FUNCNAME("ParallelDomainBase::testForMacroMesh()");

    int nMacroElements = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      TEST_EXIT(elInfo->getLevel() == 0)
	("Mesh is already refined! This does nor work with parallelization!\n");
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

    TEST_EXIT(nMacroElements >= mpiSize)
      ("The mesh has less macro elements than number of mpi processes!\n");
  }

195
196
  void ParallelDomainBase::setDofMatrix(DOFMatrix* mat, int dispMult, 
					int dispAddRow, int dispAddCol)
197
  {
198
199
200
201
    FUNCNAME("ParallelDomainBase::setDofMatrix()");

    TEST_EXIT(mat)("No DOFMatrix!\n");

202
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
203
204
205
206
207
208
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    traits::col<Matrix>::type col(mat->getBaseMatrix());
    traits::const_value<Matrix>::type value(mat->getBaseMatrix());

209
    typedef traits::range_generator<row, Matrix>::type cursor_type;
210
211
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

212
213
214
215
216
217
218
219
220
221
222
223
224
    std::vector<int> cols;
    std::vector<double> values;
    cols.reserve(300);
    values.reserve(300);

    for (cursor_type cursor = begin<row>(mat->getBaseMatrix()), 
	   cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {

      cols.clear();
      values.clear();

      int r = mapLocalGlobalDOFs[*cursor] * dispMult + dispAddRow;

225
226
      for (icursor_type icursor = begin<nz>(cursor), 
	     icend = end<nz>(cursor); icursor != icend; ++icursor)
227
	if (value(*icursor) != 0.0) {
228
229
 	  cols.push_back(mapLocalGlobalDOFs[col(*icursor)] * dispMult + dispAddCol);
 	  values.push_back(value(*icursor));
230
	}
231
232
233

      MatSetValues(petscMatrix, 1, &r, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES);
    }
234
  }
235

236

237
238
239
  void ParallelDomainBase::setDofVector(DOFVector<double>* vec, 
					int dispMult, int dispAdd)
  {
240
241
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
242
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()] * dispMult + dispAdd;
243
244
245
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
246
    }    
247
248
  }

249

250
251
  void ParallelDomainBase::fillPetscMatrix(DOFMatrix *mat, DOFVector<double> *vec)
  {
252
253
254
255
    MatCreate(PETSC_COMM_WORLD, &petscMatrix);
    MatSetSizes(petscMatrix, nRankRows, nRankRows, nOverallRows, nOverallRows);
    MatSetType(petscMatrix, MATAIJ);

256
257
258
259
260
261
262
263
264
265
266
267
    VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
    VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
    VecSetType(petscRhsVec, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &petscSolVec);
    VecSetSizes(petscSolVec, nRankRows, nOverallRows);
    VecSetType(petscSolVec, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
    VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
    VecSetType(petscTmpVec, VECMPI);

268
269
270
271
272
    setDofMatrix(mat);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
273
274
275
276
    setDofVector(vec);

    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);
277
278
279
280
281
  }

  
  void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
  {
282
283
284
285
    FUNCNAME("ParallelDomainBase::fillPetscMatrix()");

    clock_t first = clock();

286
287
288
289
290
291
292
293
294
295
296
297
    VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
    VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
    VecSetType(petscRhsVec, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &petscSolVec);
    VecSetSizes(petscSolVec, nRankRows, nOverallRows);
    VecSetType(petscSolVec, VECMPI);

    VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
    VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
    VecSetType(petscTmpVec, VECMPI);

298
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
299
300
301
302
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    int d_nnz[nRankRows];
303
304
    int o_nnz[nRankRows];

305
306
    std::map<int, std::vector< std::pair<int, int> > > sendMatrixEntry;

307
    for (int i = 0; i < nRankRows; i++) {
308
      d_nnz[i] = 0;
309
310
      o_nnz[i] = 0;
    }
311

312
313
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
314
315
316
317
318
319
 	if ((*mat)[i][j]) {
	  Matrix bmat = (*mat)[i][j]->getBaseMatrix();

	  traits::col<Matrix>::type col(bmat);
	  traits::const_value<Matrix>::type value(bmat);
	  
320
	  typedef traits::range_generator<row, Matrix>::type cursor_type;
321
322
	  typedef traits::range_generator<nz, cursor_type>::type icursor_type;
	  
323
324
325
326
327
	  for (cursor_type cursor = begin<row>(bmat), 
		 cend = end<row>(bmat); cursor != cend; ++cursor) {

	    int r = mapLocalGlobalDOFs[*cursor] * nComponents + i;

328
329
	    if (isRankDof[*cursor]) {
	      r -= rstart * nComponents;
330
331
332
333
334
335
336
337
338
339

#if (DEBUG != 0)    
	      if (r < 0 || r >= nRankRows) {
		std::cout << "ERROR in rank: " << mpiRank << std::endl;
		std::cout << "  Wrong r = " << r << " " << *cursor << " " 
			  << mapLocalGlobalDOFs[*cursor] << " " 
			  << nRankRows << std::endl;
		ERROR_EXIT("Should not happen!\n");
	      }
#endif
340
	      
341
342
343
344
	      for (icursor_type icursor = begin<nz>(cursor), 
		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
		if (value(*icursor) != 0.0) {
		  int c = mapLocalGlobalDOFs[col(*icursor)] * nComponents + j;
345

346
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
		  if (c >= rstart * nComponents && 
		      c < rstart * nComponents + nRankRows)
		    d_nnz[r]++;
		  else
		    o_nnz[r]++;		  
		}    
	      }
	    } else {
	      int sendToRank = -1;

	      for (RankToDofContainer::iterator it = recvDofs.begin();
		   it != recvDofs.end(); ++it) {
		for (DofContainer::iterator dofIt = it->second.begin();
		     dofIt != it->second.end(); ++dofIt) {
		  if (**dofIt == *cursor) {
		    sendToRank = it->first;
		    break;
		  }
		}

		if (sendToRank != -1)
		  break;
	      }

	      TEST_EXIT_DBG(sendToRank != -1)("Should not happen!\n");

	      for (icursor_type icursor = begin<nz>(cursor), 
		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
		if (value(*icursor) != 0.0) {
		  int c = mapLocalGlobalDOFs[col(*icursor)] * nComponents + j;
		  
		  sendMatrixEntry[sendToRank].push_back(std::make_pair(r, c));
		}
	      }
	    }
381
382
	  }
	}
383
384
385
386
387
388
389
390
391
392
393
394
      }
    }

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;

    std::vector<int*> sendBuffers;
    sendBuffers.reserve(recvDofs.size());

    for (RankToDofContainer::iterator it = recvDofs.begin(); 
	 it != recvDofs.end(); ++it) {
      int nSend = sendMatrixEntry[it->first].size();
395

396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
      request[requestCounter++] = mpiComm.Isend(&nSend, 1, MPI_INT, it->first, 0);
      
      if (nSend > 0) {
	int *sendbuf = new int[nSend * 2];
	for (int i = 0; i < nSend; i++) {
	  sendbuf[i * 2] = sendMatrixEntry[it->first][i].first;
	  sendbuf[i * 2 + 1] = sendMatrixEntry[it->first][i].second;
	}
	sendBuffers.push_back(sendbuf);
      } else {
	sendBuffers.push_back(NULL);
      }
    }

    std::vector<int> recvSize;
    recvSize.reserve(sendDofs.size());
    
    int i = 0;
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      request[requestCounter++] = 
	mpiComm.Irecv(&(recvSize[i++]), 1, MPI_INT, it->first, 0);

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

    requestCounter = 0;

    i = 0;
    for (RankToDofContainer::iterator it = recvDofs.begin(); 
	 it != recvDofs.end(); ++it) {
      int nSend = sendMatrixEntry[it->first].size();

428
      if (nSend > 0)
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
	request[requestCounter++] = 
	  mpiComm.Isend(sendBuffers[i], nSend * 2, MPI_INT, it->first, 0);

      i++;
    }

    std::vector<int*> recvBuffers;
    recvBuffers.reserve(sendDofs.size());
    
    i = 0;
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it) {
      if (recvSize[i] > 0) {
	int *recvbuf = new int[recvSize[i] * 2];
	recvBuffers.push_back(recvbuf);

	request[requestCounter++] =
	  mpiComm.Irecv(recvbuf, recvSize[i] * 2, MPI_INT, it->first, 0);
      } else {
	recvBuffers.push_back(NULL);
      }

      i++;
    }

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

456
457
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      if (sendBuffers[j])
Thomas Witkowski's avatar
Thomas Witkowski committed
458
 	delete [] sendBuffers[j];
459
460
461
462
463
464
465
466

    i = 0;
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it) {
      if (recvSize[i] > 0) {
	for (int j = 0; j < recvSize[i]; j++) {
	  int r = recvBuffers[i][j * 2];
	  int c = recvBuffers[i][j * 2 + 1];
467

468
	  r -= rstart * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
469

470
471
472
473
474
	  TEST_EXIT_DBG(r >= 0 && r < nRankRows)("Should not happen!\n");
	  
	  if (c < rstart * nComponents ||
	      c >= rstart * nComponents + nRankRows)
	    o_nnz[r]++;
475
476
	  else
	    d_nnz[r]++;
477
478
479
480
	}

	delete [] recvBuffers[i];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
481
482

      i++;
483
    }
484
485

    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
486
487
488
489
490
491
492
493
		    0, d_nnz, 0, o_nnz, &petscMatrix);

#if (DEBUG != 0)
    int a, b;
    MatGetOwnershipRange(petscMatrix, &a, &b);
    TEST_EXIT(a == rstart * nComponents)("Wrong matrix ownership range!\n");
    TEST_EXIT(b == rstart * nComponents + nRankRows)("Wrong matrix ownership range!\n");
#endif
494

495
496
497
498
499
500
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
501
	if ((*mat)[i][j])
502
503
504
505
506
507
508
	  setDofMatrix((*mat)[i][j], nComponents, i, j);

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

    for (int i = 0; i < nComponents; i++)
      setDofVector(vec->getDOFVector(i), nComponents, i);
509

510
511
512
    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);

513
    INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));
514
515
516
  }


517
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
518
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
519
520
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

521
522
523
524
    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
525
    KSPSetOperators(ksp, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN);
526
    KSPGetPC(ksp, &pc);
527
528
    //PCSetType(pc, PCJACOBI);
    PCSetType(pc, PCILU);
529
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
530
531
    KSPSetType(ksp, KSPBCGS);
    //KSPSetType(ksp, KSPCG);
532
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
533
534
535
536
537
    KSPSolve(ksp, petscRhsVec, petscSolVec);

    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

538
    for (int i = 0; i < nRankDofs; i++)
539
      (*vec)[mapLocalToDofIndex[i]] = vecPointer[i];
540
541
542

    VecRestoreArray(petscSolVec, &vecPointer);

543
544
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
545
546
547

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
548
549
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
550
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
551
	 sendIt != sendDofs.end(); ++sendIt, i++) {
552
553
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
554

555
      for (int j = 0; j < nSendDOFs; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
556
	sendBuffers[i][j] = (*vec)[*((sendIt->second)[j])];
557

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

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
563
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
564
	 recvIt != recvDofs.end(); ++recvIt, i++) {
565
566
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
567

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

Thomas Witkowski's avatar
Thomas Witkowski committed
572
573
574

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

575
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
576
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
577
	 recvIt != recvDofs.end(); ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
578
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
579
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
580
581
582
583

      delete [] recvBuffers[i];
    }
    
584
    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
585
      delete [] sendBuffers[i];
586
587

    MatDestroy(petscMatrix);
588
589
590
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
    VecDestroy(petscTmpVec);
591
592
  }

593

594
  void ParallelDomainBase::solvePetscMatrix(SystemVector &vec)
595
596
597
  {
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

598
    KSP solver;
599

600
601
602
    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 

603
604
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);

605
606
    KSPSetType(solver, KSPBCGS);
    KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
607
    KSPSetFromOptions(solver);
608
609

    KSPSolve(solver, petscRhsVec, petscSolVec);
610
611
612
613
614

    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

    for (int i = 0; i < nComponents; i++) {
615
      DOFVector<double> *dofvec = vec.getDOFVector(i);
616
      for (int j = 0; j < nRankDofs; j++)
617
	(*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];      
618
619
620
621
    }

    VecRestoreArray(petscSolVec, &vecPointer);

622
    synchVectors(vec);
623

624
625
626
627
628
629
630
631
632
633
634
    int iterations = 0;
    KSPGetIterationNumber(solver, &iterations);
    MSG("  Number of iterations: %d\n", iterations);
    
    double norm = 0.0;
    MatMult(petscMatrix, petscSolVec, petscTmpVec);
    VecAXPY(petscTmpVec, -1.0, petscRhsVec);
    VecNorm(petscTmpVec, NORM_2, &norm);
    MSG("  Residual norm: %e\n", norm);

    MatDestroy(petscMatrix);
635
636
637
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
    VecDestroy(petscTmpVec);
638
639
640
641
642
    KSPDestroy(solver);
  }
  
  void ParallelDomainBase::synchVectors(SystemVector &vec)
  {
643
644
645
646
647
648
649
650
651
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
    
    int i = 0;
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt, i++) {
652
653
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs * nComponents];
654
655
656

      int counter = 0;
      for (int j = 0; j < nComponents; j++) {
657
	DOFVector<double> *dofvec = vec.getDOFVector(j);
658
659
660
661
	for (int k = 0; k < nSendDOFs; k++)
	  sendBuffers[i][counter++] = (*dofvec)[*((sendIt->second)[k])];
      }

662
663
      request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs * nComponents,
						MPI_DOUBLE, sendIt->first, 0);
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
    }

    i = 0;
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt, i++) {
      int nRecvDOFs = recvIt->second.size() * nComponents;
      recvBuffers[i] = new double[nRecvDOFs];

      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
    }


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

    i = 0;
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt, i++) {
      int nRecvDOFs = recvIt->second.size();

      int counter = 0;
      for (int j = 0; j < nComponents; j++) {
686
	DOFVector<double> *dofvec = vec.getDOFVector(j);
687
688
 	for (int k = 0; k < nRecvDOFs; k++)
 	  (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++];
689
690
691
692
693
694
695
696
697
      }

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

698
699
700
701
702
703
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
744
745
746
747
748
749
750
751
752
753
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
  {
    int vecSize = data.size();
    SerUtil::serialize(out, &vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
      SerUtil::serialize(out, &dofIndex);
    }
  }


  void ParallelDomainBase::deserialize(std::istream &in, DofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    FUNCNAME("ParallelDomainBase::deserialize()");

    int vecSize = 0;
    SerUtil::deserialize(in, &vecSize);
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
      SerUtil::deserialize(in, &dofIndex);

      TEST_EXIT_DBG(dofMap.count(dofIndex) != 0)
	("Dof index could not be deserialized correctly!\n");

      data[i] = dofMap[dofIndex];
    }
  }


  void ParallelDomainBase::serialize(std::ostream &out, RankToDofContainer &data)
  {
    int mapSize = data.size();
    SerUtil::serialize(out, &mapSize);
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
      SerUtil::serialize(out, &rank);
      serialize(out, it->second);
    }
  }

  
  void ParallelDomainBase::deserialize(std::istream &in, RankToDofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    int mapSize = 0;
    SerUtil::deserialize(in, &mapSize);
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
      SerUtil::deserialize(in, &rank);
      deserialize(in, data[rank], dofMap);      
    }
  }

754

755
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
  {
    double localWeightSum = 0.0;
    int elNum = -1;

    elemWeights.clear();

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

      // get partition data
      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(element->getElementData(PARTITION_ED));

      if (partitionData && partitionData->getPartitionStatus() == IN) {
	if (partitionData->getLevel() == 0) {
	  elNum = element->getIndex();
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
776
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
777
778
779
780
781
782
783
784
785
786
787
788
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

789

790
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
791
792
793
794
795
796
797
798
799
800
801
802
803
  {
    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);
  }

804

805
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs)
806
  {
807
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
808
809
810

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

811
    TraverseStack stack;
812
813
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
814
815
816
817
818
    while (elInfo) {
      Element *element = elInfo->getElement();

      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
819

820
821
822
823
824
825
826
      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));
827

828
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
829
830
	    // We have found an element that is at an interior boundary. 

831
832
	    int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];
	    bool ranksBoundary = (mpiRank > otherElementRank);
833
834

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
835
	      (ranksBoundary ?
836
837
	       myIntBoundary.getNewAtomicBoundary(otherElementRank) :
	       otherIntBoundary.getNewAtomicBoundary(otherElementRank));
Thomas Witkowski's avatar
Thomas Witkowski committed
838

839
	    bound.rankObject.el = element;
840
	    bound.rankObject.elIndex = element->getIndex();
841
842
	    bound.rankObject.subObjAtBoundary = EDGE;
	    bound.rankObject.ithObjAtBoundary = i;
843
844
845
846
	    // Do not set a pointer to the element, because if will be deleted from
	    // mesh after partitioning and the pointer would become invalid.
	    bound.neighbourObject.el = NULL;
	    bound.neighbourObject.elIndex = elInfo->getNeighbour(i)->getIndex();
847
	    bound.neighbourObject.subObjAtBoundary = EDGE;
848
849
850
851
852
	    bound.neighbourObject.ithObjAtBoundary = elInfo->getSideOfNeighbour(i);

	    // i == 2  =>  getSideOfNeighbour(i) == 2
	    TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
	      ("Should not happen!\n");
853
854
855
856
857
858
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
859
860
861


    // === Once we have this information, we must care about their order. Eventually ===
Thomas Witkowski's avatar
Thomas Witkowski committed
862
863
    // === all the boundaries have to be in the same order on both ranks that share  ===
    // === the bounday.                                                              ===
Thomas Witkowski's avatar
Thomas Witkowski committed
864

Thomas Witkowski's avatar
Thomas Witkowski committed
865
    std::vector<int*> sendBuffers, recvBuffers;
Thomas Witkowski's avatar
Thomas Witkowski committed
866

Thomas Witkowski's avatar
Thomas Witkowski committed
867
868
869
870
    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
871
872
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
873
      int nSendInt = rankIt->second.size();
874
875
876
877
878
      int* buffer = new int[nSendInt * 2];
      for (int i = 0; i < nSendInt; i++) {
	buffer[i * 2] = (rankIt->second)[i].rankObject.elIndex;
	buffer[i * 2 + 1] = (rankIt->second)[i].rankObject.ithObjAtBoundary;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
879
880
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
881
      request[requestCounter++] =
882
	mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
883
884
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
885
886
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
887
      int nRecvInt = rankIt->second.size();
888
      int *buffer = new int[nRecvInt * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
889
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
890

Thomas Witkowski's avatar
Thomas Witkowski committed
891
      request[requestCounter++] = 
892
	mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
893
894
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
895
896
897

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

Thomas Witkowski's avatar
Thomas Witkowski committed
898
899

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
900
901
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
902
903
904
905
906
907
908

      // === 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.
909
910
	if ((rankIt->second)[j].neighbourObject.elIndex != recvBuffers[i][j * 2] ||
	    (rankIt->second)[j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
911
912
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
913
914
	    if ((rankIt->second)[k].neighbourObject.elIndex == recvBuffers[i][j * 2] && 
		(rankIt->second)[k].neighbourObject.ithObjAtBoundary == recvBuffers[i][j * 2 + 1])
Thomas Witkowski's avatar
Thomas Witkowski committed
915
916
917
	      break;

	  // The element must always be found, because the list is just in another order.
Thomas Witkowski's avatar
Thomas Witkowski committed
918
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
919
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
920
921
922
923
924
925
926
927
928
929
930
931
932

	  // 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];
933
934
935
  }


936
  void ParallelDomainBase::removeMacroElements()
937
938
939
940
941
942
943
944
945
946
947
948
  {
    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);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
949
    mesh->removeMacroElements(macrosToRemove, feSpace);
950
951
952
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
953
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
954
						      int& nRankDofs, 
Thomas Witkowski's avatar
Thomas Witkowski committed
955
						      int& nOverallDOFs)
956
  {
957
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
958
959

    // === Get rank information about DOFs. ===
960
961
962

    // Stores to each DOF pointer the set of ranks the DOF is part of.
    std::map<const DegreeOfFreedom*, std::set<int> > partitionDOFs;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
963
    DofContainer rankAllDofs;
964
    DofToRank boundaryDofs;
965

966
    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDofs, vertexDof);
967

968
    nRankDofs = rankDOFs.size();
969
970
    nOverallDOFs = partitionDOFs.size();

971

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

974
    rstart = 0;
975
976
    mpiComm.Scan(&nRankDofs, &rstart, 1, MPI_INT, MPI_SUM);
    rstart -= nRankDofs;
977

978
979
980
981
    
    // === Create for all dofs in rank new indices. The new index must start at ===
    // === index 0, must be continues and have the same order as the indices    ===
    // === had before.                                                          ===
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
982

Thomas Witkowski's avatar
Thomas Witkowski committed
983

984
985
986
    // Do not change the indices now, but create a new indexing a store it here.
    DofIndexMap rankDofsNewLocalIndex;
    isRankDof.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
987
988
989
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
990
991
992
      rankDofsNewLocalIndex[*dofIt] = i;
      // First, we set all dofs in ranks partition to be owend by the rank. Later,
      // the dofs in ranks partition that are owned by other rank are set to false.
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
993
      isRankDof[i] = true;
994
      i++;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
995
996
    }

997
998
999

    // === Create for all rank owned dofs a new global indexing.                ===

1000
1001
    // Stores for dofs in rank a new global index.
    DofIndexMap rankDofsNewGlobalIndex;
1002
1003
1004
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1005
1006
1007
    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
1008
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
1009
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1010
1011
1012
      i++;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1013
 
1014
1015
1016
1017
1018
    // === 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.
1019
    std::map<int, DofIndexMap > sendNewDofs;
1020

1021
1022
1023
1024
    // 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;
1025

1026
    for (DofToRank::iterator it = boundaryDofs.begin(); it != boundaryDofs.end(); ++it) {
1027
1028
1029
1030
1031
1032
1033
1034

      if (it->second == mpiRank) {
	// If the boundary dof is a rank dof, it must be send to other ranks.

	// 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) {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1035
	  if (*itRanks != mpiRank) {
1036
	    TEST_EXIT_DBG(rankDofsNewGlobalIndex.count(it->first) == 1)
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1037
1038
	      ("DOF Key not found!\n");

1039
	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1040
	  }
1041
1042
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1043
1044
	// 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.
1045
1046
1047
1048
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
1049
1050
1051
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1052

1053
    // === Send and receive the dof indices at boundary. ===
1054

Thomas Witkowski's avatar
Thomas Witkowski committed
1055
    std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size());
1056
1057
1058

    sendDofs.clear();
    recvDofs.clear();
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074

#if 0
    StlMpi stlMpi(mpiComm);
    stlMpi.prepareCommunication(false);
    std::map<

    for (std::map<int, DofIndexMap >::iterator sendIt = sendNewDofs.begin();
	 sendIt != sendNewDofs.end(); ++sendIt, i++)
      stlMpi.send(it->second, it->first);

    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end(); ++recvIt, i++)
      stlMpi.recv(it->first, it->second * 2); 

    stlMpi.startCommunication();
#endif
1075
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1076
1077
1078
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1079
    i = 0;
1080
1081
    for (std::map<int, DofIndexMap >::iterator sendIt = sendNewDofs.begin();
	 sendIt != sendNewDofs.end(); ++sendIt, i++) {
1082
1083
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
1084

1085
      int c = 0;
1086
1087
      for (DofIndexMap::iterator dofIt = sendIt->second.begin();
	   dofIt != sendIt->second.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1088
	sendBuffers[i][c++] = *(dofIt->first);
1089
	sendBuffers[i][c++] = dofIt->second;
1090

1091
	sendDofs[sendIt->first].push_back(dofIt->first);
1092
1093
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
1094
1095
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
1096
1097
1098
    }

    i = 0;
1099
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
1100
	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
1101
1102
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
1103

Thomas Witkowski's avatar
Thomas Witkowski committed
1104
1105
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
1106
1107
1108
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
1109
    MPI::Request::Waitall(requestCounter, request);
1110
1111

    
1112
    // === Delete send buffers. ===
1113

Thomas Witkowski's avatar
Thomas Witkowski committed
1114
1115
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
1116
1117


1118
    // === Change dof indices at boundary from other ranks. ===
1119

Thomas Witkowski's avatar
Thomas Witkowski committed
1120
1121
    // Within this small data structure we track which dof index was already changed.
    // This is used to avoid the following situation: Assume, there are two dof indices
1122
    // a and b in boundaryDofs. Then we have to change index a to b and b to c. When
Thomas Witkowski's avatar
Thomas Witkowski committed
1123
1124
1125
1126
1127
    // the second rule applies, we have to avoid that not the first b, resulted from
    // changing a to b, is set to c, but the second one. Therefore, after the first
    // rule was applied, the dof pointer is set to false in this data structure and 
    // is not allowed to be changed anymore.
    std::map<const DegreeOfFreedom*, bool> dofChanged;