ParallelDomainBase.cc 68.7 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
#include "VertexVector.h"
20
21

#include "petscksp.h"
22
23
24

namespace AMDiS {

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

    return 0;
  }

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

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

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

67
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
68
  {
69
70
71
72
73
74
75
    FUNCNAME("ParallelDomainBase::initParallelization()");

    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");

    // If the problem has been already read from a file, we do not need to do anything.
    if (deserialized)
76
77
      return;

78
79
80
81
82
    // 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();

83
84
85
86
87
88
89
    // 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
90
91
92
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
    DbgCreateElementMap(elMap);
93
94
95

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

98
    // === Create new global and local DOF numbering. ===
99

100
    // Set of all DOFs of the rank.
101
    std::vector<const DegreeOfFreedom*> rankDofs;
102
    // Number of DOFs in ranks partition that are owned by the rank.
103
    nRankDofs = 0;
104
    // Number of all DOFs in the macro mesh.
105
    int nOverallDOFs = 0;
106

107
    createLocalGlobalNumbering(rankDofs, nRankDofs, nOverallDOFs);
108

Thomas Witkowski's avatar
Thomas Witkowski committed
109
110
    // === Create interior boundary information ===

111
    createInteriorBoundaryInfo(rankDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
112

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
113
114
115
116
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

Thomas Witkowski's avatar
Thomas Witkowski committed
117
#if (DEBUG != 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
118
    DbgTestElementMap(elMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
119
120
    DbgTestInteriorBoundary();
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
121

122
    // === Reset all DOFAdmins of the mesh. ===
123

124
    updateDofAdmins();
125

126
127
128
    // === ===

    createPeriodicMap();
129

130
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
131

Thomas Witkowski's avatar
Thomas Witkowski committed
132
133
134
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

Thomas Witkowski's avatar
Thomas Witkowski committed
135
136
    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);
137

138
139
140
141
142
#if (DEBUG != 0)
      elMap.clear();
      DbgCreateElementMap(elMap);
#endif

143
      updateLocalGlobalNumbering(nRankDofs, nOverallDOFs);
144
145
146
147
148
149

      updateDofAdmins();

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

152

Thomas Witkowski's avatar
Thomas Witkowski committed
153
#if (DEBUG != 0)
154
    //    DbgTestCommonDofs(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
155
#endif
156

157
    nRankRows = nRankDofs * nComponents;
158
    nOverallRows = nOverallDOFs * nComponents;
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191

#if 0
    DOFVector<double> *tmp = new DOFVector<double>(feSpace, "tmp");
    tmp->set(0.0);
    if (mpiRank == 0)
      VtkWriter::writeFile(tmp, "test0.vtu");
    else
      VtkWriter::writeFile(tmp, "test1.vtu");

    if (mpiRank == 1) {
      for (PeriodicDofMap::iterator it = periodicDof.begin();
	   it != periodicDof.end(); ++it) {
	std::cout << "DOF MAP " << it->first << ": ";
	for (std::set<DegreeOfFreedom>::iterator dofit = it->second.begin();
	     dofit != it->second.end(); ++dofit)
	  std::cout << *dofit << " ";
	std::cout << std::endl;

	DegreeOfFreedom localdof;
	for (DofMapping::iterator dofIt = mapLocalGlobalDOFs.begin();
	     dofIt != mapLocalGlobalDOFs.end(); ++dofIt)
	  if (dofIt->second == it->first)
	    localdof = dofIt->first;

	WorldVector<double> coords;
	mesh->getDofIndexCoords(localdof, feSpace, coords);
	coords.print();
      }
    }
#endif

    if (mpiRank == 0)
      writePartitioningMesh("part.vtu");
192
193
  }

194
195

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
196
  {}
197

198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
  
  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());
    }
  }

216

217
218
219
220
221
222
223
224
225
226
  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)
227
	("Mesh is already refined! This does not work with parallelization!\n");
228
229
230
231
232
233
234
235
236
237
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

238
239
  void ParallelDomainBase::setDofMatrix(DOFMatrix* mat, int dispMult, 
					int dispAddRow, int dispAddCol)
240
  {
241
242
243
244
    FUNCNAME("ParallelDomainBase::setDofMatrix()");

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

245
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
246
247
248
249
250
251
    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());

252
    typedef traits::range_generator<row, Matrix>::type cursor_type;
253
254
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

255
256
257
258
259
260
261
262
263
264
265
    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();

266
267
      int globalRowDof = mapLocalGlobalDOFs[*cursor];
      bool periodicRow = (periodicDof.count(globalRowDof) > 0);
268

269
270
      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	   icursor != icend; ++icursor) {
271
	if (value(*icursor) != 0.0) {
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
	  int globalColDof = mapLocalGlobalDOFs[col(*icursor)];
	  int colIndex = globalColDof * dispMult + dispAddCol;

	  if (periodicRow == false && periodicDof.count(globalColDof) > 0) {
	    cols.push_back(colIndex);
	    values.push_back(value(*icursor) * 0.5);

	    for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalColDof].begin();
		 it != periodicDof[globalColDof].end(); ++it) {
	      cols.push_back(*it * dispMult + dispAddCol);
	      values.push_back(value(*icursor) * 0.5);
	    }
	  } else {
	    cols.push_back(colIndex);
	    values.push_back(value(*icursor));
	  }
288
	}
289
290
291
292
293
294
      }

      int rowIndex = globalRowDof * dispMult + dispAddRow;

      if (periodicRow) {
	int diagIndex = -1;
295

296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
	for (int i = 0; i < static_cast<int>(values.size()); i++) {
	  if (cols[i] != rowIndex)
	    values[i] *= 0.5;
	  else
	    diagIndex = i;
	}

	MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES);	

	if (diagIndex != -1)
	  values[diagIndex] = 0.0;

	for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRowDof].begin();
	     it != periodicDof[globalRowDof].end(); ++it) {
	  int perRowIndex = *it * dispMult + dispAddRow;
	  MatSetValues(petscMatrix, 1, &perRowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES);
	}

      } else {
	MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES);
      }
317
    }
318
  }
319

320

321
  void ParallelDomainBase::setDofVector(Vec& petscVec, DOFVector<double>* vec, 
322
323
					int dispMult, int dispAdd)
  {
324
325
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
326
327
328
329
330
331
332
333
334
335
336
337
      int globalRow = mapLocalGlobalDOFs[dofIt.getDOFIndex()];
      int index = globalRow * dispMult + dispAdd;

      if (periodicDof.count(globalRow) > 0) {
	double value = *dofIt * 0.5;
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);

	for (std::set<DegreeOfFreedom>::iterator it = periodicDof[globalRow].begin();
	     it != periodicDof[globalRow].end(); ++it) {
	  index = *it * dispMult + dispAdd;
	  VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
	}
338

339
340
341
342
      } else {
	double value = *dofIt;
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
      }
343
    }    
344
345
  }

346

347
348
  void ParallelDomainBase::fillPetscMatrix(DOFMatrix *mat, DOFVector<double> *vec)
  {
349
350
351
352
    MatCreate(PETSC_COMM_WORLD, &petscMatrix);
    MatSetSizes(petscMatrix, nRankRows, nRankRows, nOverallRows, nOverallRows);
    MatSetType(petscMatrix, MATAIJ);

353
354
355
356
357
358
359
360
361
362
363
364
    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);

365
366
367
368
369
    setDofMatrix(mat);

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

370
    setDofVector(petscRhsVec, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
373

    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);
374
375
376
377
378
  }

  
  void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
  {
379
380
381
382
    FUNCNAME("ParallelDomainBase::fillPetscMatrix()");

    clock_t first = clock();

383
384
385
386
387
388
389
390
391
392
393
394
    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);

395
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
396
397
398
399
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    int d_nnz[nRankRows];
400
401
    int o_nnz[nRankRows];

402
403
    std::map<int, std::vector< std::pair<int, int> > > sendMatrixEntry;

404
    for (int i = 0; i < nRankRows; i++) {
405
      d_nnz[i] = 0;
406
407
      o_nnz[i] = 0;
    }
408

409
410
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
411
412
413
414
415
416
 	if ((*mat)[i][j]) {
	  Matrix bmat = (*mat)[i][j]->getBaseMatrix();

	  traits::col<Matrix>::type col(bmat);
	  traits::const_value<Matrix>::type value(bmat);
	  
417
	  typedef traits::range_generator<row, Matrix>::type cursor_type;
418
419
	  typedef traits::range_generator<nz, cursor_type>::type icursor_type;
	  
420
421
422
423
424
	  for (cursor_type cursor = begin<row>(bmat), 
		 cend = end<row>(bmat); cursor != cend; ++cursor) {

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

425
426
	    if (isRankDof[*cursor]) {
	      r -= rstart * nComponents;
427
428
429
430
431
432
433
434
435
436

#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
437
	      
438
439
440
441
	      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;
442

443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
		  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));
		}
	      }
	    }
478
479
	  }
	}
480
481
482
483
484
485
486
487
488
489
490
491
      }
    }

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

493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
      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();

525
      if (nSend > 0)
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
	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
552

553
554
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      if (sendBuffers[j])
Thomas Witkowski's avatar
Thomas Witkowski committed
555
 	delete [] sendBuffers[j];
556
557
558
559
560
561
562
563

    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];
564

565
	  r -= rstart * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
566

567
568
	  TEST_EXIT_DBG(r >= 0 && r < nRankRows)("Should not happen!\n");
	  
569
	  if (c < rstart * nComponents || c >= rstart * nComponents + nRankRows)
570
	    o_nnz[r]++;
571
572
	  else
	    d_nnz[r]++;
573
574
575
576
	}

	delete [] recvBuffers[i];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
577
578

      i++;
579
    }
580
581

    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
582
583
584
585
586
587
588
589
		    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
590

591
592
593
594
595
596
    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++)
597
	if ((*mat)[i][j])
598
599
600
601
602
603
	  setDofMatrix((*mat)[i][j], nComponents, i, j);

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

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

606
607
608
    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);

609
    INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));
610
611
612
  }


613
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
614
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
615
616
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

617
618
619
620
    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
621
    KSPSetOperators(ksp, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN);
622
    KSPGetPC(ksp, &pc);
623
624
    //PCSetType(pc, PCJACOBI);
    PCSetType(pc, PCILU);
625
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
626
627
    KSPSetType(ksp, KSPBCGS);
    //KSPSetType(ksp, KSPCG);
628
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
629
630
631
632
633
    KSPSolve(ksp, petscRhsVec, petscSolVec);

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

634
    for (int i = 0; i < nRankDofs; i++)
635
      (*vec)[mapLocalToDofIndex[i]] = vecPointer[i];
636
637
638

    VecRestoreArray(petscSolVec, &vecPointer);

639
640
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
641
642
643

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
644
645
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
646
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
647
	 sendIt != sendDofs.end(); ++sendIt, i++) {
648
649
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
650

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

Thomas Witkowski's avatar
Thomas Witkowski committed
654
655
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
656
657
658
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
659
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
660
	 recvIt != recvDofs.end(); ++recvIt, i++) {
661
662
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
663

Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
666
667
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
668
669
670

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

671
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
672
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
673
	 recvIt != recvDofs.end(); ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
674
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
675
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
676
677
678
679

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

    MatDestroy(petscMatrix);
684
685
686
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
    VecDestroy(petscTmpVec);
687
688
  }

689

690
  void ParallelDomainBase::solvePetscMatrix(SystemVector &vec)
691
692
693
  {
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

694
695
696
697
698
699
700
701
702
#if 0
    // Set old solution to be initiual guess for petsc solver.
    for (int i = 0; i < nComponents; i++)
      setDofVector(petscSolVec, vec->getDOFVector(i), nComponents, i);

    VecAssemblyBegin(petscSolVec);
    VecAssemblyEnd(petscSolVec);
#endif

703
    KSP solver;
704

705
706
707
    KSPCreate(PETSC_COMM_WORLD, &solver);
    KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 

708
709
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);

710
711
    KSPSetType(solver, KSPBCGS);
    KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0);
712
    KSPSetFromOptions(solver);
713
714
    // Do not delete the solution vector, use it for the initial guess.
    //    KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
715
716

    KSPSolve(solver, petscRhsVec, petscSolVec);
717
718
719
720
721

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

    for (int i = 0; i < nComponents; i++) {
722
      DOFVector<double> *dofvec = vec.getDOFVector(i);
723
      for (int j = 0; j < nRankDofs; j++)
724
	(*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];      
725
726
727
728
    }

    VecRestoreArray(petscSolVec, &vecPointer);

729
    synchVectors(vec);
730

731
732
733
734
735
736
737
738
739
740
741
    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);
742
743
744
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
    VecDestroy(petscTmpVec);
745
746
747
748
749
    KSPDestroy(solver);
  }
  
  void ParallelDomainBase::synchVectors(SystemVector &vec)
  {
750
751
752
753
754
755
756
757
758
    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++) {
759
760
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs * nComponents];
761
762
763

      int counter = 0;
      for (int j = 0; j < nComponents; j++) {
764
	DOFVector<double> *dofvec = vec.getDOFVector(j);
765
766
767
768
	for (int k = 0; k < nSendDOFs; k++)
	  sendBuffers[i][counter++] = (*dofvec)[*((sendIt->second)[k])];
      }

769
770
      request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs * nComponents,
						MPI_DOUBLE, sendIt->first, 0);
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
    }

    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++) {
793
	DOFVector<double> *dofvec = vec.getDOFVector(j);
794
795
 	for (int k = 0; k < nRecvDOFs; k++)
 	  (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++];
796
797
798
799
800
801
802
803
804
      }

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

805
806
807
808
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
  {
    int vecSize = data.size();
809
    SerUtil::serialize(out, vecSize);
810
811
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
812
      SerUtil::serialize(out, dofIndex);
813
814
815
816
817
818
819
820
821
822
    }
  }


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

    int vecSize = 0;
823
    SerUtil::deserialize(in, vecSize);
824
825
826
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
827
      SerUtil::deserialize(in, dofIndex);
828
829
830
831
832
833
834
835
836
837
838
839

      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();
840
    SerUtil::serialize(out, mapSize);
841
842
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
843
      SerUtil::serialize(out, rank);
844
845
846
847
848
849
850
851
852
      serialize(out, it->second);
    }
  }

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

861

862
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
  {
    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) {
880
	if (partitionData->getLevel() == 0)
881
	  elNum = element->getIndex();
882
	
Thomas Witkowski's avatar
Thomas Witkowski committed
883
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
884
885
886
887
888
889
890
891
892
893
894
895
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

896

897
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
898
899
900
901
902
903
904
905
906
907
908
909
910
  {
    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);
  }

911

912
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDofs)
913
  {
914
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
915

916
917
918
919
920
    // To be general, we have to use elInfo->getBoundary(EDGE/FACE, i) instead of
    // elInfo->getBoundar(i).
    TEST_EXIT(mesh->getDim() == 2)
      ("getBoundary(i) returns always i-th edge, generalize for 3d!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
921
922
    // === First, create all the information about the interior boundaries. ===

923
    TraverseStack stack;
924
925
926
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
927
928
929
930
931
    while (elInfo) {
      Element *element = elInfo->getElement();

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

933
934
935
936
937
938
939
      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));
940

941
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
Thomas Witkowski's avatar
Thomas Witkowski committed
942

943
	    AtomicBoundary bound;	    	    
944
	    bound.rankObject.el = element;
945
	    bound.rankObject.elIndex = element->getIndex();
946
947
	    bound.rankObject.subObjAtBoundary = EDGE;
	    bound.rankObject.ithObjAtBoundary = i;
948
949
950
951
	    // 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();
952
	    bound.neighbourObject.subObjAtBoundary = EDGE;
953
	    bound.neighbourObject.ithObjAtBoundary = elInfo->getSideOfNeighbour(i);
954
	    
955
956
957
	    // i == 2  =>  getSideOfNeighbour(i) == 2
	    TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
	      ("Should not happen!\n");
958
959
960
961
962
963
964
965
966
967
968
969
970
971

	    int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];

	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(i))) {	      
	      // We have found an element that is at an interior, periodic boundary.
	      AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
	      b = bound;
	    } else {
	      // We have found an element that is at an interior, non-periodic boundary.
	      AtomicBoundary& b = ((mpiRank > otherElementRank) ?
				   myIntBoundary.getNewAtomic(otherElementRank) :
				   otherIntBoundary.getNewAtomic(otherElementRank));
	      b = bound;	      
	    }
972
973
974
975
976
977
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
978
979

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

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

985
986
987
    MPI::Request request[max(myIntBoundary.boundary.size() + 
			     otherIntBoundary.boundary.size(),
			     periodicBoundary.boundary.size())];
Thomas Witkowski's avatar
Thomas Witkowski committed
988
989
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
990
991
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
992
      int nSendInt = rankIt->second.size();
993
994
995
996
997
      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
998
999
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
1000
      request[requestCounter++] =
1001
	mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1002
1003
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1004
1005
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1006
      int nRecvInt = rankIt->second.size();
1007
      int *buffer = new int[nRecvInt * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
1008
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
1009

Thomas Witkowski's avatar
Thomas Witkowski committed
1010
      request[requestCounter++] = 
1011
	mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1012
1013
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1014
1015
1016

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1017
1018

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1019
1020
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1021
1022
1023
1024
1025
1026
1027

      // === 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.
1028
1029
	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
1030
1031
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
1032
1033
	    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
1034
1035
1036
	      break;

	  // The element must always be found, because the list is just in another order.
Thomas Witkowski's avatar
Thomas Witkowski committed
1037
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
1038
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
1039
1040
1041
1042
1043
1044
1045
1046

	  // Swap the current with the found element.
	  AtomicBoundary tmpBound = (rankIt->second)[k];
	  (rankIt->second)[k] = (rankIt->second)[j];
	  (rankIt->second)[j] = tmpBound;	
	}
      }

1047
1048
      delete [] recvBuffers[i];
      i++;
Thomas Witkowski's avatar
Thomas Witkowski committed
1049
1050
1051
1052
    }

    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124

    recvBuffers.clear();
    sendBuffers.clear();
 

    // === Deal with periodic boundaries. ===

    requestCounter = 0;

    for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
	 rankIt != periodicBoundary.boundary.end(); ++rankIt) {

      TEST_EXIT_DBG(rankIt->first != mpiRank)
	("It is no possible to have an interior boundary within a rank partition!\n");

      if (rankIt->first < mpiRank) {
	int nSendInt = rankIt->second.size();
	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;
	}
	sendBuffers.push_back(buffer);

	request[requestCounter++] =
	  mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
      } else {
	int nRecvInt = rankIt->second.size();
	int *buffer = new int[nRecvInt * 2];
	recvBuffers.push_back(buffer);
	
	request[requestCounter++] = 
	  mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);	
      }
    }


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


    i = 0;
    for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
	 rankIt != periodicBoundary.boundary.end(); ++rankIt) {
      if (rankIt->first > mpiRank) {

	for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) 
	  if (periodicBoundary.boundary[rankIt->first][j].neighbourObject.elIndex != recvBuffers[i][j * 2] ||
	      periodicBoundary.boundary[rankIt->first][j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) {
	    int k = j + 1;

	    for (; k < static_cast<int>(rankIt->second.size()); k++)
	      if (periodicBoundary.boundary[rankIt->first][k].neighbourObject.elIndex == recvBuffers[i][j * 2] && 
		  periodicBoundary.boundary[rankIt->first][k].neighbourObject.ithObjAtBoundary == recvBuffers[i][j * 2 + 1])
		break;

	    // The element must always be found, because the list is just in another order.
	    TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
	      ("Should never happen!\n");
	    
	    // 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];
	i++;
      }
    }

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


1128
  void ParallelDomainBase::removeMacroElements()
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
  {
    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
1141
    mesh->removeMacroElements(macrosToRemove, feSpace);
1142
1143
1144
  }


1145
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDofs,
1146
						      int& nRankDofs, 
Thomas Witkowski's avatar
Thomas Witkowski committed
1147
						      int& nOverallDOFs)
1148
  {
1149
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
1150
1151

    // === Get rank information about DOFs. ===
1152
1153
1154

    // 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
1155
    DofContainer rankAllDofs;
1156
    DofToRank boundaryDofs;
1157

1158
    createDofMemberInfo(partitionDOFs, rankDofs, rankAllDofs, boundaryDofs, vertexDof);
1159

1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
    // === BEGIN EXP

    DofIndexToBool rankAllDofIndices;
    for (DofContainer::iterator dofIt = rankAllDofs.begin(); 
	 dofIt != rankAllDofs.end(); ++dofIt) {
      rankAllDofIndices[**dofIt] = true;
    }

    for (std::map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	 it != mesh->getPeriodicAssociations().end(); ++it) {
      VertexVector *tmp = it->second;
      DOFIteratorBase it(const_cast<DOFAdmin*>(tmp->getAdmin()), USED_DOFS);
      for (it.reset(); !it.end(); ++it) {
	if (!it.isDOFFree()) {
	  if (!rankAllDofIndices[(*tmp)[it.getDOFIndex()]] &&
	      rankAllDofIndices[it.getDOFIndex()])
	    //	    (*tmp)[it.getDOFIndex()] = -1;
	    (*tmp)[it.getDOFIndex()] = it.getDOFIndex();
	}
      }
    }
    
    // === ENDE EXP

1184
    nRankDofs = rankDofs.size();
1185
1186
    nOverallDOFs = partitionDOFs.size();

1187

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

1190
    rstart = 0;
1191
1192
    mpiComm.Scan(&nRankDofs, &rstart, 1, MPI_INT, MPI_SUM);
    rstart -= nRankDofs;
1193

1194
1195
1196
1197
    
    // === 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