ParallelDomainBase.cc 53.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)
27
      std::cout << "  Iteration " << iter << ": " << rnorm << std::endl;
28
29
30
31

    return 0;
  }

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

37
38
39
40
41
  ParallelDomainBase::ParallelDomainBase(const std::string& name,
					 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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
#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);
      }
    }


    if (mpiRank == 1) {
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      while (elInfo) {
	if (elInfo->getElement()->getIndex() == 3) {
	  WorldVector<double> x;
	  mesh->getDofIndexCoords(elInfo->getElement()->getDOF(2), feSpace, x);
	  std::cout << "FOUND!" << std::endl;
	  x.print();
	}

	elInfo = stack.traverseNext(elInfo);
      }
    }

    exit(0);
#endif

168

169
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
170

Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
173
    int globalRefinement = 0;
    GET_PARAMETER(0, "testgr", "%d", &globalRefinement);

Thomas Witkowski's avatar
Thomas Witkowski committed
174
175
    if (globalRefinement > 0) {
      refinementManager->globalRefine(mesh, globalRefinement);
176

177
178
179
180
181
#if (DEBUG != 0)
      elMap.clear();
      DbgCreateElementMap(elMap);
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
182
      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
183
184
185
186
187
188

      updateDofAdmins();

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

191

Thomas Witkowski's avatar
Thomas Witkowski committed
192
#if (DEBUG != 0)
193
    DbgTestCommonDofs(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
194
#endif
195

196

197
    // === Create petsc matrix. ===
198

199
200
    int nRankRows = nRankDOFs * nComponents;
    int nOverallRows = nOverallDOFs * nComponents;
201

202
203
204
    MatCreate(PETSC_COMM_WORLD, &petscMatrix);
    MatSetSizes(petscMatrix, nRankRows, nRankRows, nOverallRows, nOverallRows);
    MatSetType(petscMatrix, MATAIJ);
205

206
207
208
209
210
211
212
    VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
    VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
    VecSetType(petscRhsVec, VECMPI);

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

215
216
217

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
  {
218
219
220
    MatDestroy(petscMatrix);
    VecDestroy(petscRhsVec);
    VecDestroy(petscSolVec);
221
  }
222

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

241

242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
  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");
  }

263
264
  void ParallelDomainBase::setDofMatrix(DOFMatrix* mat, int dispMult, 
					int dispAddRow, int dispAddCol)
265
  {
266
267
268
269
    FUNCNAME("ParallelDomainBase::setDofMatrix()");

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

270
271
272
273
274
275
276
277
278
279
280
    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;

281
282
283
284
    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)
285
	if (value(*icursor) != 0.0) {
286
287
 	  int r = mapLocalGlobalDOFs[row(*icursor)] * dispMult + dispAddRow;
 	  int c = mapLocalGlobalDOFs[col(*icursor)] * dispMult + dispAddCol;  
288
	  double v = value(*icursor);
289

290
	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
291
	}
292
  }
293

294

295
296
297
  void ParallelDomainBase::setDofVector(DOFVector<double>* vec, 
					int dispMult, int dispAdd)
  {
298
299
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
300
      int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()] * dispMult + dispAdd;
301
302
303
      double value = *dofIt;

      VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
304
    }    
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
  void ParallelDomainBase::fillPetscMatrix(DOFMatrix *mat, DOFVector<double> *vec)
  {
    setDofMatrix(mat);

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

    setDofVector(vec);   
  }

  
  void ParallelDomainBase::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
  {
    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);
  }


338
  void ParallelDomainBase::solvePetscMatrix(DOFVector<double> *vec)
339
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
340
341
    FUNCNAME("ParallelDomainBase::solvePetscMatrix()");

342
343
344
345
346
347
    KSP ksp;
    PC pc;

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN);
    KSPGetPC(ksp, &pc);
348
    //    PCSetType(pc, PCNONE);
349
    PCSetType(pc, PCILU);
350
    KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
351
352
    KSPSetType(ksp, KSPBCGS);
    //KSPSetType(ksp, KSPCG);
353
    KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0);
354
355
    KSPSolve(ksp, petscRhsVec, petscSolVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
358
359
360
361
#if (DEBUG != 0)
    int size = 0;
    VecGetLocalSize(petscSolVec, &size);
    TEST_EXIT(size == nRankDOFs)("Vector and rank DOFs does not fit together!\n");
#endif

362
363
364
    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

Thomas Witkowski's avatar
Thomas Witkowski committed
365
    for (int i = 0; i < nRankDOFs; i++)
366
      (*vec)[mapLocalToDofIndex[i]] = vecPointer[i];
367
368
369

    VecRestoreArray(petscSolVec, &vecPointer);

370
371
    std::vector<double*> sendBuffers(sendDofs.size());
    std::vector<double*> recvBuffers(recvDofs.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
374

    MPI::Request request[sendDofs.size() + recvDofs.size()];
    int requestCounter = 0;
375
376
    
    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
377
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
378
	 sendIt != sendDofs.end(); ++sendIt, i++) {
379
380
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs];
381

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

Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0);
387
388
389
    }

    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
390
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
391
	 recvIt != recvDofs.end(); ++recvIt, i++) {
392
393
      int nRecvDOFs = recvIt->second.size();
      recvBuffers[i] = new double[nRecvDOFs];
394

Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0);
397
398
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
399
400
401

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

402
    i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
403
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
404
	 recvIt != recvDofs.end(); ++recvIt, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
405
      for (int j = 0; j < static_cast<int>(recvIt->second.size()); j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
406
	(*vec)[*(recvIt->second)[j]] = recvBuffers[i][j];
407
408
409
410

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

415

416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
  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);

#if (DEBUG != 0)
    int size = 0;
    VecGetLocalSize(petscSolVec, &size);
437
438
    TEST_EXIT(size == nRankDOFs * nComponents)
      ("Vector and rank DOFs does not fit together!\n");
439
440
441
442
443
444
445
446
447
448
449
450
451
#endif

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

452

453
454
455
456
457
458
459
460
461
    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++) {
462
463
      int nSendDOFs = sendIt->second.size();
      sendBuffers[i] = new double[nSendDOFs * nComponents];
464
465
466
467
468
469
470
471

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

472
473
      request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs * nComponents,
						MPI_DOUBLE, sendIt->first, 0);
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
    }

    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);
497
498
 	for (int k = 0; k < nRecvDOFs; k++)
 	  (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++];
499
500
501
502
503
504
505
506
507
508
      }

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


509
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
  {
    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
530
	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
531
532
533
534
535
536
537
538
539
540
541
542
	if (element->isLeaf()) {
	  elemWeights[elNum] += 1.0;
	  localWeightSum += 1.0;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return localWeightSum;
  }

543

544
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
545
546
547
548
549
550
551
552
553
554
555
556
557
  {
    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);
  }

558

559
  void ParallelDomainBase::createInteriorBoundaryInfo(DofContainer& rankDOFs)
560
  {
561
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
562
563
564

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

565
    TraverseStack stack;
566
567
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
568
569
570
571
572
573
574
575
576
577
578
579
    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));
580

581
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
582
583
	    // We have found an element that is at an interior boundary. 

584
585
	    // === 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. ===
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610

	    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;

611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
#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

629
	    // === And add the part of the interior boundary. ===
630
631

	    AtomicBoundary& bound = 
Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
634
635
	      (ranksBoundary ?
	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));

636
	    bound.rankObject.el = element;
637
	    bound.rankObject.elIndex = element->getIndex();
638
639
	    bound.rankObject.subObjAtBoundary = EDGE;
	    bound.rankObject.ithObjAtBoundary = i;
640
641
642
643
	    // 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();
644
	    bound.neighbourObject.subObjAtBoundary = EDGE;
645
646
647
648
649
	    bound.neighbourObject.ithObjAtBoundary = elInfo->getSideOfNeighbour(i);

	    // i == 2  =>  getSideOfNeighbour(i) == 2
	    TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
	      ("Should not happen!\n");
650
651
652
653
654
655
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
656
657
658


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
666
667
    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
668
669
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
670
      int nSendInt = rankIt->second.size();
671
672
673
674
675
      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
676
677
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
678
      request[requestCounter++] =
679
	mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
680
681
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
684
      int nRecvInt = rankIt->second.size();
685
      int *buffer = new int[nRecvInt * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
686
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
687

Thomas Witkowski's avatar
Thomas Witkowski committed
688
      request[requestCounter++] = 
689
	mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
690
691
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
694

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

Thomas Witkowski's avatar
Thomas Witkowski committed
695
696

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
697
698
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
701
702
703
704
705

      // === 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.
706
707
	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
708
709
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
710
711
	    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
712
713
714
	      break;

	  // The element must always be found, because the list is just in another order.
Thomas Witkowski's avatar
Thomas Witkowski committed
715
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
716
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
717
718
719
720
721
722
723
724
725
726
727
728
729

	  // 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];
730
731
732
  }


733
  void ParallelDomainBase::removeMacroElements()
734
735
736
737
738
739
740
741
742
743
744
745
  {
    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
746
    mesh->removeMacroElements(macrosToRemove, feSpace);
747
748
749
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
750
751
752
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDOFs,
						      int& nRankDOFs, 
						      int& nOverallDOFs)
753
  {
754
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
755
756

    // === Get rank information about DOFs. ===
757
758
759

    // 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
760
    DofContainer rankAllDofs;
761
    DofToRank boundaryDofs;
762

763
    createDOFMemberInfo(partitionDOFs, rankDOFs, rankAllDofs, boundaryDofs);
764
765
766
767

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

768

769
    // === Get starting position for global rank dof ordering. ====
770

771
    rstart = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
772
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
773
774
    rstart -= nRankDOFs;

775
776
777
778
    
    // === 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
779

Thomas Witkowski's avatar
Thomas Witkowski committed
780

781
782
783
    // 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
784
785
786
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
787
788
789
      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
790
      isRankDof[i] = true;
791
      i++;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
792
793
    }

794
795
796

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

797
798
    // Stores for dofs in rank a new global index.
    DofIndexMap rankDofsNewGlobalIndex;
799
800
801
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
802
803
804
    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
805
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
806
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
807
808
809
      i++;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
810
 
811
812
813
814
815
816
    // === 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;
817

818
819
820
821
    // 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;
822

823
    for (DofToRank::iterator it = boundaryDofs.begin(); it != boundaryDofs.end(); ++it) {
824
825
826
827
828
829
830
831

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

836
	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
837
	  }
838
839
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
840
841
	// 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.
842
843
844
845
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
846
847
848
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
849

850
    // === Send and receive the dof indices at boundary. ===
851

Thomas Witkowski's avatar
Thomas Witkowski committed
852
    std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size());
853
854
855

    sendDofs.clear();
    recvDofs.clear();
856
    
Thomas Witkowski's avatar
Thomas Witkowski committed
857
858
859
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
860
    i = 0;
861
862
    for (std::map<int, std::map<const DegreeOfFreedom*, DegreeOfFreedom> >::iterator 
	   sendIt = sendNewDofs.begin();
863
864
	 sendIt != sendNewDofs.end();
	 ++sendIt, i++) {
865
866
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
867

868
      int c = 0;
869
870
      for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
	     dofIt = sendIt->second.begin();
871
	   dofIt != sendIt->second.end();
872
	   ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
873
	sendBuffers[i][c++] = *(dofIt->first);
874
	sendBuffers[i][c++] = dofIt->second;
875

876
877
878
879
880
#if 0
	if (mpiRank == 3 && sendIt->first == 2)
	  std::cout << "SEND DOF: " << dofIt->first << std::endl;
#endif

881
	sendDofs[sendIt->first].push_back(dofIt->first);
882
883
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
884
885
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
886
887
888
    }

    i = 0;
889
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
890
891
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {
892
893
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
894

Thomas Witkowski's avatar
Thomas Witkowski committed
895
896
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
897
898
899
    }


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

    
903
    // === Delete send buffers. ===
904

Thomas Witkowski's avatar
Thomas Witkowski committed
905
906
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
907
908


909
    // === Change dof indices at boundary from other ranks. ===
910

Thomas Witkowski's avatar
Thomas Witkowski committed
911
912
    // 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
913
    // 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
914
915
916
917
918
    // 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;
919
    for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
920
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
921
922
      dofChanged[dofIt->first] = false;

923
    i = 0;
924
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
925
926
927
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

928
      for (int j = 0; j < recvIt->second; j++) {
929

930
931
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
932

933
	bool found = false;
934

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

937
938
	for (DofToRank::iterator dofIt = boundaryDofs.begin(); 
	     dofIt != boundaryDofs.end(); ++dofIt) {
939

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

943
944
945
946
947
#if 0
	    if (mpiRank == 2 && recvIt->first == 3)
	      std::cout << "RECV DOF: " << dofIt->first << std::endl;
#endif
	
948
	    recvDofs[recvIt->first].push_back(dofIt->first);
949
	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
950
	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
951

952
	    found = true;
953
954
955
	    break;
	  }
	}
956

Thomas Witkowski's avatar
Thomas Witkowski committed
957
	TEST_EXIT_DBG(found)("Should not happen!\n");
958
959
960
961
      }

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

963
    // === Create now the local to global index and local to dof index mappings.  ===
964

965
966
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
967
968
969
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
970
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs)
971
  {
972
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
973

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

976
    // === Get all DOFs in ranks partition. ===
977

Thomas Witkowski's avatar
Thomas Witkowski committed
978
    ElementDofIterator elDofIt(feSpace);
979
    DofSet rankDOFSet;
980
981
982
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
983
      Element *element = elInfo->getElement();     
Thomas Witkowski's avatar
Thomas Witkowski committed
984
985
      elDofIt.reset(element);
      do {
986
	rankDOFSet.insert(elDofIt.getDofPtr());
Thomas Witkowski's avatar
Thomas Witkowski committed
987
      } while(elDofIt.next());
988
989
990
991

      elInfo = stack.traverseNext(elInfo);
    }

992
    DofContainer rankAllDofs;
993
994
995
996
997
998
999
    for (DofSet::iterator dofIt = rankDOFSet.begin(); dofIt != rankDOFSet.end(); ++dofIt) {
      /*      if (mpiRank == 1) {
	std::cout << "COORDs of dof = " << **dofIt << std::endl;
	WorldVector<double> x;
	mesh->getDofIndexCoords(*dofIt, feSpace, x);
	x.print();	
	}*/
1000
      rankAllDofs.push_back(*dofIt);
1001
    }
1002
1003
1004
1005
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
    DofContainer rankDOFs = rankAllDofs;


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

1009
1010
    sendDofs.clear();
    recvDofs.clear();
1011

Thomas Witkowski's avatar
Thomas Witkowski committed
1012
1013
1014
    for (RankToBoundMap::iterator it =  myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end(); ++it) {

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

1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
#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


1028
1029
	DofContainer dofs;
	DofContainer &dofsToSend = sendDofs[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
1030
1031
1032

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
1033
1034
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1035
1036
	  break;
	case 1:
1037
1038
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1039
1040
	  break;
	case 2:
1041
1042
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
Thomas Witkowski's avatar
Thomas Witkowski committed
1043
1044
1045
1046
1047
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

1048
1049
1050
1051
1052
1053
1054
#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();
1055
	}
1056
1057
1058
1059
1060
1061
#endif
	

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

1063
1064
1065
1066
1067
	dofs.clear();
	addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
			 dofs);	
  	addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
  		       dofs);
1068
	for (int i = 0; i < static_cast<int>(dofs.size()); i++)
1069
	  dofsToSend.push_back(dofs[i]);
1070
      }
1071
1072
1073
    }

    
1074

Thomas Witkowski's avatar
Thomas Witkowski committed
1075
1076
1077
    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end(); ++it) {

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

1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
#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


1091
1092
	DofContainer dofs;
	DofContainer &dofsToRecv = recvDofs[it->first];
1093
1094
1095

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
1096
1097
1098
1099
1100
1101
1102
	  if (boundIt->neighbourObject.ithObjAtBoundary == 0) {
	    dofs.push_back(boundIt->rankObject.el->getDOF(2));
	    dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  } else {
	    dofs.push_back(boundIt->rankObject.el->getDOF(1));
	    dofs.push_back(boundIt->rankObject.el->getDOF(2));
	  }
1103
1104
	  break;
	case 1:
1105
1106
1107
1108
1109
1110
1111
	  if (boundIt->neighbourObject.ithObjAtBoundary == 0) {
	    dofs.push_back(boundIt->rankObject.el->getDOF(0));
	    dofs.push_back(boundIt->rankObject.el->getDOF(2));
	  } else {
	    dofs.push_back(boundIt->rankObject.el->getDOF(2));
	    dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  }
1112
1113
	  break;
	case 2: