MeshDistributor.cc 62.5 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
      if (writePartMesh > 0) {
	debug::writeElementIndexMesh(mesh, "elementIndex.vtu");
110
	writePartitioningMesh("part.vtu");
111
      } else {
112
	MSG("Skip write part mesh!\n");
113
      }
114
    }
115

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

119

120
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
121

122
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
123

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

128

129
130
131
132
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

133

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

    removeMacroElements();
137
138
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
139

140
141
142
143
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

144

145
146
    // === If in debug mode, make some tests. ===

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

154
    debug::writeMesh(feSpace, -1, "macro_mesh");   
155
#endif
156

157

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

Thomas Witkowski's avatar
Thomas Witkowski committed
160
    createPeriodicMap();    
161

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

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

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

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

174
      mesh->dofCompress();
175
      updateLocalGlobalNumbering();
176
177

     
178
      // === Update periodic mapping, if there are periodic boundaries. ===
179
      
180
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
181
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
182

183

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

186
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
187

188

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

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

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

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

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

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

265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
      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;
281
282
283
284
285
286
287
    }

    probStat.push_back(probVec);
  }


  void MeshDistributor::exitParallelization(AdaptInfo *adaptInfo)
288
  {}
289

290
  
291
  void MeshDistributor::updateDofAdmins()
292
  {
293
    FUNCNAME("MeshDistributor::updateDofAdmins()");
294
295

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
296
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
297
298
299

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

308
309
310
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
311
312
313
    }
  }

314

315
  void MeshDistributor::testForMacroMesh()
316
  {
317
    FUNCNAME("MeshDistributor::testForMacroMesh()");
318
319
320
321
322
323
324

    int nMacroElements = 0;

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

      elInfo = stack.traverseNext(elInfo);
    }

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

336

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
361
362
363
    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];
  }
364
365


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

Thomas Witkowski's avatar
Thomas Witkowski committed
383
      stdMpi.send(sendIt->first, dofs);
384
385
386
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    stdMpi.startCommunication<double>(MPI_DOUBLE);
391
392

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

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

406

407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
  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);
      }
    }
  }


426
427
428
429
430
431
432
433
434
435
  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())
436
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
	}
	
	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
458
459
460
461
462
463
464
465
466
467
468
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


469
  void MeshDistributor::checkMeshChange()
470
  {
471
    FUNCNAME("MeshDistributor::checkMeshChange()");
472

473
474
475
476
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "before_check_mesh");
#endif

477
478
479
480
481
    // === 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);
482

483
    if (recvAllValues == 0)
484
485
      return;

486
487
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
488

489
490
491
492
493
    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.
494
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
495
496

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

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

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

508

509
      // === Check the boundaries and adapt mesh if necessary. ===
510

511
512
513
514
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

515
516
517
518
519
520
521
      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);
522
523
524
525

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

528
529

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

538
    mesh->dofCompress();
539
    updateLocalGlobalNumbering();
540
541
542
543

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

    createPeriodicMap();
544
545
546
  }

  
547
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
548
  {
549
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
550
551
552
553
554
555

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
556

557
558
559
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
560
	elCode.init(boundIt->rankObj);
561
562
563
564
	sendCodes[it->first].push_back(elCode);
      }
    }

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

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

582
583
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
584

585
586
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
587
588

	  bool b = fitElementToMeshCode(recvCodes[i], 
589
					boundIt->rankObj.el,
590
					boundIt->rankObj.subObj,
591
					boundIt->rankObj.ithObj, 
592
					boundIt->rankObj.elType,
593
594
					boundIt->neighObj.reverseMode);  

595
596
	  if (b)
	    meshFitTogether = false;
597
 	}
598

599
	i++;
600
601
602
      }
    }

603
    return meshFitTogether;
604
  }
605
606


607
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
608
609
610
611
612
					     Element *el, 
					     GeoIndex subObj,
					     int ithObj, 
					     int elType,
					     bool reverseMode)
613
  {
614
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");
615

616
617
    TEST_EXIT_DBG(el)("No element given!\n");

618
619
    if (code.empty())
      return false;
620
621
622

    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
623
624
625

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

626
    bool meshChanged = false;
627
628
629
630
631
632
633
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

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

    if (s1 != -1 && s2 != -1) {
636
      TraverseStack stack;
637
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
638
639
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
640

641
642
      TEST_EXIT_DBG(elInfo->getElement() == el)("This should not happen!\n");

643
      meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
644
645
      return meshChanged;
    }
646
647

    if (el->isLeaf()) {
648
649
650
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

651
652
653
      meshChanged = true;

      TraverseStack stack;
654
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
655
656
657
658
659
660
661

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

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

662
      el->setMark(1);
663
664
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
665
      refineManager->refineFunction(elInfo);
666
    }
667

668
669
670
671
672
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode)
      std::swap(child0, child1);

673
    if (s1 != -1) {
674
      TraverseStack stack;
675
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
676

677
      while (elInfo && elInfo->getElement() != child0) {
678
679
680
	elInfo = stack.traverseNext(elInfo);
      }

681
      meshChanged |= 
682
	fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
683
    } else {
684
      TraverseStack stack;
685
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
686

687
      while (elInfo && elInfo->getElement() != child1) {
688
689
690
	elInfo = stack.traverseNext(elInfo);
      }

691
      meshChanged |= 
692
	fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
693
694
695
    }

    return meshChanged;
696
697
  }

698

699
  bool MeshDistributor::fitElementToMeshCode2(MeshStructure &code, 
700
701
702
703
704
					      TraverseStack &stack,
					      GeoIndex subObj,
					      int ithObj, 
					      int elType,
					      bool reverseMode)
705
  {
706
    FUNCNAME("MeshDistributor::fitElementToMeshCode2()");
707

708
    ElInfo *elInfo = stack.getElInfo();
709
    bool value = false;
710
711
712
713
714
715
716
717
718
719
720
    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);
721

722
723
      return value;
    }
724

725
726
    if (!elInfo)
      return value;
727

728
    if (el->isLeaf()) {
729
      el->setMark(1);
730
731
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
732
      refineManager->refineFunction(elInfo);
733
734
735
      value = true;
    }

736
737
    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
738
739
740
741
742
743
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
      std::swap(s1, s2);
      std::swap(child0, child1);
    }
744

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

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

758
    if (s2 != -1) {
759
      code.nextElement();
760
      value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
761
762
    } else {
      int level = elInfo->getLevel();
763

764
765
766
767
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
768

769
770
771
    return value;
  }

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


784
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
785
				    std::map<int, const DegreeOfFreedom*> &dofMap)
786
  {
787
    FUNCNAME("MeshDistributor::deserialize()");
788
789

    int vecSize = 0;
790
    SerUtil::deserialize(in, vecSize);
791
    data.clear();
792
793
794
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
795
      SerUtil::deserialize(in, dofIndex);
796
797
798
799
800
801
802
803
804

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

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


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

  
817
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
818
				    std::map<int, const DegreeOfFreedom*> &dofMap)
819
  {
820
821
    data.clear();

822
    int mapSize = 0;
823
    SerUtil::deserialize(in, mapSize);
824
825
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
826
      SerUtil::deserialize(in, rank);
827
828
829
830
      deserialize(in, data[rank], dofMap);      
    }
  }

831

832
  double MeshDistributor::setElemWeights(AdaptInfo *adaptInfo) 
833
  {
834
    FUNCNAME("MeshDistributor::setElemWeights()");
835

836
    double localWeightSum = 0.0;
837
    elemWeights.clear();
838
839
840
841
842
843
844
845
846
847
848
849
850
851
      
    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;
852

853
854
855
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
856

857
      infile.close();
858
    } else {           
859
      TraverseStack stack;