ParallelDomainBase.cc 67.4 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
#include <algorithm>
2
3
#include <iostream>
#include <fstream>
Thomas Witkowski's avatar
Thomas Witkowski committed
4
#include <boost/lexical_cast.hpp>
5
#include "parallel/ParallelDomainBase.h"
6
#include "parallel/ParallelDomainDbg.h"
7
#include "parallel/StdMpi.h"
8
9
10
11
12
13
14
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
15
16
#include "DOFMatrix.h"
#include "DOFVector.h"
17
#include "SystemVector.h"
18
#include "VtkWriter.h"
19
#include "ElementDofIterator.h"
20
21
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
22
#include "ElementFileWriter.h"
23
#include "VertexVector.h"
24
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
25
26
#include "ProblemVec.h"
#include "ProblemInstat.h"
27
#include "Debug.h"
28

29
30
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
31
32
  using boost::lexical_cast;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
38
39
40
41
42
43
44
45
46
47
  ParallelDomainBase::ParallelDomainBase(ProblemVec *problemStat,
					 ProblemInstatVec *problemInstat)
    : iterationIF(problemStat),
      timeIF(problemInstat),
      probStat(problemStat),
      name(problemStat->getName()),
      feSpace(problemStat->getFESpace(0)),
      mesh(feSpace->getMesh()),
      refineManager(problemStat->getRefinementManager(0)),
      info(problemStat->getInfo()),
48
      initialPartitionMesh(true),
49
      nRankDofs(0),
50
      rstart(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
51
      nComponents(problemStat->getNumComponents()),
52
53
      deserialized(false),
      lastMeshChangeIndex(0)
54
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
55
56
57
58
59
60
61
    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");

62
63
64
65
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90

    // Test if all fe spaces are equal. Yet, different fe spaces are not supported for
    // domain parallelization.
    const FiniteElemSpace *fe = probStat->getFESpace(0);
    for (int i = 0; i < nComponents; i++)
      TEST_EXIT(fe == probStat->getFESpace(i))
	("Parallelization does not supported different FE spaces!\n");

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
    GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
    if (writeSerialization)
      probStat->getFileWriterList().push_back(new Serializer<ParallelDomainBase>(this));

    int readSerialization = 0;
    GET_PARAMETER(0, name + "->input->read serialization", "%d", &readSerialization);
    if (readSerialization) {
      std::string filename = "";
      GET_PARAMETER(0, name + "->input->serialization filename", &filename);
      filename += ".p" + lexical_cast<std::string>(mpiRank);
      MSG("Start serialization with %s\n", filename.c_str());
      std::ifstream in(filename.c_str());
      deserialize(in);
      in.close();
    }
91
92
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
93

94
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
95
  {
96
97
98
99
100
101
102
    FUNCNAME("ParallelDomainBase::initParallelization()");

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

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

105
106
107
108
109
    // 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();

110
111
112
113
114
115
116
    // 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
117
#if (DEBUG != 0)
118
119
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
120
121
122
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
123

124
125
126
127
128
      if (writePartMesh > 0)
	writePartitioningMesh("part.vtu");
      else 
	MSG("Skip write part mesh!\n");
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
129
#endif
130

131
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
132

133
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
134

135
136
137
138
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
139
140
141
142
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

143
144
145
146
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

147
148
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
149
#if (DEBUG != 0)
150
    MSG("AMDiS runs in debug mode, so make some test ...\n");
151
152
153
    debug::testSortedDofs(mesh, elMap);
    ParallelDomainDbg::testInteriorBoundary(*this);
    ParallelDomainDbg::testCommonDofs(*this, true);
154
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
155

156
157
    debug::writeMesh(feSpace, -1, "macromesh");   
#endif
158

159
160
161
    // === Create periodic dof mapping, if there are periodic boundaries. ===

    createPeriodicMap();
162

163
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
164

Thomas Witkowski's avatar
Thomas Witkowski committed
165
    int globalRefinement = 0;
166
    GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinement);
Thomas Witkowski's avatar
Thomas Witkowski committed
167

Thomas Witkowski's avatar
Thomas Witkowski committed
168
    if (globalRefinement > 0) {
169
      refineManager->globalRefine(mesh, globalRefinement);
170

171
172
173
174
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "grmesh");
#endif

175
      updateLocalGlobalNumbering();
176

177
      // === Update periodic mapping, if there are periodic boundaries. ===
178

179
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
180
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219

    /// === Set DOF rank information to all matrices and vectors. ===

    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
 	if (probStat->getSystemMatrix(i, j))
 	  probStat->getSystemMatrix(i, j)->setRankDofs(isRankDof);      

      TEST_EXIT_DBG(probStat->getRHS()->getDOFVector(i))("No rhs vector!\n");
      TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))("No solution vector!\n");

      probStat->getRHS()->getDOFVector(i)->setRankDofs(isRankDof);
      probStat->getSolution()->getDOFVector(i)->setRankDofs(isRankDof);
    }

    // === Remove periodic boundary conditions in sequential problem definition. ===

    // Remove periodic boundaries in boundary manager on matrices and vectors.
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	DOFMatrix* mat = probStat->getSystemMatrix(i, j);
 	if (mat && mat->getBoundaryManager())
	  removeBoundaryCondition(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
      }

      if (probStat->getSolution()->getDOFVector(i)->getBoundaryManager())
	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat->getSolution()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));

      if (probStat->getRHS()->getDOFVector(i)->getBoundaryManager())
	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat->getRHS()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));
    }

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh,  -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
220
221
  }

222
223

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
224
  {}
225

226
227
228
  
  void ParallelDomainBase::updateDofAdmins()
  {
229
230
231
    FUNCNAME("ParallelDomainBase::updateDofAdmins()");

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
232
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
233
234
235

      // There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise,
      // it is not possible to define a first DOF hole.
236
      if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDofs.size())
237
	admin.enlargeDOFLists();
238
239
240
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
241
      for (unsigned int j = 0; j < mapLocalGlobalDofs.size(); j++)
242
243
 	admin.setDOFFree(j, false);

244
245
246
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
247
248
249
    }
  }

250

251
252
253
254
255
256
257
258
259
260
  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)
261
	("Mesh is already refined! This does not work with parallelization!\n");
262
263
264
265
266
267
268
269
270
271
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

272

Thomas Witkowski's avatar
Thomas Witkowski committed
273
  void ParallelDomainBase::synchVector(DOFVector<double> &vec)
274
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
275
    StdMpi<std::vector<double> > stdMpi(mpiComm);
276
277

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
278
	 sendIt != sendDofs.end(); ++sendIt) {
279
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
282
      
Thomas Witkowski's avatar
Thomas Witkowski committed
283
284
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
285
286
287
288

      stdMpi.send(sendIt->first, dofs);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
289
290
291
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
292

Thomas Witkowski's avatar
Thomas Witkowski committed
293
    stdMpi.startCommunication<double>(MPI_DOUBLE);
294

Thomas Witkowski's avatar
Thomas Witkowski committed
295
296
297
298
299
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      for (unsigned int i = 0; i < recvIt->second.size(); i++)
	vec[*(recvIt->second)[i]] = stdMpi.getRecvData(recvIt->first)[i];
  }
300
301


Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
304
305
306
307
308
309
310
311
312
313
314
315
  void ParallelDomainBase::synchVector(SystemVector &vec)
  {
    StdMpi<std::vector<double> > stdMpi(mpiComm);

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt) {
      std::vector<double> dofs;
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nComponents * nSendDofs);
      
      for (int i = 0; i < nComponents; i++) {
	DOFVector<double> *dofvec = vec.getDOFVector(i);
	for (int j = 0; j < nSendDofs; j++)
	  dofs.push_back((*dofvec)[*((sendIt->second)[j])]);
316
317
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
318
      stdMpi.send(sendIt->first, dofs);
319
320
321
    }

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size() * nComponents);
324

Thomas Witkowski's avatar
Thomas Witkowski committed
325
    stdMpi.startCommunication<double>(MPI_DOUBLE);
326
327

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
328
329
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
330
331

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
334
335
336
      for (int i = 0; i < nComponents; i++) {
	DOFVector<double> *dofvec = vec.getDOFVector(i);
 	for (int j = 0; j < nRecvDofs; j++)
	  (*dofvec)[*(recvIt->second)[j]] = 
	    stdMpi.getRecvData(recvIt->first)[counter++];
337
338
339
340
      }
    }
  }

341

Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
344
345
346
347
348
349
350
351
352
353
  void ParallelDomainBase::removeBoundaryCondition(BoundaryIndexMap& boundaryMap)
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


354
355
  void ParallelDomainBase::checkMeshChange()
  {
356
357
    FUNCNAME("ParallelDomainBase::checkMeshChange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
358
359
    int dim = mesh->getDim();

360
361
362
363
364
    // === If mesh has not been changed on all ranks, return. ===

    int recvAllValues = 0;
    int sendValue = static_cast<int>(mesh->getChangeIndex() != lastMeshChangeIndex);
    mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
365

366
    if (recvAllValues == 0)
367
368
      return;

369
370
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
371

372
373
374
375
376
    clock_t first = clock();
    
    do {
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
377
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
378
379

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
380
381
	if (it->rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
382
383

      for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
384
385
	if (it->rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
386

387

388
      // === Check the boundaries and adapt mesh if necessary. ===
389

390
391
392
393
394
395
396
397
      bool meshChanged = checkAndAdaptBoundary(allBound);

      // === Check on all ranks if at least one rank's mesh has changed. ===

      int sendValue = static_cast<int>(!meshChanged);
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
    } while (recvAllValues != 0);
398

399
400

#if (DEBUG != 0)
401
    debug::writeMesh(feSpace, -1, "mesh");
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
#endif

    INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n", 
		  TIME_USED(first, clock()));

    // === Because the mesh has been changed, update the DOF numbering and mappings. ===
   
    updateLocalGlobalNumbering();
  }

  
  bool ParallelDomainBase::checkAndAdaptBoundary(RankToBoundMap &allBound)
  {
    FUNCNAME("ParallelDomainBase::checkAndAdaptBoundary()");

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
422
423
424
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
425
	elCode.init(boundIt->rankObj);
426
427
428
429
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
430
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
431
    stdMpi.send(sendCodes);
432
433
434
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
435
    // === Compare received mesh structure codes. ===
436
    
437
438
    bool meshFitTogether = true;

439
440
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
441
442
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
443
      
444
445
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
446

447
448
449
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
450
451
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
452

453
454
	  //	  MSG("START WITH I = %d\n", i);

455
	  bool b = fitElementToMeshCode(recvCodes[i], 
456
457
					boundIt->rankObj.el,
					boundIt->rankObj.ithObj, 
458
					boundIt->rankObj.elType);  
459

460
461
	  if (b)
	    meshFitTogether = false;
462
	}
463
	
464
	i++;
465
466
467
      }
    }

468
    return meshFitTogether;
469
  }
470
471
472
473
474


  bool ParallelDomainBase::fitElementToMeshCode(MeshStructure &code, 
						Element *el, 
						int ithSide, 
475
						int elType)
476
  {
477
    FUNCNAME("ParallelDomainBase::fitElementToMeshCode()");
478

479
480
    TEST_EXIT_DBG(el)("No element given!\n");

481
482
483
484
485
486
487
488
    if (code.empty())
      return false;
    
    int s1 = el->getSideOfChild(0, ithSide, elType);
    int s2 = el->getSideOfChild(1, ithSide, elType);

    TEST_EXIT_DBG(s1 != -1 || s2 != -1)("This should not happen!\n");

489
490
491
    bool meshChanged = false;

    if (s1 != -1 && s2 != -1) {
492
493
494
495
496
497
498
499
500
501
      TraverseStack stack;
      ElInfo *elInfo = 
	stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
      while (elInfo && elInfo->getElement() != el) {
	elInfo = stack.traverseNext(elInfo);
      }

      //      MSG("------- START 0 ----------\n");

      meshChanged = fitElementToMeshCode2(code, stack, elInfo, ithSide, elType);
502
503
      return meshChanged;
    }
504
505

    if (el->isLeaf()) {
506
507
508
509
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

      ERROR_EXIT("NOT YET WORKING!\n");
510
      
511
512
      ElInfo *elInfo = el->getMesh()->createNewElInfo();    
      //      el->getMesh()->fillElInfo(elInfo, el, Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
513
      el->setMark(1);
514
      refineManager->refineFunction(elInfo);
515
516
      delete elInfo;

517
      meshChanged = true;
518
    }
519
520

    if (s1 != -1) {
521
522
523
524
525
526
527
528
529
530
531
      //      MSG("------- START 1 ----------\n");

      TraverseStack stack;
      ElInfo *elInfo = 
	stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);

      while (elInfo && elInfo->getElement() != el->getFirstChild()) {
	elInfo = stack.traverseNext(elInfo);
      }

      meshChanged |= fitElementToMeshCode2(code, stack, elInfo, s1, el->getChildType(elType));
532
    } else {
533
534
535
536
537
538
539
540
541
542
543
      //      MSG("------- START 2 ---------   %d \n", el->getIndex());

      TraverseStack stack;
      ElInfo *elInfo = 
	stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);

      while (elInfo && elInfo->getElement() != el->getSecondChild()) {
	elInfo = stack.traverseNext(elInfo);
      }

      meshChanged |= fitElementToMeshCode2(code, stack, elInfo, s2, el->getChildType(elType));
544
545
    }

546
    //    MSG("-------- ENDE ---------\n");
547
548

    return meshChanged;
549
550
  }

551

552
  bool ParallelDomainBase::fitElementToMeshCode2(MeshStructure &code, 
553
554
						 TraverseStack &stack,
						 ElInfo *aelInfo,
555
						 int ithSide, 
556
						 int elType)
557
  {
558
559
    FUNCNAME("ParallelDomainBase::fitElementToMeshCode2()");

560
561
562
    ElInfo *elInfo = stack.getElInfo();

    //    MSG("START EL WITH LEVEL = %d\n", elInfo->getLevel());
563
564

    bool value = false;
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
    if (!elInfo)
      return value;

    Element *el = elInfo->getElement();

    if (code.isLeafElement()) {
      int level = elInfo->getLevel();

      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);

      //      MSG("RETURN CODE LEAF: %d\n", elInfo->getLevel());
      
      return value;
    }
581

582
583
    if (!elInfo)
      return value;
584

585
586
    if (el->isLeaf()) {
      //      MSG("REFINE CODE NO LEAF!\n");
587
      el->setMark(1);
588
589
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
590
      refineManager->refineFunction(elInfo);
591
592
593
594
595
      value = true;
    }

    int s1 = el->getSideOfChild(0, ithSide, elType);
    int s2 = el->getSideOfChild(1, ithSide, elType);
596

597
    if (s1 != -1) {
598
599
      //      MSG("STEP LEFT 1\n");
      stack.traverseNext(elInfo);
600
      code.nextElement();
601
602
603
604
605
606
607
608
      value |= fitElementToMeshCode2(code, stack, elInfo, s1, el->getChildType(elType));
      elInfo = stack.getElInfo();
      //      MSG("STEP LEFT 2\n");
    } else {
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getElement() != el->getSecondChild());      
      //      MSG("STEP LEFT OMMIT\n");
609
610
    }  

611
612
613
    TEST_EXIT_DBG(elInfo->getElement() == el->getSecondChild())
      ("This should not happen!\n");

614
    if (s2 != -1) {
615
616
617
618
619
620
      //      MSG("STEP RIGHT 1\n");
      code.nextElement();
      value |= fitElementToMeshCode2(code, stack, elInfo, s2, el->getChildType(elType));
      //      MSG("STEP RIGHT 2\n");
    } else {
      int level = elInfo->getLevel();
621

622
623
624
625
626
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
      //      MSG("STEP RIGHT OMMIT\n");
    }
627

628
629
630
    return value;
  }

631
632
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
633
  {    
634
    int vecSize = data.size();
635
    SerUtil::serialize(out, vecSize);
636
637
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
638
      SerUtil::serialize(out, dofIndex);
639
640
641
642
643
644
645
646
647
648
    }
  }


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

    int vecSize = 0;
649
    SerUtil::deserialize(in, vecSize);
650
651
652
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
653
      SerUtil::deserialize(in, dofIndex);
654
655
656
657
658
659
660
661
662
663
664
665

      TEST_EXIT_DBG(dofMap.count(dofIndex) != 0)
	("Dof index could not be deserialized correctly!\n");

      data[i] = dofMap[dofIndex];
    }
  }


  void ParallelDomainBase::serialize(std::ostream &out, RankToDofContainer &data)
  {
    int mapSize = data.size();
666
    SerUtil::serialize(out, mapSize);
667
668
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
669
      SerUtil::serialize(out, rank);
670
671
672
673
674
675
676
677
678
      serialize(out, it->second);
    }
  }

  
  void ParallelDomainBase::deserialize(std::istream &in, RankToDofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    int mapSize = 0;
679
    SerUtil::deserialize(in, mapSize);
680
681
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
682
      SerUtil::deserialize(in, rank);
683
684
685
686
      deserialize(in, data[rank], dofMap);      
    }
  }

687

688
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
689
  {
690
    FUNCNAME("ParallelDomainBase::setElemWeights()");
691

692
    double localWeightSum = 0.0;
693
    elemWeights.clear();
694
695
696
697
698
699
700
701
702
703
704
705
706
707
      
    std::string filename = "";
    GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

      std::ifstream infile;
      infile.open(filename.c_str(), std::ifstream::in);
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;
708

709
710
711
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
712

713
714
715
716
717
718
719
720
721
722
723
      infile.close();
    } else {
           
      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));
724
	
725
726
727
728
729
730
731
732
733
734
	if (partitionData && partitionData->getPartitionStatus() == IN) {
	  int elNum = -1;
	  if (partitionData->getLevel() == 0)
	    elNum = element->getIndex();
	  
	  TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
	  if (element->isLeaf()) {
	    elemWeights[elNum] += 1.0;
	    localWeightSum += 1.0;
	  }
735
	}
736
737
	
	elInfo = stack.traverseNext(elInfo);
738
739
740
741
742
743
      }
    }

    return localWeightSum;
  }

744

745
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
746
  {
747
748
    FUNCNAME("ParallelDomainBase::partitionMesh()");

749
750
751
752
753
754
755
756
757
758
759
760
    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);
  }

761

Thomas Witkowski's avatar
Thomas Witkowski committed
762
763
764
765
766
767
768
769
770
  void ParallelDomainBase::solveInitialProblem(AdaptInfo *adaptInfo)
  {     
    if (timeIF)
      timeIF->solveInitialProblem(adaptInfo);
    
    synchVector(*(probStat->getSolution()));
  }


771
  void ParallelDomainBase::createInteriorBoundaryInfo()
772
  {
773
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
774

Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
    // === First, create the interior boundaries based on macro element's  ===
    // === neighbour informations.                                         ===

    createMacroElementInteriorBoundaryInfo();

    // === Second, search the whole mesh for interior boundaries that consists of ===
    // === substructures of the macro elements.                                   ===

    createSubstructureInteriorBoundaryInfo();


    // === Once we have this information, we must care about the order of the atomic ===
    // === bounds in the three boundary handling object. Eventually all the bound-   ===
    // === aries have to be in the same order on both ranks that share the bounday.  ===

    StdMpi<std::vector<AtomicBoundary> > stdMpi(mpiComm);
    stdMpi.send(myIntBoundary.boundary);
    stdMpi.recv(otherIntBoundary.boundary);
    stdMpi.startCommunication<int>(MPI_INT);
   
    // === The information about all neighbouring boundaries has been received. So ===
    // === the rank tests if its own atomic boundaries are in the same order. If   ===
    // === not, the atomic boundaries are swaped to the correct order.             ===

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

      // === 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.

	BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj;

	if ((rankIt->second)[j].neighObj != recvedBound) {
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
 	    if ((rankIt->second)[k].neighObj == recvedBound)
	      break;

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

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

    // === Do the same for the periodic boundaries. ===

    if (periodicBoundary.boundary.size() > 0) {
      stdMpi.clear();

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

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

	if (rankIt->first < mpiRank)
	  sendBounds.boundary[rankIt->first] = rankIt->second;
	else
	  recvBounds.boundary[rankIt->first] = rankIt->second;
      }

      stdMpi.send(sendBounds.boundary);
      stdMpi.recv(recvBounds.boundary);
      stdMpi.startCommunication<int>(MPI_INT);

      for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
	   rankIt != periodicBoundary.boundary.end(); ++rankIt) {
	if (rankIt->first <= mpiRank)
	  continue;
	  
	for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
	  
	  BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj;
	  
	  if (periodicBoundary.boundary[rankIt->first][j].neighObj != recvedBound) {    
	    int k = j + 1;	    
	    for (; k < static_cast<int>(rankIt->second.size()); k++)
	      if (periodicBoundary.boundary[rankIt->first][k].neighObj == recvedBound)
		break;
	    
	    // The element must always be found, because the list is just in 
	    // another order.
	    TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
	      ("Should never happen!\n");
	    
	    // Swap the current with the found element.
	    AtomicBoundary tmpBound = (rankIt->second)[k];
	    (rankIt->second)[k] = (rankIt->second)[j];
	    (rankIt->second)[j] = tmpBound;	
	  }  	  
	}
      }     
    } // periodicBoundary.boundary.size() > 0
  }


  void ParallelDomainBase::createMacroElementInteriorBoundaryInfo()
  {
    FUNCNAME("ParallelDomainBase::createMacroElementInteriorBoundaryInfo()");

886
887
    int nNeighbours = mesh->getGeo(NEIGH);
    int dim = mesh->getDim();
888
    GeoIndex subObj = INDEX_OF_DIM(dim - 1, dim);
889
890
891

    // === First, traverse the mesh and search for all elements that have an  ===
    // === boundary with an element within another partition.                 ===
Thomas Witkowski's avatar
Thomas Witkowski committed
892

893
    TraverseStack stack;
894
895
896
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
897
898
899
900
    while (elInfo) {
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
901

902
      // Check, if the element is within rank's partition.
903
      if (partitionData->getPartitionStatus() == IN) {
904
	for (int i = 0; i < nNeighbours; i++) {
905
906
907
908
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
909
910
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->
						getElementData(PARTITION_ED));
911

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

914
915
916
917
918
	    // We have found an element that is in rank's partition, but has a 
	    // neighbour outside of the rank's partition.

	    // === Create information about the boundary between the two elements. ===

919
	    AtomicBoundary bound;	    	    
920
921
	    bound.rankObj.el = element;
	    bound.rankObj.elIndex = element->getIndex();
922
	    bound.rankObj.elType = elInfo->getType();
923
	    bound.rankObj.subObj = subObj;
924
	    bound.rankObj.ithObj = i;
925
926

	    bound.neighObj.el = elInfo->getNeighbour(i);
927
	    bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex();
928
	    bound.neighObj.elType = -1;
929
	    bound.neighObj.subObj = subObj;
930
	    bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i);
931
	    
932
	    // Get rank number of the neighbouring element.
933
934
	    int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];

935
936
937
938
939
940
941
942
	    // === Add the boundary information object to the corresponding overall ===
	    // === boundary. There are three possibilities: if the boundary is a    ===
	    // === periodic boundary, just add it to \ref periodicBounadry. Here it ===
	    // === does not matter which rank is responsible for this boundary.     ===
	    // === Otherwise, the boundary is added either to \ref myIntBoundary or ===
	    // === to \ref otherIntBoundary. It dependes on which rank is respon-   ===
	    // === sible for this boundary.                                         ===

943
	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(subObj, i))) {	      
944
945
946
947
948
	      // We have found an element that is at an interior, periodic boundary.
	      AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
	      b = bound;
	    } else {
	      // We have found an element that is at an interior, non-periodic boundary.
949
950
951
952
953
954
955
956
957
	      if (mpiRank > otherElementRank) {
		AtomicBoundary& b = myIntBoundary.getNewAtomic(otherElementRank);
		b = bound;
		b.rankObj.setReverseMode(b.neighObj, feSpace);
	      } else {
		AtomicBoundary& b = otherIntBoundary.getNewAtomic(otherElementRank);
		b = bound;	 
		b.neighObj.setReverseMode(b.rankObj, feSpace);
	      }	      
958
	    }
959
960
961
962
963
964
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
965
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
966
967


Thomas Witkowski's avatar
Thomas Witkowski committed
968
969
970
  void ParallelDomainBase::createSubstructureInteriorBoundaryInfo()
  {
    FUNCNAME("ParallelDomainBase::createSubstructureInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
971

972
973
974
    // === Seach for all vertices/edges, which are part of an interior boundary,  ===
    // === but are not part of the interior boundaries that are created based on  ===
    // === the information of macro elements.                                     ===
Thomas Witkowski's avatar
Thomas Witkowski committed
975
976

    int dim = mesh->getDim();
977
    const BasisFunction *basFcts = feSpace->getBasisFcts();
978
    std::vector<DegreeOfFreedom> localIndices(basFcts->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
979
980
981

    // Maps each DOF in the whole domain to the rank that ownes this DOF.
    std::map<DegreeOfFreedom, int> dofOwner;
982
    // Maps each DOF in ranks partition to the element object that contains this DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
983
    std::map<DegreeOfFreedom, BoundaryObject> rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
984

985
986
    // Maps each edge in the whole domain to the rank that ownes this edge.
    std::map<GlobalEdge, int> edgeOwner;
987
    // Maps each edge in ranks partition to the element object that contains this edge.
988
    std::map<GlobalEdge, BoundaryObject> rankEdges;
989

990

Thomas Witkowski's avatar
Thomas Witkowski committed
991
    // === Traverse whole mesh and fill the maps defined above.                    ===
992

Thomas Witkowski's avatar
Thomas Witkowski committed
993
994
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
995
996
997
998
999
1000
    while (elInfo) {
      Element *el = elInfo->getElement();
      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);

      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(el->getElementData(PARTITION_ED));
1001
1002
          
      // In 2d and 3d, traverse all vertices of current element.
1003
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
	DegreeOfFreedom dof0 = localIndices[i];
	// Update the owner of the current vertex DOF.
	dofOwner[dof0] = max(dofOwner[dof0], partitionVec[el->getIndex()]);

	// If the DOF is part of an element that is within ranks partition, add it
	// to the set of all ranks vertex DOFs.
	if (partitionData && partitionData->getPartitionStatus() == IN)
	  rankDofs[dof0] = BoundaryObject(el, elInfo->getType(), VERTEX, i);
      }      

Thomas Witkowski's avatar
Thomas Witkowski committed
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
      // In 3d, traverse all edges of current element.
      if (dim == 3) {
	for (int i = 0; i < 6; i++) {
	  DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(i, 0)];
	  DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(i, 1)];
	  GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));

	  // Update the owner of the current edge.
	  edgeOwner[edge] = max(edgeOwner[edge], partitionVec[el->getIndex()]);
	  
1024
1025
	  // If the edge is part of an element that is within ranks partition, add it
	  // to the set of all ranks edges.
Thomas Witkowski's avatar
Thomas Witkowski committed
1026
1027
	  if (partitionData && partitionData->getPartitionStatus() == IN)
	    rankEdges[edge] = BoundaryObject(el, elInfo->getType(), EDGE, i);
1028
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
1029
1030
      }

1031
1032
1033
      elInfo = stack.traverseNext(elInfo);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1034

1035
    // === Create a set of all edges and vertices at ranks interior boundaries. ===
1036
1037
1038

    // Stores all edges at rank's interior boundaries.
    std::set<GlobalEdge> rankBoundaryEdges;
Thomas Witkowski's avatar
Thomas Witkowski committed
1039
    // Stores all vertices at rank's interior boundaries.
1040
1041
    std::set<DegreeOfFreedom> rankBoundaryDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
1042
1043
   
    // First, traverse the rank owned elements at the interior boundaries.
Thomas Witkowski's avatar
Thomas Witkowski committed
1044
    for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
1045
      Element *el = it->rankObj.el;
Thomas Witkowski's avatar
Thomas Witkowski committed
1046
      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
1047
      
1048
1049
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
	int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, i);
1050
1051
1052
1053
1054
1055
1056
1057
	DegreeOfFreedom dof = localIndices[dofNo];
	
	if (dofOwner[dof] > mpiRank)
	  it->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
	else
	  rankBoundaryDofs.insert(dof);
      }    

Thomas Witkowski's avatar
Thomas Witkowski committed
1058
      if (dim == 3) {
1059
1060
	for (int i = 0; i < 3; i++) {
	  int edgeNo = el->getEdgeOfFace(it->rankObj.ithObj, i);
Thomas Witkowski's avatar
Thomas Witkowski committed
1061
1062
1063
1064
1065
1066
1067
1068
1069
	  DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)];
	  DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)];
	  GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));	
	  
	  // If the edge at rank's interior boundarie has a higher owner rank, than
	  // we have to remove this edge from the corresponding boundary element.
	  // Otherwise, it is part of the interior boundary and we add it to the set
	  // rankBoundaryEdges.
	  if (edgeOwner[edge] > mpiRank)
1070
	    it->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
1071
	  else
Thomas Witkowski's avatar
Thomas Witkowski committed
1072
	    rankBoundaryEdges.insert(edge);
1073
	}
1074
      }     
Thomas Witkowski's avatar
Thomas Witkowski committed
1075
1076
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1077
1078
    // Now, do the same with all other elements at the interior boundaries.
    for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) {
1079
      Element *el = it->rankObj.el;
Thomas Witkowski's avatar
Thomas Witkowski committed
1080
1081
      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
      
1082
1083
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
	int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, i);
1084
1085
1086
1087
1088
1089
1090
1091
1092

	DegreeOfFreedom dof = localIndices[dofNo];
	
	if (dofOwner[dof] > it.getRank())
	  it->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
	else
	  rankBoundaryDofs.insert(dof);
      }	      

Thomas Witkowski's avatar
Thomas Witkowski committed
1093
      if (dim == 3) {
1094
1095
	for (int i = 0; i < 3; i++) {
	  int edgeNo = el->getEdgeOfFace(it->rankObj.ithObj, i);
1096
1097
1098
	  DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)];
	  DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)];
	  GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));	
1099
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
1100
	  if (edgeOwner[edge] > it.getRank())
1101
	    it->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
1102
	  else
Thomas Witkowski's avatar
Thomas Witkowski committed
1103
1104
	    rankBoundaryEdges.insert(edge);	    
	}
1105
      }     
1106
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1107

1108
    // === Create the new interior boundaries consisting only of edges. These   ===
1109
    // === boundaries are created on that ranks, which do not own the boundary  ===
1110
1111
    // === but are on the other side of the edge. Then, these ranks inform the  ===
    // === owner of the edge.                                                   ===
1112
1113
1114
1115
1116

    // Maps that stores to each rank number the global edges this rank has an 
    // interior boundary with. 
    std::map<int, std::vector<GlobalEdge> > recvEdges;
    // Stores to the map above the corresponding element objects that include the edge.
Thomas Witkowski's avatar
Thomas Witkowski committed
1117
    std::map<int, std::vector<BoundaryObject> > recvObjs;