MeshDistributor.cc 72.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
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
241
242
243
244
245
246
247
248
249
250
251
252
253
254
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
      std::string f = "";
      GET_PARAMETER(0, probVec->getName() + "->output->serialization filename", &f);
      path myPath(f);
      std::string meshFilename = 
	myPath.parent_path().directory_string() + "/meshDistributor.ser";


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

    int readSerialization = 0;
257
258
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", 
		  &readSerialization);
259
260
261
262
    if (readSerialization) {
      std::string filename = "";
      GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
      filename += ".p" + lexical_cast<std::string>(mpiRank);
263
      MSG("Start deserialization with %s\n", filename.c_str());
264
      std::ifstream in(filename.c_str());
265
      probVec->deserialize(in);
266
      in.close();
267
268
269
270
271
272
273
274
275
276
277
278
279
280
      MSG("Deserialization from file: %s\n", filename.c_str());

      if (!deserialized) {
	GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
	path myPath(filename);
	std::string meshFilename = 
	  myPath.parent_path().directory_string() + "/meshDistributor.ser.p" + 
	  lexical_cast<std::string>(mpiRank);
	std::ifstream in(meshFilename.c_str());
	deserialize(in);
	in.close();
	MSG("Deserializtion of mesh distributor from file: %s\n", meshFilename.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
426
427
428
429
430
431
432
433
434
435
436
437
438
  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
439
440
441
442
443
444
445
446
447
448
449
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


450
  void MeshDistributor::checkMeshChange()
451
  {
452
    FUNCNAME("MeshDistributor::checkMeshChange()");
453

454
455
456
457
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "before_check_mesh");
#endif

458
459
460
461
462
    // === 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);
463

464
    if (recvAllValues == 0)
465
466
      return;

467
468
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
469

470
471
472
473
474
    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.
475
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
476
477

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

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

485
486
487
488
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
	  allBound[it.getRank()].push_back(*it);	

489

490
      // === Check the boundaries and adapt mesh if necessary. ===
491

492
493
494
495
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

496
497
498
499
500
501
502
      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);
503
504
505
506

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

509
510

#if (DEBUG != 0)
511
    debug::writeMesh(feSpace, -1, "mesh");
512
513
514
515
516
517
#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. ===
518

519
    mesh->dofCompress();
520
    updateLocalGlobalNumbering();
521
522
523
524

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

    createPeriodicMap();
525
526
527
  }

  
528
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
529
  {
530
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
531
532
533
534
535
536

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
537
538
539
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
540
	elCode.init(boundIt->rankObj);
541
542
543
544
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
545
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
546
    stdMpi.send(sendCodes);
547
548
549
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
550
    // === Compare received mesh structure codes. ===
551
    
552
553
    bool meshFitTogether = true;

554
555
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
556
557
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
558
      
559
560
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
561

562
563
564
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
565
566
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
567
568

	  bool b = fitElementToMeshCode(recvCodes[i], 
569
					boundIt->rankObj.el,
570
					boundIt->rankObj.subObj,
571
					boundIt->rankObj.ithObj, 
572
573
					boundIt->rankObj.elType,
					boundIt->rankObj.reverseMode);  
574

575
576
	  if (b)
	    meshFitTogether = false;
577
	}
578
	
579
	i++;
580
581
582
      }
    }

583
    return meshFitTogether;
584
  }
585
586


587
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
588
589
590
591
592
					     Element *el, 
					     GeoIndex subObj,
					     int ithObj, 
					     int elType,
					     bool reverseMode)
593
  {
594
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");
595

596
597
    TEST_EXIT_DBG(el)("No element given!\n");

598
599
    if (code.empty())
      return false;
600
601
602

    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
603
604
605

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

606
    bool meshChanged = false;
607
608
609
610
611
612
613
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

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

    if (s1 != -1 && s2 != -1) {
616
      TraverseStack stack;
617
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
618
619
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
620

621
      meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
622
623
      return meshChanged;
    }
624
625

    if (el->isLeaf()) {
626
627
628
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

629
630
631
      meshChanged = true;

      TraverseStack stack;
632
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
633
634
635
636
637
638
639

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

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

640
      el->setMark(1);
641
642
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
643
      refineManager->refineFunction(elInfo);
644
    }
645

646
647
648
649
650
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode)
      std::swap(child0, child1);

651
    if (s1 != -1) {
652
      TraverseStack stack;
653
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
654

655
      while (elInfo && elInfo->getElement() != child0) {
656
657
658
	elInfo = stack.traverseNext(elInfo);
      }

659
      meshChanged |= 
660
	fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
661
    } else {
662
      TraverseStack stack;
663
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
664

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

669
      meshChanged |= 
670
	fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
671
672
673
    }

    return meshChanged;
674
675
  }

676

677
  bool MeshDistributor::fitElementToMeshCode2(MeshStructure &code, 
678
679
680
681
682
					      TraverseStack &stack,
					      GeoIndex subObj,
					      int ithObj, 
					      int elType,
					      bool reverseMode)
683
  {
684
    FUNCNAME("MeshDistributor::fitElementToMeshCode2()");
685

686
687
    ElInfo *elInfo = stack.getElInfo();

688
    bool value = false;
689
690
691
692
693
694
695
696
697
698
699
    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);
700
     
701
702
      return value;
    }
703

704
705
    if (!elInfo)
      return value;
706

707
    if (el->isLeaf()) {
708
      el->setMark(1);
709
710
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
711
      refineManager->refineFunction(elInfo);
712
713
714
      value = true;
    }

715
716
    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
717
718
719
720
721
722
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
      std::swap(s1, s2);
      std::swap(child0, child1);
    }
723

724
    if (s1 != -1) {
725
      stack.traverseNext(elInfo);
726
      code.nextElement();
727
      value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
728
729
730
731
      elInfo = stack.getElInfo();
    } else {
      do {
	elInfo = stack.traverseNext(elInfo);
732
      } while (elInfo && elInfo->getElement() != child1);      
733
734
    }  

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

737
    if (s2 != -1) {
738
      code.nextElement();
739
      value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
740
741
    } else {
      int level = elInfo->getLevel();
742

743
744
745
746
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
747

748
749
750
    return value;
  }

751
  
752
  void MeshDistributor::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
753
  {    
754
    int vecSize = data.size();
755
    SerUtil::serialize(out, vecSize);
756
757
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
758
      SerUtil::serialize(out, dofIndex);
759
760
761
762
    }
  }


763
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
764
				    std::map<int, const DegreeOfFreedom*> &dofMap)
765
  {
766
    FUNCNAME("MeshDistributor::deserialize()");
767
768

    int vecSize = 0;
769
    SerUtil::deserialize(in, vecSize);
770
771
772
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
773
      SerUtil::deserialize(in, dofIndex);
774
775
776
777
778
779
780
781
782

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

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


783
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)
784
785
  {
    int mapSize = data.size();
786
    SerUtil::serialize(out, mapSize);
787
788
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
789
      SerUtil::serialize(out, rank);
790
791
792
793
794
      serialize(out, it->second);
    }
  }

  
795
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
796
				    std::map<int, const DegreeOfFreedom*> &dofMap)
797
798
  {
    int mapSize = 0;
799
    SerUtil::deserialize(in, mapSize);
800
801
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
802
      SerUtil::deserialize(in, rank);
803
804
805
806
      deserialize(in, data[rank], dofMap);      
    }
  }

807

808
  double MeshDistributor::setElemWeights(AdaptInfo *adaptInfo) 
809
  {
810
    FUNCNAME("MeshDistributor::setElemWeights()");
811

812
    double localWeightSum = 0.0;
813
    elemWeights.clear();
814
815
816
817
818
819
820
821
822
823
824
825
826
827
      
    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;
828

829
830
831
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
832

833
      infile.close();
834
    } else {           
835
      TraverseStack stack;
836
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
837
      while (elInfo) {
838
839
840
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
	localWeightSum++;

841
	elInfo = stack.traverseNext(elInfo);
842
843
844
845
846
847
      }
    }

    return localWeightSum;
  }

848

849
  void MeshDistributor::partitionMesh(AdaptInfo *adaptInfo)
850
  {
851
    FUNCNAME("MeshDistributor::partitionMesh()");
852

853
854
855
856
857
858
859
860
861
862
863
864
    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);
  }

865

866
  void MeshDistributor::createInteriorBoundaryInfo()
867
  {
868
    FUNCNAME("MeshDistributor::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
869

Thomas Witkowski's avatar
Thomas Witkowski committed
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
    // === First, create the interior boundaries based on macro element's  ===
    // === neighbour informations.                                         ===

    createMacroElementInteriorBoundaryInfo();

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

    createSubstructureInteriorBoundaryInfo();


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

    StdMpi<std::vector<AtomicBoundary> > stdMpi(mpiComm);
    stdMpi.send(myIntBoundary.boundary);
    stdMpi.recv(otherIntBoundary.boundary);
    stdMpi.startCommunication<int>(MPI_INT);
889

Thomas Witkowski's avatar
Thomas Witkowski committed
890
891
892
893
894
895
896
897
    // === 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 ===
898