MeshDistributor.cc 72.7 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
79
    // If the problem has been already read from a file, we need only to remove
    // periodic boundary conditions (if there are some).
    if (deserialized) {
      removePeriodicBoundaryConditions();
80
      return;
81
    }
82
    
83

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

89

90
91
92
93
94
95
96
    // 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);   

97

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

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

111
    ParallelDebug::testAllElements(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
112
#endif
113

114

115
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
116

117
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
118

119
#if (DEBUG != 0)
120
    ParallelDebug::printBoundaryInfo(*this);
121
122
#endif

123

124
125
126
127
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

128

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
129
130
131
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
132
133
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
134

135
136
137
138
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

139

140
141
    // === If in debug mode, make some tests. ===

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

149
    debug::writeMesh(feSpace, -1, "macro_mesh");   
150
#endif
151

152

153
154
155
    // === Create periodic dof mapping, if there are periodic boundaries. ===

    createPeriodicMap();
156

157

158
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
159

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

Thomas Witkowski's avatar
Thomas Witkowski committed
163
    if (globalRefinement > 0) {
164
      refineManager->globalRefine(mesh, globalRefinement);
165

166
#if (DEBUG != 0)
167
      debug::writeMesh(feSpace, -1, "gr_mesh");
168
169
#endif

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

177

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

180
181
182
183
184
185
    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
186

187
188
189
190
191
192
	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
193
194
    }

195

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

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

      int tsModulo = -1;
      GET_PARAMETER(0, probVec->getName() + "->output->write every i-th timestep", 
		    "%d", &tsModulo);
      
251
      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
252
253
      writeSerializationFile = true;
    }    
254
255

    int readSerialization = 0;
256
257
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", 
		  &readSerialization);
258
259
260
261
    if (readSerialization) {
      std::string filename = "";
      GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
      filename += ".p" + lexical_cast<std::string>(mpiRank);
262
      MSG("Start deserialization with %s\n", filename.c_str());
263
      std::ifstream in(filename.c_str());
264
265
266
267

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

268
      probVec->deserialize(in);
269
      in.close();
270
271
272
      MSG("Deserialization from file: %s\n", filename.c_str());

      if (!deserialized) {
273
274
275
276
277
278
279
280
281
282
283
284
	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);
	std::ifstream in(rankFilename.c_str());

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

285
286
	deserialize(in);
	in.close();
287
	MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
288
289
	deserialized = true;
      }
290
291
292
293
294
295
296
    }

    probStat.push_back(probVec);
  }


  void MeshDistributor::exitParallelization(AdaptInfo *adaptInfo)
297
  {}
298

299
  
300
  void MeshDistributor::updateDofAdmins()
301
  {
302
    FUNCNAME("MeshDistributor::updateDofAdmins()");
303
304

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
305
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
306
307
308

      // There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise,
      // it is not possible to define a first DOF hole.
309
      if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDofs.size())
310
	admin.enlargeDOFLists();
311
312
313
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
314
      for (unsigned int j = 0; j < mapLocalGlobalDofs.size(); j++)
315
316
 	admin.setDOFFree(j, false);

317
318
319
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
320
321
322
    }
  }

323

324
  void MeshDistributor::testForMacroMesh()
325
  {
326
    FUNCNAME("MeshDistributor::testForMacroMesh()");
327
328
329
330
331
332
333

    int nMacroElements = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      TEST_EXIT(elInfo->getLevel() == 0)
334
	("Mesh is already refined! This does not work with parallelization!\n");
335
336
337
338
339
340
341
342
343
344
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

345

346
  void MeshDistributor::synchVector(DOFVector<double> &vec)
347
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
348
    StdMpi<std::vector<double> > stdMpi(mpiComm);
349
350

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
351
	 sendIt != sendDofs.end(); ++sendIt) {
352
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
353
354
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
355
      
Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
358
359
360
361

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

Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
364
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
365

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

Thomas Witkowski's avatar
Thomas Witkowski committed
368
369
370
371
372
    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];
  }
373
374


375
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
376
  {
377
    int nComponents = vec.getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
378
379
380
381
382
383
384
385
386
387
388
389
    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])]);
390
391
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
392
      stdMpi.send(sendIt->first, dofs);
393
394
395
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
399
    stdMpi.startCommunication<double>(MPI_DOUBLE);
400
401

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
402
403
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
404
405

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
406
407
408
409
410
      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++];
411
412
413
414
      }
    }
  }

415

416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
  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())
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));	  
	}
	
	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
448
449
450
451
452
453
454
455
456
457
458
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


459
  void MeshDistributor::checkMeshChange()
460
  {
461
    FUNCNAME("MeshDistributor::checkMeshChange()");
462

463
464
465
466
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "before_check_mesh");
#endif

467
468
469
470
471
    // === 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);
472

473
    if (recvAllValues == 0)
474
475
      return;

476
477
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
478

479
480
481
482
483
    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.
484
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486

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

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

494
495
496
497
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
	  allBound[it.getRank()].push_back(*it);	

498

499
      // === Check the boundaries and adapt mesh if necessary. ===
500

501
502
503
504
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

505
506
507
508
509
510
511
      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);
512
513
514
515

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

518
519

#if (DEBUG != 0)
520
    debug::writeMesh(feSpace, -1, "mesh");
521
522
523
524
525
526
#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. ===
527

528
    mesh->dofCompress();
529
    updateLocalGlobalNumbering();
530
531
532
533

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

    createPeriodicMap();
534
535
536
  }

  
537
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
538
  {
539
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
540
541
542
543
544
545

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
546
547
548
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
549
	elCode.init(boundIt->rankObj);
550
551
552
553
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
554
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
555
    stdMpi.send(sendCodes);
556
557
558
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
559
    // === Compare received mesh structure codes. ===
560
    
561
562
    bool meshFitTogether = true;

563
564
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
565
566
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
567
      
568
569
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
570

571
572
573
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
574
575
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
576
577

	  bool b = fitElementToMeshCode(recvCodes[i], 
578
					boundIt->rankObj.el,
579
					boundIt->rankObj.subObj,
580
					boundIt->rankObj.ithObj, 
581
582
					boundIt->rankObj.elType,
					boundIt->rankObj.reverseMode);  
583

584
585
	  if (b)
	    meshFitTogether = false;
586
	}
587
	
588
	i++;
589
590
591
      }
    }

592
    return meshFitTogether;
593
  }
594
595


596
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
597
598
599
600
601
					     Element *el, 
					     GeoIndex subObj,
					     int ithObj, 
					     int elType,
					     bool reverseMode)
602
  {
603
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");
604

605
606
    TEST_EXIT_DBG(el)("No element given!\n");

607
608
    if (code.empty())
      return false;
609
610
611

    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
612
613
614

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

615
    bool meshChanged = false;
616
617
618
619
620
621
622
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

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

    if (s1 != -1 && s2 != -1) {
625
      TraverseStack stack;
626
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
627
628
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
629

630
      meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
631
632
      return meshChanged;
    }
633
634

    if (el->isLeaf()) {
635
636
637
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

638
639
640
      meshChanged = true;

      TraverseStack stack;
641
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
642
643
644
645
646
647
648

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

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

649
      el->setMark(1);
650
651
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
652
      refineManager->refineFunction(elInfo);
653
    }
654

655
656
657
658
659
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode)
      std::swap(child0, child1);

660
    if (s1 != -1) {
661
      TraverseStack stack;
662
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
663

664
      while (elInfo && elInfo->getElement() != child0) {
665
666
667
	elInfo = stack.traverseNext(elInfo);
      }

668
      meshChanged |= 
669
	fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
670
    } else {
671
      TraverseStack stack;
672
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
673

674
      while (elInfo && elInfo->getElement() != child1) {
675
676
677
	elInfo = stack.traverseNext(elInfo);
      }

678
      meshChanged |= 
679
	fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
680
681
682
    }

    return meshChanged;
683
684
  }

685

686
  bool MeshDistributor::fitElementToMeshCode2(MeshStructure &code, 
687
688
689
690
691
					      TraverseStack &stack,
					      GeoIndex subObj,
					      int ithObj, 
					      int elType,
					      bool reverseMode)
692
  {
693
    FUNCNAME("MeshDistributor::fitElementToMeshCode2()");
694

695
696
    ElInfo *elInfo = stack.getElInfo();

697
    bool value = false;
698
699
700
701
702
703
704
705
706
707
708
    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);
709
     
710
711
      return value;
    }
712

713
714
    if (!elInfo)
      return value;
715

716
    if (el->isLeaf()) {
717
      el->setMark(1);
718
719
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
720
      refineManager->refineFunction(elInfo);
721
722
723
      value = true;
    }

724
725
    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
726
727
728
729
730
731
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
      std::swap(s1, s2);
      std::swap(child0, child1);
    }
732

733
    if (s1 != -1) {
734
      stack.traverseNext(elInfo);
735
      code.nextElement();
736
      value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
737
738
739
740
      elInfo = stack.getElInfo();
    } else {
      do {
	elInfo = stack.traverseNext(elInfo);
741
      } while (elInfo && elInfo->getElement() != child1);      
742
743
    }  

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

746
    if (s2 != -1) {
747
      code.nextElement();
748
      value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
749
750
    } else {
      int level = elInfo->getLevel();
751

752
753
754
755
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
756

757
758
759
    return value;
  }

760
  
761
  void MeshDistributor::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
762
  {    
763
    int vecSize = data.size();
764
    SerUtil::serialize(out, vecSize);
765
766
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
767
      SerUtil::serialize(out, dofIndex);
768
769
770
771
    }
  }


772
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
773
				    std::map<int, const DegreeOfFreedom*> &dofMap)
774
  {
775
    FUNCNAME("MeshDistributor::deserialize()");
776
777

    int vecSize = 0;
778
    SerUtil::deserialize(in, vecSize);
779
780
781
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
782
      SerUtil::deserialize(in, dofIndex);
783
784
785
786
787
788
789
790
791

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

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


792
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)
793
794
  {
    int mapSize = data.size();
795
    SerUtil::serialize(out, mapSize);
796
797
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
798
      SerUtil::serialize(out, rank);
799
800
801
802
803
      serialize(out, it->second);
    }
  }

  
804
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
805
				    std::map<int, const DegreeOfFreedom*> &dofMap)
806
807
  {
    int mapSize = 0;
808
    SerUtil::deserialize(in, mapSize);
809
810
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
811
      SerUtil::deserialize(in, rank);
812
813
814
815
      deserialize(in, data[rank], dofMap);      
    }
  }

816

817
  double MeshDistributor::setElemWeights(AdaptInfo *adaptInfo) 
818
  {
819
    FUNCNAME("MeshDistributor::setElemWeights()");
820

821
    double localWeightSum = 0.0;
822
    elemWeights.clear();
823
824
825
826
827
828
829
830
831
832
833
834
835
836
      
    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;
837

838
839
840
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
841

842
      infile.close();
843
    } else {           
844
      TraverseStack stack;
845
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
846
      while (elInfo) {
847
848
849
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
	localWeightSum++;

850
	elInfo = stack.traverseNext(elInfo);
851
852
853
854
855
856
      }
    }

    return localWeightSum;
  }

857

858
  void MeshDistributor::partitionMesh(AdaptInfo *adaptInfo)
859
  {
860
    FUNCNAME("MeshDistributor::partitionMesh()");
861

862
863
864
865
866
867
868
869
870
871
872
873
    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);
  }

874

875
  void MeshDistributor::createInteriorBoundaryInfo()
876
  {