MeshDistributor.cc 63.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
110
111
112
      if (writePartMesh > 0)
	writePartitioningMesh("part.vtu");
      else 
	MSG("Skip write part mesh!\n");
    }
113

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

117

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

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

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

126

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

    createLocalGlobalNumbering();

131

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

    removeMacroElements();
135
136
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
137

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

    updateDofAdmins();

142

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

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

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

155

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

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

160

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

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

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

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

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

181

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

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

186

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

189
    removePeriodicBoundaryConditions();
190
191
  }

192

193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
  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;
229
230
231
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
232
233
234
235
236
      std::string filename = "";
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
237
238
239
240
241

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

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

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

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

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

    probStat.push_back(probVec);
  }


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

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

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

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

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

312

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

    int nMacroElements = 0;

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

      elInfo = stack.traverseNext(elInfo);
    }

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

334

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

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

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

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

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

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


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

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

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

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

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

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

404

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


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


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

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

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

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

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

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

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

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

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

506

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

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

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

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

526
527

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

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

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

    createPeriodicMap();
542
543
544
  }

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

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

      //      MSG("HAVE %d BOUNDARIES WITH RANK %d!\n", it->second.size(), it->first);

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
576
577
578
579
580
581
582
//     if (mpiRank == 0) {
//       debug::highlightElementIndexMesh(mesh, 139, "rank0_139.vtu");
//       debug::highlightElementIndexMesh(mesh, 161, "rank0_161.vtu");
//       debug::printElementHierarchie(mesh, 139);
//     } else {
//       debug::highlightElementIndexMesh(mesh, 137, "rank1_137.vtu");
//       debug::highlightElementIndexMesh(mesh, 163, "rank1_163.vtu");
//     }

583
584
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
585
586
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
587
      
588
589
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
590

591
	MeshStructure elCode;	
592
593
594
// 	if (i == 13) 
// 	  elCode.setDebugMode(true);

595
596
	elCode.init(boundIt->rankObj);
	
597
598
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
599

600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
	  /*
	  MSG("----------------------------------------------------------------\n");
 	  MSG("CODE %d: MESH SUB OBJ = %d   INDEX = %d   ITH = %d   REVERSE = %d\n", 
	      i,
 	      boundIt->rankObj.subObj,
 	      boundIt->rankObj.el->getIndex(),
	      boundIt->rankObj.ithObj,
	      boundIt->rankObj.reverseMode);

	  MSG("RANK CODE: \n");
	  elCode.print();
	  elCode.reset();

	  MSG("OTHER CODE: \n");
	  recvCodes[i].print();
	  recvCodes[i].reset();
	  */
617
	  bool b = fitElementToMeshCode(recvCodes[i], 
618
					boundIt->rankObj.el,
619
					boundIt->rankObj.subObj,
620
					boundIt->rankObj.ithObj, 
621
					boundIt->rankObj.elType,
622
623
624
					boundIt->neighObj.reverseMode);  

	  //	  MSG("Mesh changed in bound from %d to %d: %d\n", mpiRank, it->first, b);
625

626
627
	  if (b)
	    meshFitTogether = false;
628
629
630
631
632
633
634
635
636
637
638
639
640
641
 	} else {
	  /*
	  MSG("----------------------------------------------------------------\n");
 	  MSG("CODE %d is equal!\n", i);

	  MSG("RANK CODE: \n");
	  elCode.print();
	  elCode.reset();

	  MSG("OTHER CODE: \n");
	  recvCodes[i].print();
	  recvCodes[i].reset();
	  */
 	}
642
	
643
	i++;
644
645
646
      }
    }

647
    return meshFitTogether;
648
  }
649
650


651
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
652
653
654
655
656
					     Element *el, 
					     GeoIndex subObj,
					     int ithObj, 
					     int elType,
					     bool reverseMode)
657
  {
658
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");
659

660
661
    TEST_EXIT_DBG(el)("No element given!\n");

662
663
    if (code.empty())
      return false;
664

665
666
667
668
669
670
671
//     MSG("FIT EL %d!\n", el->getIndex());
//     if (mpiRank == 0) {
//       debug::highlightElementIndexMesh(mesh, 139, "test.vtu");
//     }

//    MSG("REVERSE MODE = %d\n", reverseMode);

672
673
    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
674
675
676

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

677
    bool meshChanged = false;
678
679
680
681
682
683
684
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

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

686
687
    //    MSG("SUB-OBJ (SWAPED): %d %d\n", s1, s2);

688
    if (s1 != -1 && s2 != -1) {
689
      //      MSG("GO TO FIT 1\n");
690
      TraverseStack stack;
691
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
692
693
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
694

695
696
697
698
      TEST_EXIT_DBG(elInfo->getElement() == el)("This should not happen!\n");

      //      MSG("START WITH EL = %d\n", elInfo->getElement()->getIndex());

699
      meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
700
701
      return meshChanged;
    }
702
703

    if (el->isLeaf()) {
704
705
706
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

707
708
      //      MSG("GO TO FIT 2\n");

709
710
711
      meshChanged = true;

      TraverseStack stack;
712
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
713
714
715
716
717
718
719

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

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

720
      el->setMark(1);
721
722
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
723
      refineManager->refineFunction(elInfo);
724
    }
725

726
727
728
729
730
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode)
      std::swap(child0, child1);

731
    if (s1 != -1) {
732
733
      //      MSG("GO TO FIT 3\n");

734
      TraverseStack stack;
735
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
736

737
      while (elInfo && elInfo->getElement() != child0) {
738
739
740
	elInfo = stack.traverseNext(elInfo);
      }

741
      meshChanged |= 
742
	fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
743
    } else {
744
745
      //      MSG("GO TO FIT 4\n");

746
      TraverseStack stack;
747
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
748

749
      while (elInfo && elInfo->getElement() != child1) {
750
751
752
	elInfo = stack.traverseNext(elInfo);
      }

753
      meshChanged |= 
754
	fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
755
756
757
    }

    return meshChanged;
758
759
  }

760

761
  bool MeshDistributor::fitElementToMeshCode2(MeshStructure &code, 
762
763
764
765
766
					      TraverseStack &stack,
					      GeoIndex subObj,
					      int ithObj, 
					      int elType,
					      bool reverseMode)
767
  {
768
    FUNCNAME("MeshDistributor::fitElementToMeshCode2()");
769

770
771
    ElInfo *elInfo = stack.getElInfo();

772
773
774
775
776
    //    MSG("IN %d\n", elInfo->getElement()->getIndex());

//     MeshStructure c = code;
//     c.print(false);
    
777
    bool value = false;
778
779
780
781
782
783
    if (!elInfo)
      return value;

    Element *el = elInfo->getElement();

    if (code.isLeafElement()) {
784
      //      MSG("HERE 1\n");
785
786
787
788
789
      int level = elInfo->getLevel();

      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
790
791
792
793
794
795

//       if (elInfo) {
// 	MSG("NEXT ELEMENT IS: %d\n", elInfo->getElement()->getIndex());
//       } else {
// 	MSG("NO NEXT ELEMENT!\n");
//       }
796
     
797
798
      return value;
    }
799

800
801
    if (!elInfo)
      return value;
802

803
    if (el->isLeaf()) {
804
805
806
807
//       MSG("HERE 2: %d\n", el->getIndex());
//       if (mpiRank == 0) {
// 	debug::highlightElementIndexMesh(mesh, 1592, "test_1592.vtu");
//       }
808
      el->setMark(1);
809
810
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
811
      refineManager->refineFunction(elInfo);
812
813
814
      value = true;
    }

815
816
    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
817
818
819
820
821
822
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
      std::swap(s1, s2);
      std::swap(child0, child1);
    }
823

824
    if (s1 != -1) {
825
      //      MSG("HERE 3\n");
826
      stack.traverseNext(elInfo);
827
      code.nextElement();
828
      value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
829
830
      elInfo = stack.getElInfo();
    } else {
831
      //      MSG("HERE 4\n");
832
833
      do {
	elInfo = stack.traverseNext(elInfo);
834
835
      } while (elInfo && elInfo->getElement() != child1);  
      //      MSG("S1 == -1   -> SWITCHED TO EL: %d\n", elInfo->getElement()->getIndex());
836
837
    }  

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

840
    if (s2 != -1) {
841
      //      MSG("HERE 5\n");
842
      code.nextElement();
843
      value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
844
    } else {
845
      //      MSG("HERE 6\n");
846
      int level = elInfo->getLevel();
847

848
849
850
851
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
852

853
854
855
    return value;
  }

856
  
857
  void MeshDistributor::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
858
  {    
859
    int vecSize = data.size();
860
    SerUtil::serialize(out, vecSize);
861
862
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
863
      SerUtil::serialize(out, dofIndex);
864
865
866
867
    }
  }


868
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
869
				    std::map<int, const DegreeOfFreedom*> &dofMap)
870
  {
871
    FUNCNAME("MeshDistributor::deserialize()");
872
873

    int vecSize = 0;
874
    SerUtil::deserialize(in, vecSize);
875
    data.clear();
876
877
878
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
879
      SerUtil::deserialize(in, dofIndex);
880
881
882
883
884
885
886
887
888

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

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


889
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)
890
891
  {
    int mapSize = data.size();
892
    SerUtil::serialize(out, mapSize);
893
894
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
895
      SerUtil::serialize(out, rank);
896
897
898
899
900
      serialize(out, it->second);
    }
  }

  
901
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
Thomas Witkowski's avatar