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

38

39
40
41
42
43
44
45
46
  MeshDistributor::MeshDistributor(std::string str)
    : probStat(0),
      name(str),
      feSpace(NULL),
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
47
      initialPartitionMesh(true),
48
      nRankDofs(0),
49
      nOverallDofs(0),
50
      rstart(0),
51
      deserialized(false),
52
53
      lastMeshChangeIndex(0),
      macroElementStructureConsisten(false)
54
  {
55
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
56

57
58
59
60
61
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
62

63
  void MeshDistributor::initParallelization(AdaptInfo *adaptInfo)
64
  {
65
    FUNCNAME("MeshDistributor::initParallelization()");
66
67
68
69

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

70
71
    TEST_EXIT(mesh)("No mesh has been defined for mesh distribution!\n");

72
73
    // If the problem has been already read from a file, we do not need to do anything.
    if (deserialized)
74
      return;
75
    
76

77
78
79
80
81
    // 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();

82

83
84
85
86
87
88
89
    // 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);   

90

Thomas Witkowski's avatar
Thomas Witkowski committed
91
#if (DEBUG != 0)
92
93
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
94
95
96
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
97

98
99
100
101
102
      if (writePartMesh > 0)
	writePartitioningMesh("part.vtu");
      else 
	MSG("Skip write part mesh!\n");
    }
103
104

    ParallelDomainDbg::testAllElements(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
105
#endif
106

107

108
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
109

110
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
111

112
113
114
115
#if (DEBUG != 0)
    ParallelDomainDbg::printBoundaryInfo(*this);
#endif

116

117
118
119
120
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

121

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

    removeMacroElements();
125
126
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
127

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

    updateDofAdmins();

132

133
134
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
135
#if (DEBUG != 0)
136
    MSG("AMDiS runs in debug mode, so make some test ...\n");
137
138
139
    debug::testSortedDofs(mesh, elMap);
    ParallelDomainDbg::testInteriorBoundary(*this);
    ParallelDomainDbg::testCommonDofs(*this, true);
140
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
141

142
    debug::writeMesh(feSpace, -1, "macro_mesh");   
143
#endif
144

145

146
147
148
    // === Create periodic dof mapping, if there are periodic boundaries. ===

    createPeriodicMap();
149

150

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

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

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

159
#if (DEBUG != 0)
160
      debug::writeMesh(feSpace, -1, "gr_mesh");
161
162
#endif

163
      updateLocalGlobalNumbering();
164
      
165
      // === Update periodic mapping, if there are periodic boundaries. ===
166
      
167
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
168
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
169

170

Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
    /// === Set DOF rank information to all matrices and vectors. ===

173
174
175
176
177
178
    for (unsigned int i = 0; i < probStat.size(); i++) {
      int nComponents = probStat[i]->getNumComponents();
      for (int j = 0; j < nComponents; j++) {
	for (int k = 0; k < nComponents; k++)
	  if (probStat[i]->getSystemMatrix(j, k))
	    probStat[i]->getSystemMatrix(j, k)->setRankDofs(isRankDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
179

180
181
182
183
184
185
	TEST_EXIT_DBG(probStat[i]->getRhs()->getDOFVector(j))("No RHS vector!\n");
	TEST_EXIT_DBG(probStat[i]->getSolution()->getDOFVector(j))("No solution vector!\n");
	
	probStat[i]->getRhs()->getDOFVector(j)->setRankDofs(isRankDof);
	probStat[i]->getSolution()->getDOFVector(j)->setRankDofs(isRankDof);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
186
187
    }

188

Thomas Witkowski's avatar
Thomas Witkowski committed
189
190
191
    // === Remove periodic boundary conditions in sequential problem definition. ===

    // Remove periodic boundaries in boundary manager on matrices and vectors.
192
193
194
    for (unsigned int i = 0; i < probStat.size(); i++) {
      int nComponents = probStat[i]->getNumComponents();

Thomas Witkowski's avatar
Thomas Witkowski committed
195
      for (int j = 0; j < nComponents; j++) {
196
197
198
199
200
201
202
203
204
	for (int k = 0; k < nComponents; k++) {
	  DOFMatrix* mat = probStat[i]->getSystemMatrix(j, k);
	  if (mat && mat->getBoundaryManager())
	    removeBoundaryCondition(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
	}
	
	if (probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager())
	  removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
	
205
	if (probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager())
206
	  removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
Thomas Witkowski's avatar
Thomas Witkowski committed
207
208
209
210
211
212
213
214
215
216
      }
    }

    // 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);
    }    
217
218
  }

219

220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

    if (feSpace != NULL) {
      std::vector<FiniteElemSpace*> feSpaces = probVec->getFeSpaces();
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	TEST_EXIT(feSpace == feSpaces[i])
	  ("Parallelizaton is not supported for multiple FE spaces!\n");
      }
    } else {
      feSpace = probVec->getFeSpace(0);
      mesh = feSpace->getMesh();
      info = probVec->getInfo();
      
      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");
      
      switch (mesh->getDim()) {
      case 2:
	refineManager = new RefinementManager2d();
	break;
      case 3:
	refineManager = new RefinementManager3d();
	break;
      default:
	ERROR_EXIT("This should not happen for dim = %d!\n", mesh->getDim());
      }

      partitioner = new ParMetisPartitioner(mesh, &mpiComm);
    }

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

    int readSerialization = 0;
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", &readSerialization);
    if (readSerialization) {
      ERROR_EXIT("Must be reimplemented!\n");
#if 0      
      std::string filename = "";
      GET_PARAMETER(0, probVec->getName() + "->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();
#endif
    }

    probStat.push_back(probVec);
  }


  void MeshDistributor::exitParallelization(AdaptInfo *adaptInfo)
280
  {}
281

282
  
283
  void MeshDistributor::updateDofAdmins()
284
  {
285
    FUNCNAME("MeshDistributor::updateDofAdmins()");
286
287

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
288
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
289
290
291

      // There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise,
      // it is not possible to define a first DOF hole.
292
      if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDofs.size())
293
	admin.enlargeDOFLists();
294
295
296
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
297
      for (unsigned int j = 0; j < mapLocalGlobalDofs.size(); j++)
298
299
 	admin.setDOFFree(j, false);

300
301
302
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
303
304
305
    }
  }

306

307
  void MeshDistributor::testForMacroMesh()
308
  {
309
    FUNCNAME("MeshDistributor::testForMacroMesh()");
310
311
312
313
314
315
316

    int nMacroElements = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      TEST_EXIT(elInfo->getLevel() == 0)
317
	("Mesh is already refined! This does not work with parallelization!\n");
318
319
320
321
322
323
324
325
326
327
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

328

329
  void MeshDistributor::synchVector(DOFVector<double> &vec)
330
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
331
    StdMpi<std::vector<double> > stdMpi(mpiComm);
332
333

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
334
	 sendIt != sendDofs.end(); ++sendIt) {
335
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
338
      
Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
341
342
343
344

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

Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
347
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
348

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

Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353
354
355
    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];
  }
356
357


358
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
359
  {
360
    int nComponents = vec.getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
363
364
365
366
367
368
369
370
371
372
    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])]);
373
374
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
375
      stdMpi.send(sendIt->first, dofs);
376
377
378
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
382
    stdMpi.startCommunication<double>(MPI_DOUBLE);
383
384

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
387
388

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
391
392
393
      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++];
394
395
396
397
      }
    }
  }

398

399
  void MeshDistributor::removeBoundaryCondition(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
402
403
404
405
406
407
408
409
410
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


411
  void MeshDistributor::checkMeshChange()
412
  {
413
    FUNCNAME("MeshDistributor::checkMeshChange()");
414

415
416
417
418
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "before_check_mesh");
#endif

419
420
421
422
423
    // === 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);
424

425
    if (recvAllValues == 0)
426
427
      return;

428
429
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
430

431
432
433
434
435
    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.
436
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
437
438

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
439
	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
440
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
441
442

      for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
443
	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
444
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
445

446

447
      // === Check the boundaries and adapt mesh if necessary. ===
448

449
450
451
452
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

453
454
455
456
457
458
459
      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);
460
461
462
463

#if (DEBUG != 0)
      MSG("Mesh changed on %d ranks!\n", recvAllValues);
#endif
464
    } while (recvAllValues != 0);
465

466
467

#if (DEBUG != 0)
468
    debug::writeMesh(feSpace, -1, "mesh");
469
470
471
472
473
474
#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. ===
475

476
    mesh->dofCompress();
477
478
479
480
    updateLocalGlobalNumbering();
  }

  
481
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
482
  {
483
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
484
485
486
487
488
489

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
490
491
492
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
493
	elCode.init(boundIt->rankObj);
494
495
496
497
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
498
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
499
    stdMpi.send(sendCodes);
500
501
502
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
503
    // === Compare received mesh structure codes. ===
504
    
505
506
    bool meshFitTogether = true;

507
508
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
509
510
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
511
      
512
513
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
514

515
516
517
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
518
519
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
520
521

	  bool b = fitElementToMeshCode(recvCodes[i], 
522
					boundIt->rankObj.el,
523
					boundIt->rankObj.subObj,
524
					boundIt->rankObj.ithObj, 
525
526
					boundIt->rankObj.elType,
					boundIt->rankObj.reverseMode);  
527

528
529
	  if (b)
	    meshFitTogether = false;
530
	}
531
	
532
	i++;
533
534
535
      }
    }

536
    return meshFitTogether;
537
  }
538
539


540
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
541
542
543
544
545
					     Element *el, 
					     GeoIndex subObj,
					     int ithObj, 
					     int elType,
					     bool reverseMode)
546
  {
547
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");
548

549
550
    TEST_EXIT_DBG(el)("No element given!\n");

551
552
    if (code.empty())
      return false;
553
554
555

    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
556
557
558

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

559
    bool meshChanged = false;
560
561
562
563
564
565
566
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

    if (reverseMode) {
      std::swap(s1, s2);
      traverseFlag |= Mesh::CALL_REVERSE_MODE;
    }
567
568

    if (s1 != -1 && s2 != -1) {
569
      TraverseStack stack;
570
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
571
572
573
574
      while (elInfo && elInfo->getElement() != el) {
	elInfo = stack.traverseNext(elInfo);
      }

575
      meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
576
577
      return meshChanged;
    }
578
579

    if (el->isLeaf()) {
580
581
582
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

583
584
585
      meshChanged = true;

      TraverseStack stack;
586
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
587
588
589
590
591
592
593

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

      TEST_EXIT_DBG(elInfo)("This should not happen!\n");

594
      el->setMark(1);
595
596
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
597
      refineManager->refineFunction(elInfo);
598
    }
599

600
601
602
603
604
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode)
      std::swap(child0, child1);

605
    if (s1 != -1) {
606
      TraverseStack stack;
607
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
608

609
      while (elInfo && elInfo->getElement() != child0) {
610
611
612
	elInfo = stack.traverseNext(elInfo);
      }

613
      meshChanged |= 
614
	fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
615
    } else {
616
      TraverseStack stack;
617
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
618

619
      while (elInfo && elInfo->getElement() != child1) {
620
621
622
	elInfo = stack.traverseNext(elInfo);
      }

623
      meshChanged |= 
624
	fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
625
626
627
    }

    return meshChanged;
628
629
  }

630

631
  bool MeshDistributor::fitElementToMeshCode2(MeshStructure &code, 
632
633
634
635
636
					      TraverseStack &stack,
					      GeoIndex subObj,
					      int ithObj, 
					      int elType,
					      bool reverseMode)
637
  {
638
    FUNCNAME("MeshDistributor::fitElementToMeshCode2()");
639

640
641
    ElInfo *elInfo = stack.getElInfo();

642
    bool value = false;
643
644
645
646
647
648
649
650
651
652
653
    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);
654
     
655
656
      return value;
    }
657

658
659
    if (!elInfo)
      return value;
660

661
    if (el->isLeaf()) {
662
      el->setMark(1);
663
664
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
665
      refineManager->refineFunction(elInfo);
666
667
668
      value = true;
    }

669
670
    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
671
672
673
674
675
676
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
      std::swap(s1, s2);
      std::swap(child0, child1);
    }
677

678
    if (s1 != -1) {
679
      stack.traverseNext(elInfo);
680
      code.nextElement();
681
      value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
682
683
684
685
      elInfo = stack.getElInfo();
    } else {
      do {
	elInfo = stack.traverseNext(elInfo);
686
      } while (elInfo && elInfo->getElement() != child1);      
687
688
    }  

689
    TEST_EXIT_DBG(elInfo->getElement() == child1)("This should not happen!\n");
690

691
    if (s2 != -1) {
692
      code.nextElement();
693
      value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
694
695
    } else {
      int level = elInfo->getLevel();
696

697
698
699
700
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
701

702
703
704
    return value;
  }

705
  
706
  void MeshDistributor::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
707
  {    
708
    int vecSize = data.size();
709
    SerUtil::serialize(out, vecSize);
710
711
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
712
      SerUtil::serialize(out, dofIndex);
713
714
715
716
    }
  }


717
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
718
				    std::map<int, const DegreeOfFreedom*> &dofMap)
719
  {
720
    FUNCNAME("MeshDistributor::deserialize()");
721
722

    int vecSize = 0;
723
    SerUtil::deserialize(in, vecSize);
724
725
726
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
727
      SerUtil::deserialize(in, dofIndex);
728
729
730
731
732
733
734
735
736

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

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


737
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)
738
739
  {
    int mapSize = data.size();
740
    SerUtil::serialize(out, mapSize);
741
742
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
743
      SerUtil::serialize(out, rank);
744
745
746
747
748
      serialize(out, it->second);
    }
  }

  
749
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
750
				    std::map<int, const DegreeOfFreedom*> &dofMap)
751
752
  {
    int mapSize = 0;
753
    SerUtil::deserialize(in, mapSize);
754
755
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
756
      SerUtil::deserialize(in, rank);
757
758
759
760
      deserialize(in, data[rank], dofMap);      
    }
  }

761

762
  double MeshDistributor::setElemWeights(AdaptInfo *adaptInfo) 
763
  {
764
    FUNCNAME("MeshDistributor::setElemWeights()");
765

766
    double localWeightSum = 0.0;
767
    elemWeights.clear();
768
769
770
771
772
773
774
775
776
777
778
779
780
781
      
    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;
782

783
784
785
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
786

787
      infile.close();
788
    } else {           
789
      TraverseStack stack;
790
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
791
      while (elInfo) {
792
793
794
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
	localWeightSum++;

795
	elInfo = stack.traverseNext(elInfo);
796
797
798
799
800
801
      }
    }

    return localWeightSum;
  }

802

803
  void MeshDistributor::partitionMesh(AdaptInfo *adaptInfo)
804
  {
805
    FUNCNAME("MeshDistributor::partitionMesh()");
806

807
808
809
810
811
812
813
814
815
816
817
818
    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);
  }

819

820
  void MeshDistributor::createInteriorBoundaryInfo()
821
  {
822
    FUNCNAME("MeshDistributor::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
823

Thomas Witkowski's avatar
Thomas Witkowski committed
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
    // === 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);
843

Thomas Witkowski's avatar
Thomas Witkowski committed
844
845
846
847
848
849
850
851
    // === 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 ===
852
      // === indices. Now, we have to sort the corresponding list in this rank to   ===
Thomas Witkowski's avatar
Thomas Witkowski committed
853
      // === get the same order.                                                    ===
854
     
Thomas Witkowski's avatar
Thomas Witkowski committed
855
      for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
856

Thomas Witkowski's avatar
Thomas Witkowski committed
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
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
	// 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
  }


932
  void MeshDistributor::createMacroElementInteriorBoundaryInfo()
Thomas Witkowski's avatar
Thomas Witkowski committed
933
  {
934
    FUNCNAME("MeshDistributor::createMacroElementInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
935

936
937
    int nNeighbours = mesh->getGeo(NEIGH);
    int dim = mesh->getDim();
938
    GeoIndex subObj = INDEX_OF_DIM(dim - 1, dim);
939
940
941

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

943
    TraverseStack stack;
944
945
946
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
947
948
949
    while (elInfo) {
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
950
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
951

952
      // Check, if the element is within rank's partition.
953
      if (partitionData->getPartitionStatus() == IN) {
954
	for (int i = 0; i < nNeighbours; i++) {
955
956
957
958
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
959
960
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->
						getElementData(PARTITION_ED));
961

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

964
965
966
967
968
	    // 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. ===

969
	    AtomicBoundary bound;	    	    
970
971
	    bound.rankObj.el = element;
	    bound.rankObj.elIndex = element->getIndex();