ParallelDomainBase.cc 54.3 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
248
249
250
251
252
253
254
255
256
257
    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

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

    typedef traits::range_generator<major, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

258
259
260
261
    for (cursor_type cursor = begin<major>(mat->getBaseMatrix()), 
	   cend = end<major>(mat->getBaseMatrix()); cursor != cend; ++cursor)
      for (icursor_type icursor = begin<nz>(cursor), 
	     icend = end<nz>(cursor); icursor != icend; ++icursor)
262
	if (value(*icursor) != 0.0) {
263
264
 	  int r = mapLocalGlobalDOFs[row(*icursor)] * dispMult + dispAddRow;
 	  int c = mapLocalGlobalDOFs[col(*icursor)] * dispMult + dispAddCol;  
265
	  double v = value(*icursor);
266

267
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
268
	}
269
  }
270

271

272
273
274
  void ParallelDomainBase::setDofVector(DOFVector<double>* vec, 
					int dispMult, int dispAdd)
  {
275
276
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
277
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()] * dispMult + dispAdd;
278
279
280
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
281
    }    
282
283
  }

284

285
286
  void ParallelDomainBase::fillPetscMatrix(DOFMatrix *mat, DOFVector<double> *vec)
  {
287
288
289
290
    MatCreate(PETSC_COMM_WORLD, &petscMatrix);
    MatSetSizes(petscMatrix, nRankRows, nRankRows, nOverallRows, nOverallRows);
    MatSetType(petscMatrix, MATAIJ);

291
292
293
294
295
    setDofMatrix(mat);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
298
299
    setDofVector(vec);

    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);
300
301
302
303
304
  }

  
  void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
  {
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
    FUNCNAME("ParallelDomainBase::fillPetscMatrix()");

    clock_t first = clock();

    //    MatCreate(PETSC_COMM_WORLD, &petscMatrix);
    //    MatSetSizes(petscMatrix, nRankRows, nRankRows, nOverallRows, nOverallRows);
    //MatSetType(petscMatrix, MATAIJ);

    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];
    for (int i = 0; i < nRankRows; i++)
      d_nnz[i] = 0;

    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;

		if (r < nRankRows)
		  d_nnz[r]++;
	      }
	  }
	}

    INFO(info, 8)("t1 petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));

    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
		    0, d_nnz, 0, NULL, &petscMatrix);

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

    INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock()));
368
369
370
  }


371
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
372
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

375
376
377
378
379
380
    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
    KSPGetPC(ksp, &pc);
381
382
    //PCSetType(pc, PCJACOBI);
    PCSetType(pc, PCILU);
383
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
384
385
    KSPSetType(ksp, KSPBCGS);
    //KSPSetType(ksp, KSPCG);
386
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
387
388
389
390
391
    KSPSolve(ksp, petscRhsVec, petscSolVec);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
392
    for (int i = 0; i < nRankDOFs; i++)
393
      (*vec)[mapLocalToDofIndex[i]] = vecPointer[i];
394
395
396

    VecRestoreArray(petscSolVec, &vecPointer);

397
398
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
399
400
401

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
402
403
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
404
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
405
	 sendIt != sendDofs.end(); ++sendIt, i++) {
406
407
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
408

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

Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
414
415
416
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
417
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
418
	 recvIt != recvDofs.end(); ++recvIt, i++) {
419
420
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
421

Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
424
425
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
426
427
428

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

429
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
430
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
431
	 recvIt != recvDofs.end(); ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
432
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
433
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
434
435
436
437

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

    MatDestroy(petscMatrix);
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
  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);
    //    PCSetType(pc, PCNONE);
    PCSetType(pc, PCJACOBI);
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
    //KSPSetType(ksp, KSPCG);
    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);

474

475
476
477
478
479
480
481
482
483
    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++) {
484
485
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs * nComponents];
486
487
488
489
490
491
492
493

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

494
495
      request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs * nComponents,
						MPI_DOUBLE, sendIt->first, 0);
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
    }

    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);
519
520
 	for (int k = 0; k < nRecvDOFs; k++)
 	  (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++];
521
522
523
524
525
526
527
      }

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

    MatDestroy(petscMatrix);
530
531
532
  }


533
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
  {
    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
554
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
555
556
557
558
559
560
561
562
563
564
565
566
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

567

568
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
569
570
571
572
573
574
575
576
577
578
579
580
581
  {
    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);
  }

582

583
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs)
584
  {
585
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
586
587
588

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

589
    TraverseStack stack;
590
591
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
592
593
594
595
596
597
598
599
600
601
602
603
    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));
604

605
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
606
607
	    // We have found an element that is at an interior boundary. 

608
609
	    // === 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. ===
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634

	    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;

635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
#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

653
	    // === And add the part of the interior boundary. ===
654
655

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
656
657
658
659
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

660
	    bound.rankObject.el = element;
661
	    bound.rankObject.elIndex = element->getIndex();
662
663
	    bound.rankObject.subObjAtBoundary = EDGE;
	    bound.rankObject.ithObjAtBoundary = i;
664
665
666
667
	    // 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();
668
	    bound.neighbourObject.subObjAtBoundary = EDGE;
669
670
671
672
673
	    bound.neighbourObject.ithObjAtBoundary = elInfo->getSideOfNeighbour(i);

	    // i == 2  =>  getSideOfNeighbour(i) == 2
	    TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
	      ("Should not happen!\n");
674
675
676
677
678
679
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
680
681
682


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
688
689
690
691
    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
694
      int nSendInt = rankIt->second.size();
695
696
697
698
699
      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
700
701
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
702
      request[requestCounter++] =
703
	mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
704
705
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
706
707
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
708
      int nRecvInt = rankIt->second.size();
709
      int *buffer = new int[nRecvInt * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
710
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
711

Thomas Witkowski's avatar
Thomas Witkowski committed
712
      request[requestCounter++] = 
713
	mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
718

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

Thomas Witkowski's avatar
Thomas Witkowski committed
719
720

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
721
722
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
723
724
725
726
727
728
729

      // === 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.
730
731
	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
732
733
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
734
735
	    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
736
737
738
	      break;

	  // The element must always be found, because the list is just in another order.
Thomas Witkowski's avatar
Thomas Witkowski committed
739
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
740
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
741
742
743
744
745
746
747
748
749
750
751
752
753

	  // 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];
754
755
756
  }


757
  void ParallelDomainBase::removeMacroElements()
758
759
760
761
762
763
764
765
766
767
768
769
  {
    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
770
    mesh->removeMacroElements(macrosToRemove, feSpace);
771
772
773
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
774
775
776
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
						      int& nRankDOFs, 
						      int& nOverallDOFs)
777
  {
778
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
779
780

    // === Get rank information about DOFs. ===
781
782
783

    // 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
784
    DofContainer rankAllDofs;
785
    DofToRank boundaryDofs;
786

787
    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDofs);
788
789
790
791

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

792

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

795
    rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
796
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
797
798
    rstart -= nRankDOFs;

799
800
801
802
    
    // === 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
803

Thomas Witkowski's avatar
Thomas Witkowski committed
804

805
806
807
    // 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
808
809
810
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
811
812
813
814
815
816
817
818
819
820

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

821
822
823
      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
824
      isRankDof[i] = true;
825
      i++;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
826
827
    }

828
829
830

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

831
832
    // Stores for dofs in rank a new global index.
    DofIndexMap rankDofsNewGlobalIndex;
833
834
835
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
836
837
838
    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
839
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
840
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
841
842
843
      i++;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
844
 
845
846
847
848
849
850
    // === 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;
851

852
853
854
855
    // 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;
856

857
    for (DofToRank::iterator it = boundaryDofs.begin(); it != boundaryDofs.end(); ++it) {
858
859
860
861
862
863
864
865

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

870
	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
871
	  }
872
873
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
874
875
	// 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.
876
877
878
879
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
880
881
882
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
883

884
    // === Send and receive the dof indices at boundary. ===
885

Thomas Witkowski's avatar
Thomas Witkowski committed
886
    std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size());
887
888
889

    sendDofs.clear();
    recvDofs.clear();
890
    
Thomas Witkowski's avatar
Thomas Witkowski committed
891
892
893
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
894
    i = 0;
895
896
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
897
898
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
899
900
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
901

902
      int c = 0;
903
904
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
905
	   dofIt != sendIt->second.end();
906
	   ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
907
	sendBuffers[i][c++] = *(dofIt->first);
908
	sendBuffers[i][c++] = dofIt->second;
909

910
911
912
913
914
#if 0
	if (mpiRank == 3 && sendIt->first == 2)
	  std::cout << "SEND DOF: " << dofIt->first << std::endl;
#endif

915
	sendDofs[sendIt->first].push_back(dofIt->first);
916
917
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
918
919
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
920
921
922
    }

    i = 0;
923
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
924
925
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
926
927
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
928

Thomas Witkowski's avatar
Thomas Witkowski committed
929
930
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
931
932
933
    }


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

    
937
    // === Delete send buffers. ===
938

Thomas Witkowski's avatar
Thomas Witkowski committed
939
940
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
941
942


943
    // === Change dof indices at boundary from other ranks. ===
944

Thomas Witkowski's avatar
Thomas Witkowski committed
945
946
    // 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
947
    // 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
948
949
950
951
952
    // 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;
953
    for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
954
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
955
956
      dofChanged[dofIt->first] = false;

957
    i = 0;
958
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
959
960
961
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

962
      for (int j = 0; j < recvIt->second; j++) {
963

964
965
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
966

967
	bool found = false;
968

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

971
972
	for (DofToRank::iterator dofIt = boundaryDofs.begin(); 
	     dofIt != boundaryDofs.end(); ++dofIt) {
973

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

977
978
979
980
981
#if 0
	    if (mpiRank == 2 && recvIt->first == 3)
	      std::cout << "RECV DOF: " << dofIt->first << std::endl;
#endif
	
982
	    recvDofs[recvIt->first].push_back(dofIt->first);
983
	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
984
	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
985

986
	    found = true;
987
988
989
	    break;
	  }
	}
990

Thomas Witkowski's avatar
Thomas Witkowski committed
991
	TEST_EXIT_DBG(found)("Should not happen!\n");
992
993
994
995
      }

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

997
    // === Create now the local to global index and local to dof index mappings.  ===
998

999
1000
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
1001
1002
1003
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
1004
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs)
1005
  {
1006
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
1007

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

1010
    // === Get all DOFs in ranks partition. ===
1011

Thomas Witkowski's avatar
Thomas Witkowski committed
1012
    ElementDofIterator elDofIt(feSpace);
1013
    DofSet rankDOFSet;
1014
1015
1016
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1017
      Element *element = elInfo->getElement();     
Thomas Witkowski's avatar
Thomas Witkowski committed
1018
1019
      elDofIt.reset(element);
      do {
1020
	rankDOFSet.insert(elDofIt.getDofPtr());
Thomas Witkowski's avatar
Thomas Witkowski committed
1021
      } while(elDofIt.next());
1022
1023
1024
1025

      elInfo = stack.traverseNext(elInfo);
    }

1026
    DofContainer rankAllDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
1027
1028
    for (DofSet::iterator dofIt = rankDOFSet.begin(); dofIt != rankDOFSet.end(); ++dofIt)
      rankAllDofs.push_back(*dofIt);    
1029
1030
1031
1032
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
    DofContainer rankDOFs = rankAllDofs;


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

1036
1037
    sendDofs.clear();
    recvDofs.clear();
1038

Thomas Witkowski's avatar
Thomas Witkowski committed
1039
1040
1041
    for (RankToBoundMap::iterator it =  myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end(); ++it) {

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

1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
#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


1055
1056
	DofContainer dofs;
	DofContainer &dofsToSend = sendDofs[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
1057
1058
1059

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
1060
1061
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1062
1063
	  break;
	case 1:
1064
1065
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1066
1067
	  break;
	case 2:
1068
1069
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
Thomas Witkowski's avatar
Thomas Witkowski committed
1070
1071
1072
1073
1074
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

1075
1076
1077
1078
1079
1080
1081
#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();
1082
	}
1083
1084
1085
1086
1087
1088
#endif
	

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

1090
1091
1092
1093
1094
	dofs.clear();
	addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
			 dofs);	
  	addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
  		       dofs);
1095
	for (int i = 0; i < static_cast<int>(dofs.size()); i++)
1096
	  dofsToSend.push_back(dofs[i]);
1097
      }
1098
1099
1100
    }

    
1101

Thomas Witkowski's avatar
Thomas Witkowski committed
1102
1103
1104
    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end(); ++it) {

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

1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
#if 0
	if (mpiRank == 2 && it->first == 3) 
	  std::cout << "GO ON OT BOUND: " << boundIt->rankObject.elIndex 
		    << "/" << boundIt->rankObject.ithObjAtBoundary << " with "
		    << boundIt->neighbourObject.elIndex << "/"
		    << boundIt->neighbourObject.ithObjAtBoundary
		    << std::endl;
#endif