MeshDistributor.cc 65 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
6
#include <boost/filesystem.hpp>

7
8
#include "parallel/MeshDistributor.h"
#include "parallel/ParallelDebug.h"
9
#include "parallel/StdMpi.h"
10
11
12
13
14
15
16
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
17
18
#include "DOFMatrix.h"
#include "DOFVector.h"
19
#include "SystemVector.h"
20
#include "VtkWriter.h"
21
#include "ElementDofIterator.h"
22
23
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
24
#include "ElementFileWriter.h"
25
#include "VertexVector.h"
26
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
27
28
#include "ProblemVec.h"
#include "ProblemInstat.h"
29
#include "Debug.h"
30

31
32
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
33
  using boost::lexical_cast;
34
  using namespace boost::filesystem;
Thomas Witkowski's avatar
Thomas Witkowski committed
35

36
37
38
39
40
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

41

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

61
62
63
64
65
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
66

67
  void MeshDistributor::initParallelization(AdaptInfo *adaptInfo)
68
  {
69
    FUNCNAME("MeshDistributor::initParallelization()");
70
71
72
73

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

74
75
    TEST_EXIT(mesh)("No mesh has been defined for mesh distribution!\n");

76
77
78
    // If the problem has been already read from a file, we need only to set
    // isRankDofs to all matrices and rhs vector and to remove periodic 
    // boundary conditions (if there are some).
79
    if (deserialized) {
80
      setRankDofs();
81
      removePeriodicBoundaryConditions();
82

83
      return;
84
    }
85
    
86

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

100

Thomas Witkowski's avatar
Thomas Witkowski committed
101
#if (DEBUG != 0)
102
103
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
104
105
106
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
107

108
109
110
111
112
      if (writePartMesh > 0)
	writePartitioningMesh("part.vtu");
      else 
	MSG("Skip write part mesh!\n");
    }
113

114
    ParallelDebug::testAllElements(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
115
#endif
116

117

118
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
119

120
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
121

122
#if (DEBUG != 0)
123
    ParallelDebug::printBoundaryInfo(*this);
124
125
#endif

126

127
128
129
130
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

131

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
132
133
134
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
135
136
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
137

138
139
140
141
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

142

143
144
    // === If in debug mode, make some tests. ===

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

152
    debug::writeMesh(feSpace, -1, "macro_mesh");   
153
#endif
154

155

156
157
    // === Create periodic dof mapping, if there are periodic boundaries. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
158
    createPeriodicMap();    
159

160

161
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
162

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

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

169
#if (DEBUG != 0)
170
      debug::writeMesh(feSpace, -1, "gr_mesh");
171
172
#endif

173
      updateLocalGlobalNumbering();
174
      
175
      // === Update periodic mapping, if there are periodic boundaries. ===
176
      
177
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
178
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
179

180

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

183
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
184

185

Thomas Witkowski's avatar
Thomas Witkowski committed
186
187
    // === Remove periodic boundary conditions in sequential problem definition. ===

188
    removePeriodicBoundaryConditions();
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
220
221
222
223
224
225
226
227
  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;
228
229
230
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
231
232
233
234
235
      std::string filename = "";
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
236
237
238
239
240

      int tsModulo = -1;
      GET_PARAMETER(0, probVec->getName() + "->output->write every i-th timestep", 
		    "%d", &tsModulo);
      
241
      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
242
243
      writeSerializationFile = true;
    }    
244
245

    int readSerialization = 0;
246
247
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", 
		  &readSerialization);
248
249
250
251
    if (readSerialization) {
      std::string filename = "";
      GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
      filename += ".p" + lexical_cast<std::string>(mpiRank);
252
      MSG("Start deserialization with %s\n", filename.c_str());
253
      std::ifstream in(filename.c_str());
254
255
256
257

      TEST_EXIT(!in.fail())("Could not open deserialization file: %s\n",
			    filename.c_str());

258
      probVec->deserialize(in);
259
      in.close();
260
261
      MSG("Deserialization from file: %s\n", filename.c_str());

262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
      filename = "";
      GET_PARAMETER(0, name + "->input->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
      std::string rankFilename = filename + ".p" + lexical_cast<std::string>(mpiRank);
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
278
279
280
281
282
283
284
    }

    probStat.push_back(probVec);
  }


  void MeshDistributor::exitParallelization(AdaptInfo *adaptInfo)
285
  {}
286

287
  
288
  void MeshDistributor::updateDofAdmins()
289
  {
290
    FUNCNAME("MeshDistributor::updateDofAdmins()");
291
292

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
293
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
294
295
296

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

305
306
307
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
308
309
310
    }
  }

311

312
  void MeshDistributor::testForMacroMesh()
313
  {
314
    FUNCNAME("MeshDistributor::testForMacroMesh()");
315
316
317
318
319
320
321

    int nMacroElements = 0;

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

      elInfo = stack.traverseNext(elInfo);
    }

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

333

334
  void MeshDistributor::synchVector(DOFVector<double> &vec)
335
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    StdMpi<std::vector<double> > stdMpi(mpiComm);
337
338

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
350
351
352
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
353

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

Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
358
359
360
    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];
  }
361
362


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

Thomas Witkowski's avatar
Thomas Witkowski committed
380
      stdMpi.send(sendIt->first, dofs);
381
382
383
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
387
    stdMpi.startCommunication<double>(MPI_DOUBLE);
388
389

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
390
391
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
392
393

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
394
395
396
397
398
      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++];
399
400
401
402
      }
    }
  }

403

404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
  void MeshDistributor::setRankDofs()
  {
    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);

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


423
424
425
426
427
428
429
430
431
432
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
    // Remove periodic boundaries in boundary manager on matrices and vectors.
    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++) {
	  DOFMatrix* mat = probStat[i]->getSystemMatrix(j, k);
	  if (mat && mat->getBoundaryManager())
433
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
	}
	
	if (probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
	
	if (probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat[i]->getRhs()->getDOFVector(j)->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);
    }    
  }


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
455
456
457
458
459
460
461
462
463
464
465
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


466
  void MeshDistributor::checkMeshChange()
467
  {
468
    FUNCNAME("MeshDistributor::checkMeshChange()");
469

470
471
472
473
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "before_check_mesh");
#endif

474
475
476
477
478
    // === 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);
479

480
    if (recvAllValues == 0)
481
482
      return;

483
484
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
485

486
487
488
489
490
    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.
491
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
492
493

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

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

501
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
502
503
 	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
 	  allBound[it.getRank()].push_back(*it);	
504

505

506
      // === Check the boundaries and adapt mesh if necessary. ===
507

508
509
510
511
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

512
513
514
515
516
517
518
      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);
519
520
521
522

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

525
526

#if (DEBUG != 0)
527
    debug::writeMesh(feSpace, -1, "mesh");
528
529
530
531
532
533
#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. ===
534

535
    mesh->dofCompress();
536
    updateLocalGlobalNumbering();
537
538
539
540

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

    createPeriodicMap();
541
542
543
  }

  
544
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
545
  {
546
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
547
548
549
550
551
552

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
553
554
555
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
556
	elCode.init(boundIt->rankObj);
557
558
559
560
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
561
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
562
    stdMpi.send(sendCodes);
563
564
565
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
566
    // === Compare received mesh structure codes. ===
567
    
568
569
    bool meshFitTogether = true;

570
571
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
572
573
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
574
      
575
576
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
577

578
579
580
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
581
582
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
583
584

	  bool b = fitElementToMeshCode(recvCodes[i], 
585
					boundIt->rankObj.el,
586
					boundIt->rankObj.subObj,
587
					boundIt->rankObj.ithObj, 
588
589
					boundIt->rankObj.elType,
					boundIt->rankObj.reverseMode);  
590

591
592
	  if (b)
	    meshFitTogether = false;
593
	}
594
	
595
	i++;
596
597
598
      }
    }

599
    return meshFitTogether;
600
  }
601
602


603
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
604
605
606
607
608
					     Element *el, 
					     GeoIndex subObj,
					     int ithObj, 
					     int elType,
					     bool reverseMode)
609
  {
610
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");
611

612
613
    TEST_EXIT_DBG(el)("No element given!\n");

614
615
    if (code.empty())
      return false;
616
617
618

    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
619
620
621

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

622
    bool meshChanged = false;
623
624
625
626
627
628
629
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

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

    if (s1 != -1 && s2 != -1) {
632
      TraverseStack stack;
633
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
634
635
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
636

637
      meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
638
639
      return meshChanged;
    }
640
641

    if (el->isLeaf()) {
642
643
644
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

645
646
647
      meshChanged = true;

      TraverseStack stack;
648
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
649
650
651
652
653
654
655

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

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

656
      el->setMark(1);
657
658
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
659
      refineManager->refineFunction(elInfo);
660
    }
661

662
663
664
665
666
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode)
      std::swap(child0, child1);

667
    if (s1 != -1) {
668
      TraverseStack stack;
669
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
670

671
      while (elInfo && elInfo->getElement() != child0) {
672
673
674
	elInfo = stack.traverseNext(elInfo);
      }

675
      meshChanged |= 
676
	fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
677
    } else {
678
      TraverseStack stack;
679
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
680

681
      while (elInfo && elInfo->getElement() != child1) {
682
683
684
	elInfo = stack.traverseNext(elInfo);
      }

685
      meshChanged |= 
686
	fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
687
688
689
    }

    return meshChanged;
690
691
  }

692

693
  bool MeshDistributor::fitElementToMeshCode2(MeshStructure &code, 
694
695
696
697
698
					      TraverseStack &stack,
					      GeoIndex subObj,
					      int ithObj, 
					      int elType,
					      bool reverseMode)
699
  {
700
    FUNCNAME("MeshDistributor::fitElementToMeshCode2()");
701

702
703
    ElInfo *elInfo = stack.getElInfo();

704
    bool value = false;
705
706
707
708
709
710
711
712
713
714
715
    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);
716
     
717
718
      return value;
    }
719

720
721
    if (!elInfo)
      return value;
722

723
    if (el->isLeaf()) {
724
      el->setMark(1);
725
726
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
727
      refineManager->refineFunction(elInfo);
728
729
730
      value = true;
    }

731
732
    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
733
734
735
736
737
738
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
      std::swap(s1, s2);
      std::swap(child0, child1);
    }
739

740
    if (s1 != -1) {
741
      stack.traverseNext(elInfo);
742
      code.nextElement();
743
      value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
744
745
746
747
      elInfo = stack.getElInfo();
    } else {
      do {
	elInfo = stack.traverseNext(elInfo);
748
      } while (elInfo && elInfo->getElement() != child1);      
749
750
    }  

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

753
    if (s2 != -1) {
754
      code.nextElement();
755
      value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
756
757
    } else {
      int level = elInfo->getLevel();
758

759
760
761
762
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
763

764
765
766
    return value;
  }

767
  
768
  void MeshDistributor::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
769
  {    
770
    int vecSize = data.size();
771
    SerUtil::serialize(out, vecSize);
772
773
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
774
      SerUtil::serialize(out, dofIndex);
775
776
777
778
    }
  }


779
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
780
				    std::map<int, const DegreeOfFreedom*> &dofMap)
781
  {
782
    FUNCNAME("MeshDistributor::deserialize()");
783
784

    int vecSize = 0;
785
    SerUtil::deserialize(in, vecSize);
786
    data.clear();
787
788
789
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
790
      SerUtil::deserialize(in, dofIndex);
791
792
793
794
795
796
797
798
799

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

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


800
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)
801
802
  {
    int mapSize = data.size();
803
    SerUtil::serialize(out, mapSize);
804
805
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
806
      SerUtil::serialize(out, rank);
807
808
809
810
811
      serialize(out, it->second);
    }
  }

  
812
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
813
				    std::map<int, const DegreeOfFreedom*> &dofMap)
814
  {
815
816
    data.clear();

817
    int mapSize = 0;
818
    SerUtil::deserialize(in, mapSize);
819
820
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
821
      SerUtil::deserialize(in, rank);
822
823
824
825
      deserialize(in, data[rank], dofMap);      
    }
  }

826

827
  double MeshDistributor::setElemWeights(AdaptInfo *adaptInfo) 
828
  {
829
    FUNCNAME("MeshDistributor::setElemWeights()");
830

831
    double localWeightSum = 0.0;
832
    elemWeights.clear();
833
834
835
836
837
838
839
840
841
842
843
844
845
846
      
    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;
847

848
849
850
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
851

852
      infile.close();
853
    } else {           
854
      TraverseStack stack;
855
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
856
      while (elInfo) {
857
858
859
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
	localWeightSum++;

860
	elInfo = stack.traverseNext(elInfo);
861
862
863
864
865
866
      }
    }

    return localWeightSum;
  }

867

868
  void MeshDistributor::partitionMesh(AdaptInfo *adaptInfo)
869
  {
870
    FUNCNAME("MeshDistributor::partitionMesh()");
871

872
873
874
875
876
877
878
879
880
881
882
883
    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);
  }

884

885
  void MeshDistributor::createInteriorBoundaryInfo()
886
  {
887
    FUNCNAME("MeshDistributor::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
888

889
    createBoundaryDataStructure();
Thomas Witkowski's avatar
Thomas Witkowski committed
890
891
892
893
894
895
896
897
898

    // === 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);
899

Thomas Witkowski's avatar
Thomas Witkowski committed
900
901
902
903
904
905
906
907
    // === 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 ===
908
      // === indices. Now, we have to sort the corresponding list in this rank to   ===
Thomas Witkowski's avatar
Thomas Witkowski committed
909
      // === get the same order.                                                    ===
910
     
911
      for (unsigned int j = 0; j < rankIt->second.size(); j++) {
912

Thomas Witkowski's avatar
Thomas Witkowski committed
913
914
915
916
917
	// 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) {
918
	  unsigned int k = j + 1;
919

920
	  for (; k < rankIt->second.size(); k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
921
922
923
924
 	    if ((rankIt->second)[k].neighObj == recvedBound)
	      break;

	  // The element must always be found, because the list is just in another order.
925
	  TEST_EXIT_DBG(k < rankIt->second.size())("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
926
927
928
929
930
931
932
933
934

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

935

Thomas Witkowski's avatar
Thomas Witkowski committed
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
    // === 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;
	  
963
	for (unsigned int