ParallelDomainBase.cc 54.6 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);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
37
  ParallelDomainBase::ParallelDomainBase(std::string name,
38
39
40
41
					 ProblemIterationInterface *iIF,
					 ProblemTimeInterface *tIF,
					 FiniteElemSpace *fe,
					 RefinementManager *refineManager)
42
43
    : iterationIF(iIF),
      timeIF(tIF),
44
45
      feSpace(fe),
      mesh(fe->getMesh()),
46
      refinementManager(refineManager),
47
      initialPartitionMesh(true),
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
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
75
76
77
78
79
80
81
82
83
84
85
86
#if 0
    if (mpiRank == 0) {
      std::map<int, double> vec;
      
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					   Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
      while (elInfo) {
	vec[elInfo->getElement()->getIndex()] = static_cast<double>(elInfo->getElement()->getIndex());
	elInfo = stack.traverseNext(elInfo);
      }

      ElementFileWriter::writeFile(vec, feSpace, "test.vtu");
    }
#endif


87
88
89
90
91
    // 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();

92
93
94
95
96
97
98
    // 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
99
100
101
102
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
    DbgCreateElementMap(elMap);
#endif
103

104
    // === Create new global and local DOF numbering. ===
105

106
107
108
    // Set of all DOFs of the rank.
    std::vector<const DegreeOfFreedom*> rankDOFs;
    // Number of DOFs in ranks partition that are owned by the rank.
Thomas Witkowski's avatar
Thomas Witkowski committed
109
    nRankDOFs = 0;
110
    // Number of all DOFs in the macro mesh.
111
    int nOverallDOFs = 0;
112

113
    createLocalGlobalNumbering(rankDOFs, nRankDOFs, nOverallDOFs);
114

Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
    // === Create interior boundary information ===

117
    createInteriorBoundaryInfo(rankDOFs);
Thomas Witkowski's avatar
Thomas Witkowski committed
118

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
119
120
121
122
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

Thomas Witkowski's avatar
Thomas Witkowski committed
123
#if (DEBUG != 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
124
    DbgTestElementMap(elMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
125
126
    DbgTestInteriorBoundary();
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
127

128

129
    // === Reset all DOFAdmins of the mesh. ===
130

131
    updateDofAdmins();
132

133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#if 0
    if (mpiRank == 0) {
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      while (elInfo) {
	if (elInfo->getElement()->getIndex() == 4) {
	  WorldVector<double> x;
	  mesh->getDofIndexCoords(elInfo->getElement()->getDOF(0), feSpace, x);
	  std::cout << "FOUND!" << std::endl;
	  x.print();
	}

	elInfo = stack.traverseNext(elInfo);
      }
    }
#endif

150

151
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
152

Thomas Witkowski's avatar
Thomas Witkowski committed
153
154
155
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

Thomas Witkowski's avatar
Thomas Witkowski committed
156
157
    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);
158

159
160
161
162
163
#if (DEBUG != 0)
      elMap.clear();
      DbgCreateElementMap(elMap);
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
164
      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
165
166
167
168
169
170

      updateDofAdmins();

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

173

Thomas Witkowski's avatar
Thomas Witkowski committed
174
#if (DEBUG != 0)
175
    DbgTestCommonDofs(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
176
#endif
177

178

179
    // === Create petsc matrix. ===
180

181
182
183
    nRankRows = nRankDOFs * nComponents;
    nOverallRows = nOverallDOFs * nComponents;
    
184
185
186
187
188
189
190
    VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
    VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
    VecSetType(petscRhsVec, VECMPI);

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

193
194
195

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
  {
196
197
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
198
  }
199

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

218

219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
  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");
  }

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

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

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

254
    typedef traits::range_generator<row, Matrix>::type cursor_type;
255
256
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

257
258
259
260
261
262
263
264
265
266
267
268
269
    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;

270
271
      for (icursor_type icursor = begin<nz>(cursor), 
	     icend = end<nz>(cursor); icursor != icend; ++icursor)
272
	if (value(*icursor) != 0.0) {
273
274
	  cols.push_back(mapLocalGlobalDOFs[col(*icursor)] * dispMult + dispAddCol);
	  values.push_back(value(*icursor));
275
	}
276
277
278
279

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

280
  }
281

282

283
284
285
  void ParallelDomainBase::setDofVector(DOFVector<double>* vec, 
					int dispMult, int dispAdd)
  {
286
287
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
288
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()] * dispMult + dispAdd;
289
290
291
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
292
    }    
293
294
  }

295

296
297
  void ParallelDomainBase::fillPetscMatrix(DOFMatrix *mat, DOFVector<double> *vec)
  {
298
299
300
301
    MatCreate(PETSC_COMM_WORLD, &petscMatrix);
    MatSetSizes(petscMatrix, nRankRows, nRankRows, nOverallRows, nOverallRows);
    MatSetType(petscMatrix, MATAIJ);

302
303
304
305
306
    setDofMatrix(mat);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
307
308
309
310
    setDofVector(vec);

    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);
311
312
313
314
315
  }

  
  void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
  {
316
317
318
319
320
321
322
323
324
    FUNCNAME("ParallelDomainBase::fillPetscMatrix()");

    clock_t first = clock();

    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    int d_nnz[nRankRows];
325
326
327
    int o_nnz[nRankRows];

    for (int i = 0; i < nRankRows; i++) {
328
      d_nnz[i] = 0;
329
330
      o_nnz[i] = 0;
    }
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351

    for (int i = 0; i < nComponents; i++) 
      for (int j = 0; j < nComponents; j++) 
 	if ((*mat)[i][j]) {
	  Matrix bmat = (*mat)[i][j]->getBaseMatrix();

	  traits::row<Matrix>::type row(bmat);
	  traits::col<Matrix>::type col(bmat);
	  traits::const_value<Matrix>::type value(bmat);
	  
	  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>(bmat), 
		 cend = end<major>(bmat); cursor != cend; ++cursor) {
	    for (icursor_type icursor = begin<nz>(cursor), 
		   icend = end<nz>(cursor); icursor != icend; ++icursor)
	      if (value(*icursor) != 0.0) {
		int r = mapLocalGlobalDOFs[row(*icursor)] * nComponents + i;
		r -= rstart * nComponents;

352
353
354
355
356
357
358
359
360
361
		int c = mapLocalGlobalDOFs[col(*icursor)] * nComponents + j;

		if (r < nRankRows) {

		  if (c >= rstart * nComponents && 
		      c < rstart * nComponents + nRankRows)
		    d_nnz[r]++;
		  else
		    o_nnz[r]++;		  
		}
362
363
364
365
366
367
	      }
	  }
	}


    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
368
369
370
371
372
373
374
375
		    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
376

377
378
379
380
381
382
383
384
385
386
387
388
389
390
    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++)
	if ((*mat)[i][j]) 
	  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);
391
392

    INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));
393
394
395
  }


396
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
397
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

400
401
402
403
404
405
    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
    KSPGetPC(ksp, &pc);
406
407
    //PCSetType(pc, PCJACOBI);
    PCSetType(pc, PCILU);
408
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
409
410
    KSPSetType(ksp, KSPBCGS);
    //KSPSetType(ksp, KSPCG);
411
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
412
413
414
415
416
    KSPSolve(ksp, petscRhsVec, petscSolVec);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
417
    for (int i = 0; i < nRankDOFs; i++)
418
      (*vec)[mapLocalToDofIndex[i]] = vecPointer[i];
419
420
421

    VecRestoreArray(petscSolVec, &vecPointer);

422
423
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
424
425
426

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
427
428
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
429
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
430
	 sendIt != sendDofs.end(); ++sendIt, i++) {
431
432
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
433

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

Thomas Witkowski's avatar
Thomas Witkowski committed
437
438
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
439
440
441
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
442
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
443
	 recvIt != recvDofs.end(); ++recvIt, i++) {
444
445
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
446

Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
449
450
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
453

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

454
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
455
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
456
	 recvIt != recvDofs.end(); ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
457
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
458
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
459
460
461
462

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

    MatDestroy(petscMatrix);
467
468
  }

469

470
471
472
473
474
475
476
477
478
479
  void ParallelDomainBase::solvePetscMatrix(SystemVector *vec)
  {
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
    KSPGetPC(ksp, &pc);
480
481
    PCSetType(pc, PCKSP);
    KSPSetTolerances(ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
    KSPSetType(ksp, KSPBCGS);
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
    KSPSolve(ksp, petscRhsVec, petscSolVec);

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

    for (int i = 0; i < nComponents; i++) {
      DOFVector<double> *dofvec = vec->getDOFVector(i);
      for (int j = 0; j < nRankDOFs; j++)
	(*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i];
    }

    VecRestoreArray(petscSolVec, &vecPointer);

497

498
499
500
501
502
503
504
505
506
    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++) {
507
508
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs * nComponents];
509
510
511
512
513
514
515
516

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

517
518
      request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs * nComponents,
						MPI_DOUBLE, sendIt->first, 0);
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
    }

    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++) {
	DOFVector<double> *dofvec = vec->getDOFVector(j);
542
543
 	for (int k = 0; k < nRecvDOFs; k++)
 	  (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++];
544
545
546
547
548
549
550
      }

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

    MatDestroy(petscMatrix);
553
554
555
  }


556
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
  {
    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
577
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
578
579
580
581
582
583
584
585
586
587
588
589
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

590

591
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
592
593
594
595
596
597
598
599
600
601
602
603
604
  {
    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);
  }

605

606
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs)
607
  {
608
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
609
610
611

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

612
    TraverseStack stack;
613
614
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
615
616
617
618
619
620
621
622
623
624
625
626
    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));
627

628
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
629
630
	    // We have found an element that is at an interior boundary. 

631
632
	    // === 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. ===
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657

	    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;

658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
#if 0
	    if (mpiRank == 3 && ranksBoundary && 
		partitionVec[elInfo->getNeighbour(i)->getIndex()] == 2) {
	      std::cout << "ADD MY BOUND " << element->getIndex() << "/" << i 
			<< "  with  " 
			<< elInfo->getNeighbour(i)->getIndex() << "/"
			<< elInfo->getSideOfNeighbour(i) << std::endl;
	    }

	    if (mpiRank == 2 && !ranksBoundary && 
		partitionVec[elInfo->getNeighbour(i)->getIndex()] == 3) {
	      std::cout << "ADD OT BOUND " << element->getIndex() << "/" << i 
			<< "  with  " 
			<< elInfo->getNeighbour(i)->getIndex() << "/"
			<< elInfo->getSideOfNeighbour(i) << std::endl;
	    }
#endif

676
	    // === And add the part of the interior boundary. ===
677
678

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
681
682
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

683
	    bound.rankObject.el = element;
684
	    bound.rankObject.elIndex = element->getIndex();
685
686
	    bound.rankObject.subObjAtBoundary = EDGE;
	    bound.rankObject.ithObjAtBoundary = i;
687
688
689
690
	    // 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();
691
	    bound.neighbourObject.subObjAtBoundary = EDGE;
692
693
694
695
696
	    bound.neighbourObject.ithObjAtBoundary = elInfo->getSideOfNeighbour(i);

	    // i == 2  =>  getSideOfNeighbour(i) == 2
	    TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
	      ("Should not happen!\n");
697
698
699
700
701
702
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
703
704
705


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
711
712
713
714
    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
715
716
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
717
      int nSendInt = rankIt->second.size();
718
719
720
721
722
      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
723
724
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
725
      request[requestCounter++] =
726
	mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
727
728
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
729
730
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
731
      int nRecvInt = rankIt->second.size();
732
      int *buffer = new int[nRecvInt * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
733
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
734

Thomas Witkowski's avatar
Thomas Witkowski committed
735
      request[requestCounter++] = 
736
	mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
737
738
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
741

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

Thomas Witkowski's avatar
Thomas Witkowski committed
742
743

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
746
747
748
749
750
751
752

      // === 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.
753
754
	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
755
756
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
757
758
	    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
759
760
761
	      break;

	  // The element must always be found, because the list is just in another order.
Thomas Witkowski's avatar
Thomas Witkowski committed
762
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
763
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
764
765
766
767
768
769
770
771
772
773
774
775
776

	  // 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];
777
778
779
  }


780
  void ParallelDomainBase::removeMacroElements()
781
782
783
784
785
786
787
788
789
790
791
792
  {
    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
793
    mesh->removeMacroElements(macrosToRemove, feSpace);
794
795
796
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
797
798
799
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
						      int& nRankDOFs, 
						      int& nOverallDOFs)
800
  {
801
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
802
803

    // === Get rank information about DOFs. ===
804
805
806

    // 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
807
    DofContainer rankAllDofs;
808
    DofToRank boundaryDofs;
809

810
    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDofs);
811
812
813
814

    nRankDOFs = rankDOFs.size();
    nOverallDOFs = partitionDOFs.size();

815

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

818
    rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
819
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
820
821
    rstart -= nRankDOFs;

822
823
824
825
    
    // === 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
826

Thomas Witkowski's avatar
Thomas Witkowski committed
827

828
829
830
    // 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
831
832
833
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
834
835
836
837
838
839
840
841
842
843

#if 0
      if (mpiRank == 0) {
	std::cout << "COORDs of dof = " << i << std::endl;
	WorldVector<double> x;
	mesh->getDofIndexCoords(*dofIt, feSpace, x);
	x.print();	
      }
#endif

844
845
846
      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
847
      isRankDof[i] = true;
848
      i++;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
849
850
    }

851
852
853

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

854
855
    // Stores for dofs in rank a new global index.
    DofIndexMap rankDofsNewGlobalIndex;
856
857
858
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
859
860
861
    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
862
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
863
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
864
865
866
      i++;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
867
 
868
869
870
871
872
873
    // === 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;
874

875
876
877
878
    // 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;
879

880
    for (DofToRank::iterator it = boundaryDofs.begin(); it != boundaryDofs.end(); ++it) {
881
882
883
884
885
886
887
888

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

893
	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
894
	  }
895
896
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
897
898
	// 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.
899
900
901
902
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
903
904
905
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
906

907
    // === Send and receive the dof indices at boundary. ===
908

Thomas Witkowski's avatar
Thomas Witkowski committed
909
    std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size());
910
911
912

    sendDofs.clear();
    recvDofs.clear();
913
    
Thomas Witkowski's avatar
Thomas Witkowski committed
914
915
916
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
917
    i = 0;
918
919
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
920
921
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
922
923
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
924

925
      int c = 0;
926
927
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
928
	   dofIt != sendIt->second.end();
929
	   ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
930
	sendBuffers[i][c++] = *(dofIt->first);
931
	sendBuffers[i][c++] = dofIt->second;
932

933
934
935
936
937
#if 0
	if (mpiRank == 3 && sendIt->first == 2)
	  std::cout << "SEND DOF: " << dofIt->first << std::endl;
#endif

938
	sendDofs[sendIt->first].push_back(dofIt->first);
939
940
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
941
942
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
943
944
945
    }

    i = 0;
946
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
947
948
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
949
950
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
951

Thomas Witkowski's avatar
Thomas Witkowski committed
952
953
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
954
955
956
    }


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

    
960
    // === Delete send buffers. ===
961

Thomas Witkowski's avatar
Thomas Witkowski committed
962
963
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
964
965


966
    // === Change dof indices at boundary from other ranks. ===
967

Thomas Witkowski's avatar
Thomas Witkowski committed
968
969
    // 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
970
    // 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
971
972
973
974
975
    // 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;
976
    for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
977
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
978
979
      dofChanged[dofIt->first] = false;

980
    i = 0;
981
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
982
983
984
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

985
      for (int j = 0; j < recvIt->second; j++) {
986

987
988
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
989

990
	bool found = false;
991

Thomas Witkowski's avatar
Thomas Witkowski committed
992
993
	// Iterate over all boundary dofs to find the dof, which index we have to change.

994
995
	for (DofToRank::iterator dofIt = boundaryDofs.begin(); 
	     dofIt != boundaryDofs.end(); ++dofIt) {
996

Thomas Witkowski's avatar
Thomas Witkowski committed
997
998
999
	  if (*(dofIt->first) == oldDof && !dofChanged[dofIt->first]) {
	    dofChanged[dofIt->first] = true;

1000
1001
1002
1003
1004
#if 0
	    if (mpiRank == 2 && recvIt->first == 3)
	      std::cout << "RECV DOF: " << dofIt->first << std::endl;
#endif
	
1005
	    recvDofs[recvIt->first].push_back(dofIt->first);
1006
	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1007
	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
1008

1009
	    found = true;
1010
1011
1012
	    break;
	  }
	}
1013

Thomas Witkowski's avatar
Thomas Witkowski committed
1014
	TEST_EXIT_DBG(found)("Should not happen!\n");
1015
1016
1017
1018
      }

      delete [] recvBuffers[i];
    }
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1019

1020
    // === Create now the local to global index and local to dof index mappings.  ===
1021

1022
1023
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
1024
1025
1026
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
1027
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs)
1028
  {
1029
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
1030

1031
    typedef std::set<const DegreeOfFreedom*> DofSet;
Thomas Witkowski's avatar
Thomas Witkowski committed
1032

1033
    // === Get all DOFs in ranks partition. ===
1034

Thomas Witkowski's avatar
Thomas Witkowski committed
1035
    ElementDofIterator elDofIt(feSpace);
1036
    DofSet rankDOFSet;
1037
1038
1039
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1040
      Element *element = elInfo->getElement();     
Thomas Witkowski's avatar
Thomas Witkowski committed
1041
1042
      elDofIt.reset(element);
      do {
1043
	rankDOFSet.insert(elDofIt.getDofPtr());
Thomas Witkowski's avatar
Thomas Witkowski committed
1044
      } while(elDofIt.next());
1045
1046
1047
1048

      elInfo = stack.traverseNext(elInfo);
    }

1049
    DofContainer rankAllDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
1050
1051
    for (DofSet::iterator dofIt = rankDOFSet.begin(); dofIt != rankDOFSet.end(); ++dofIt)
      rankAllDofs.push_back(*dofIt);    
1052
1053
1054
1055
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
    DofContainer rankDOFs = rankAllDofs;


Thomas Witkowski's avatar
Thomas Witkowski committed
1056
    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
1057
    // === rankDOFs to boundaryDOFs.                                               ===
1058

1059
1060
    sendDofs.clear();
    recvDofs.clear();
1061

Thomas Witkowski's avatar
Thomas Witkowski committed
1062
1063
1064
    for (RankToBoundMap::iterator it =  myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end(); ++it) {

Thomas Witkowski's avatar
Thomas Witkowski committed
1065
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1066
1067
	   boundIt != it->second.end(); ++boundIt) {

1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
#if 0
	if (mpiRank == 3 && it->first == 2) 
	  std::cout << "GO ON MY BOUND: " << boundIt->rankObject.elIndex 
		    << "/" << boundIt->rankObject.ithObjAtBoundary << " with "
		    << boundIt->neighbourObject.elIndex << "/"
		    << boundIt->neighbourObject.ithObjAtBoundary
		    << std::endl;
#endif


1078
1079
	DofContainer dofs;
	DofContainer &dofsToSend = sendDofs[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
1080
1081
1082

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
1083
1084
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1085
1086
	  break;
	case 1:
1087
1088
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1089
1090
	  break;
	case 2:
1091
1092
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
Thomas Witkowski's avatar
Thomas Witkowski committed
1093
1094
1095
1096
1097
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

1098
1099
1100
1101
1102
1103
1104
#if 0
	if (mpiRank == 3 && it->first == 2) {
	  WorldVector<double> x;
	  mesh->getDofIndexCoords(dofs[0], feSpace, x);
	  x.print();
	  mesh->getDofIndexCoords(dofs[1], feSpace, x);
	  x.print();
1105
	}
1106
1107
1108
1109
1110
1111
#endif
	

	for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt)
	  if (find(dofsToSend.begin(), dofsToSend.end(), *dofIt) == dofsToSend.end())
	    dofsToSend.push_back(*dofIt);	
1112

1113
1114
1115