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

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

#include "petscksp.h"
21
22
23

namespace AMDiS {

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

    return 0;
  }

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

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

313
    int nnz = 0;
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
    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]++;
	      }
	  }
	}

344
345
346
    for (int i = 0; i < nRankRows; i++)
      if (d_nnz[i] > nnz)
	nnz = d_nnz[i];
347
348

    MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows,
349
		    nnz, NULL, nnz / 10, NULL, &petscMatrix);
350

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

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


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

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

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

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

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

    VecRestoreArray(petscSolVec, &vecPointer);

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
427

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

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

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

    MatDestroy(petscMatrix);
441
442
  }

443

444
445
446
447
448
449
450
451
452
453
454
  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);
455
    PCSetType(pc, PCSOR);
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
    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);

473

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

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

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

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

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

    MatDestroy(petscMatrix);
529
530
531
  }


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

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

566

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

581

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

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

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

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

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

	    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;

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

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

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

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

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

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


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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
715
716
717

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

Thomas Witkowski's avatar
Thomas Witkowski committed
718
719

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

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

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

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


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


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

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

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

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

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

791

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
803

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

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

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

827
828
829

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
882

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

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

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

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

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

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

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

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

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

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


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

    
936
    // === Delete send buffers. ===
937

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


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

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

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

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

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

966
	bool found = false;
967

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

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

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

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

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

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

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

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

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


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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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


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

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

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

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

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


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

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

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

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

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

    
1100

Thomas Witkowski's avatar
Thomas Witkowski committed
1101
1102
1103
    for (RankToBoundMap::