MeshDistributor.cc 65.2 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(feSpace)("No FE space has been defined for the mesh distributor!\n");
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
76

77
78
79
    // 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).
80
    if (deserialized) {
81
      setRankDofs();
82
      removePeriodicBoundaryConditions();
83

84
      return;
85
    }
86

87
   
88
89
90
91
92
    // 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();

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

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

107
108
      if (writePartMesh > 0) {
	debug::writeElementIndexMesh(mesh, "elementIndex.vtu");
109
	writePartitioningMesh("part.vtu");
110
      } else {
111
	MSG("Skip write part mesh!\n");
112
      }
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
    ParallelDebug::testGlobalIndexByCoords(*this);
151
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
152

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

156

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

Thomas Witkowski's avatar
Thomas Witkowski committed
159
    createPeriodicMap();    
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

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

182

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

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

187

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

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

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

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

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

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

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

    probStat.push_back(probVec);
  }


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

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

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

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

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

313

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

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
328
329
330
331
332
333
334
335
336
337
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

338

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

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

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

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

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

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


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

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

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

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

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

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

408

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


428
429
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
430
431
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

432
433
434
435
436
437
438
439
    // 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())
440
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
	}
	
	if (probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
	
	if (probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
      }
    }

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh,  -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
458
459
460
461
462
463

    // Remove periodic vertex associations
    for (std::map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	 it != mesh->getPeriodicAssociations().end(); ++it)
      const_cast<DOFAdmin&>(mesh->getDOFAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
    mesh->getPeriodicAssociations().clear();
464
465
466
467
  }


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
468
469
470
471
472
473
474
475
476
477
478
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


479
  void MeshDistributor::checkMeshChange()
480
  {
481
    FUNCNAME("MeshDistributor::checkMeshChange()");    
482

483
484
485
486
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "before_check_mesh");
#endif

487
488
489
490
491
    // === 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);
492

493
    if (recvAllValues == 0)
494
495
      return;

496
497
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
498

499
500
501
502
503
    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.
504
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
505
506

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

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

514
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
515
	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
516
 	  allBound[it.getRank()].push_back(*it);	
517

518

519
      // === Check the boundaries and adapt mesh if necessary. ===
520
      MSG("TEST 1: %d\n", mesh->getDOFAdmin(0).getUsedSize());
521
522
523
524
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

525
526
527
528
529
530
531
      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);
532
533
534
535

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

538
    MSG("TEST 2: %d\n", mesh->getDOFAdmin(0).getUsedSize());
539
540

#if (DEBUG != 0)
541
    debug::writeMesh(feSpace, -1, "mesh");
542
543
544
545
546
547
#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. ===
548

549
    mesh->dofCompress();
550
    updateLocalGlobalNumbering();
551
552
553
554

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

    createPeriodicMap();
555
556
557
  }

  
558
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
559
  {
560
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
561
562
563
564
565
566

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

568
569
570
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
571
	elCode.init(boundIt->rankObj);
572
573
574
575
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
576
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
577
    stdMpi.send(sendCodes);
578
579
580
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
581
    // === Compare received mesh structure codes. ===
582
    
583
584
    bool meshFitTogether = true;

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

593
594
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
595

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

599
600
601
602
603
604
	  bool b = startFitElementToMeshCode(recvCodes[i], 
					     boundIt->rankObj.el,
					     boundIt->rankObj.subObj,
					     boundIt->rankObj.ithObj, 
					     boundIt->rankObj.elType,
					     boundIt->rankObj.reverseMode);
605

606
607
	  if (b)
	    meshFitTogether = false;
608
 	}
609

610
	i++;
611
612
613
      }
    }

614
    return meshFitTogether;
615
  }
616
617


618
619
620
621
622
623
  bool MeshDistributor::startFitElementToMeshCode(MeshStructure &code, 
						  Element *el, 
						  GeoIndex subObj,
						  int ithObj, 
						  int elType,
						  bool reverseMode)
624
  {
625
    FUNCNAME("MeshDistributor::startFitElementToMeshCode()");
626

627
628
    TEST_EXIT_DBG(el)("No element given!\n");

629
630
    // If the code is empty, the element does not matter and the function can
    // return without chaning the mesh.
631
632
    if (code.empty())
      return false;
633

634
635
636
637
638
    // s0 and s1 are the number of the edge/face in both child of the element,
    // which contain the edge/face the function has to traverse through. If the
    // edge/face is not contained in one of the children, s0 or s1 is -1.
    int s0 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s1 = el->getSubObjOfChild(1, subObj, ithObj, elType);
639

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

642
    bool meshChanged = false;
643
644
645
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

646
647
648
649
650
    // Test for reverse mode, in which the left and right children of elements
    // are flipped.
    if (reverseMode)
      traverseFlag |= Mesh::CALL_REVERSE_MODE;    

651

652
653
654
655
656
    // === If the edge/face is contained in both children. ===

    if (s0 != -1 && s1 != -1) {
      // Create traverse stack and traverse within the mesh until the element,
      // which should be fitted to the mesh structure code, is reached.
657
      TraverseStack stack;
658
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
659
660
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
661

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

664
      meshChanged = fitElementToMeshCode(code, stack, subObj, ithObj, reverseMode);
665
666
      return meshChanged;
    }
667

668
669
670

    // === The edge/face is contained in only on of the both children. ===

671
    if (el->isLeaf()) {
672
673

      // If element is leaf and code contains only one leaf element, we are finished.
674
675
676
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

677
      // Create traverse stack and traverse the mesh to the element.
678
      TraverseStack stack;
679
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
680
681
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
682
683
684

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

685
      // Code is not leaf, therefore refine the element.
686
      el->setMark(1);
687
688
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
689
      refineManager->refineFunction(elInfo);
690
      meshChanged = true;
691
    }
692

693
694
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
695
696
697
698
    if (reverseMode) {
      std::swap(s0, s1);
      std::swap(child0, child1);    
    }
699

700
701
702
    // === We know that the edge/face is contained in only one of the children. ===
    // === Therefore, traverse the mesh to this children and fit this element   ===
    // === To the mesh structure code.                                          ===
703

704
705
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
706

707
708
709
    if (s0 != -1) {
      while (elInfo && elInfo->getElement() != child0)
	elInfo = stack.traverseNext(elInfo);     
710

711
712
713
714
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
    } else {
      while (elInfo && elInfo->getElement() != child1) 
	elInfo = stack.traverseNext(elInfo);      
715

716
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
717
718
    }

719

720
    return meshChanged;
721
722
  }

723

724
725
726
727
728
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
					     TraverseStack &stack,
					     GeoIndex subObj,
					     int ithObj, 
					     bool reverseMode)
729
  {
730
731
732
733
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");


    // === Test if there are more elements in stack to check with the code. ===
734

735
736
    ElInfo *elInfo = stack.getElInfo();
    if (!elInfo)
737
      return false;
738

739
740
741
742

    // === Test if code contains a leaf element. If this is the case, the ===
    // === current element is finished. Traverse the mesh to the next     ===
    // === coarser element.                                               ===
743
744
745
746
747
748
749

    if (code.isLeafElement()) {
      int level = elInfo->getLevel();

      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
750

751
      return false;
752
    }
753

754
755
756
757
758
759

    bool meshChanged = false;
    Element *el = elInfo->getElement();


    // === If element is leaf (and the code is not), refine the element. ===
760

761
    if (el->isLeaf()) {
762
763
      TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");

764
      el->setMark(1);
765
766
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
767
      refineManager->refineFunction(elInfo);
768
      meshChanged = true;
769
770
    }

771
772
773
774
775
776

    // === Continue fitting the mesh structure code to the children of the ===
    // === current element.                                                ===

    int s0 = el->getSubObjOfChild(0, subObj, ithObj, elInfo->getType());
    int s1 = el->getSubObjOfChild(1, subObj, ithObj, elInfo->getType());
777
778
779
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
780
      std::swap(s0, s1);
781
782
      std::swap(child0, child1);
    }
783

784
785
786
787
788
789
790
    
    // === Traverse left child. ===

    if (s0 != -1) {
      // The edge/face is contained in the left child, therefore fit this
      // child to the mesh structure code.

791
      stack.traverseNext(elInfo);
792
      code.nextElement();
793
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
794
795
      elInfo = stack.getElInfo();
    } else {
796
797
798
799
      // The edge/face is not contained in the left child. Hence we need
      // to traverse through all subelements of the left child until we get
      // the second child of the current element.

800
801
      do {
	elInfo = stack.traverseNext(elInfo);
802
803
804
      } while (elInfo && elInfo->getElement() != child1); 

      TEST_EXIT_DBG(elInfo != NULL)("This should not happen!\n");
805
806
    }  

807
808
809
810
811
812
813
814
815
816
    TEST_EXIT_DBG(elInfo->getElement() == child1)
      ("Got wrong child with idx = %d! Searched for child idx = %d\n",
       elInfo->getElement()->getIndex(), child1->getIndex());


    // === Traverse right child. ===

    if (s1 != -1) {
      // The edge/face is contained in the right child, therefore fit this
      // child to the mesh structure code.
817
818

      code.nextElement();
819
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
820
    } else {
821
822
823
824
      // The edge/face is not contained in the right child. Hence we need
      // to traverse through all subelements of the right child until we have
      // finished traversing the current element with all its subelements.

825
      int level = elInfo->getLevel();
826

827
828
829
830
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
831

832
833

    return meshChanged;
834
835
  }

836
  
837
  void MeshDistributor::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
838
  {    
839
    int vecSize = data.size();
840
    SerUtil::serialize(out, vecSize);
841
842
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
843
      SerUtil::serialize(out, dofIndex);
844
845
846
847
    }
  }


848
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
849
				    std::map<int, const DegreeOfFreedom*> &dofMap)
850
  {
851
    FUNCNAME("MeshDistributor::deserialize()");
852
853

    int vecSize = 0;
854
    SerUtil::deserialize(in, vecSize);
855
    data.clear();
856
857
858
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
859
      SerUtil::deserialize(in, dofIndex);
860
861
862
863
864
865
866
867
868

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

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


869
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)
870
871
  {
    int mapSize = data.size();
872
    SerUtil::serialize(out, mapSize);
873
874
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
875
      SerUtil::serialize(out, rank);
876
877
878
879
880
      serialize(out, it->second);
    }
  }

  
881
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,