MeshDistributor.cc 74.1 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
#include <algorithm>
2
3
#include <iostream>
#include <fstream>
4
5
#include <limits>
#include <stdint.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
6
#include <boost/lexical_cast.hpp>
7
8
#include <boost/filesystem.hpp>

9
10
#include "parallel/MeshDistributor.h"
#include "parallel/ParallelDebug.h"
11
#include "parallel/StdMpi.h"
12
#include "parallel/ParMetisPartitioner.h"
13
14
15
16
17
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
18
19
#include "DOFMatrix.h"
#include "DOFVector.h"
20
#include "SystemVector.h"
21
#include "VtkWriter.h"
22
#include "ElementDofIterator.h"
23
24
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
25
#include "ElementFileWriter.h"
26
#include "VertexVector.h"
27
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28
29
#include "ProblemVec.h"
#include "ProblemInstat.h"
30
#include "Debug.h"
31
#include "MacroInfo.h"
32

33
34
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
35
  using boost::lexical_cast;
36
  using namespace boost::filesystem;
Thomas Witkowski's avatar
Thomas Witkowski committed
37

38
39
40
41
42
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

43

44
45
46
47
48
49
50
51
  MeshDistributor::MeshDistributor(std::string str)
    : probStat(0),
      name(str),
      feSpace(NULL),
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
52
      nRankDofs(0),
53
      nOverallDofs(0),
54
      rstart(0),
55
      deserialized(false),
56
      writeSerializationFile(false),
57
58
      lastMeshChangeIndex(0),
      macroElementStructureConsisten(false)
59
  {
60
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
61

62
63
64
65
66
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
67

68
  void MeshDistributor::initParallelization()
69
  {
70
    FUNCNAME("MeshDistributor::initParallelization()");
71
72
73
74

    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");

75
76
    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");
77

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

85
      return;
86
    }
87

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

94
95
96
97
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

98
99
    // create an initial partitioning of the mesh
    partitioner->createPartitionData();
100

101
    // set the element weights, which are 1 at the very first begin
102
    setInitialElementWeights();
103
104
105
106
107
108

    // and now partition the mesh    
    partitioner->fillCoarsePartitionVec(&oldPartitionVec);
    partitioner->partition(&elemWeights, INITIAL);
    partitioner->fillCoarsePartitionVec(&partitionVec);

109

Thomas Witkowski's avatar
Thomas Witkowski committed
110
#if (DEBUG != 0)
111
112
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
113
114
115
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
116

117
118
      if (writePartMesh > 0) {
	debug::writeElementIndexMesh(mesh, "elementIndex.vtu");
119
	writePartitioningMesh("part.vtu");
120
      } else {
121
	MSG("Skip write part mesh!\n");
122
      }
123
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
124
#endif
125

126

127
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
128

129
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
130

131
#if (DEBUG != 0)
132
    ParallelDebug::printBoundaryInfo(*this);
133
134
#endif

135

136
137
138
139
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

140

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
141
142
143
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
144
145
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
146

147
148
149
150
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

151

152
153
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
154
#if (DEBUG != 0)
155
    MSG("AMDiS runs in debug mode, so make some test ...\n");
156
    ParallelDebug::testAllElements(*this);
157
    debug::testSortedDofs(mesh, elMap);
158
159
    ParallelDebug::testInteriorBoundary(*this);
    ParallelDebug::testCommonDofs(*this, true);
160
    ParallelDebug::testGlobalIndexByCoords(*this);
161
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
162

163
    debug::writeMesh(feSpace, -1, "macro_mesh");   
164
#endif
165

166

167
168
    // === Create periodic dof mapping, if there are periodic boundaries. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
169
    createPeriodicMap();    
170

171
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
172

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

Thomas Witkowski's avatar
Thomas Witkowski committed
176
    if (globalRefinement > 0) {
177
      refineManager->globalRefine(mesh, globalRefinement);
178

179
#if (DEBUG != 0)
180
      debug::writeMesh(feSpace, -1, "gr_mesh");
181
182
#endif

183
      mesh->dofCompress();
184
      updateLocalGlobalNumbering();
185
186

     
187
      // === Update periodic mapping, if there are periodic boundaries. ===
188
      
189
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
190
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
191

192

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

195
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
196

197

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

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

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

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

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

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

274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
      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;
290
291
292
293
294
295
    }

    probStat.push_back(probVec);
  }


296
  void MeshDistributor::exitParallelization()
297
  {}
298

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

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

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

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

323

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

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
338
339
340
341
342
343
344
345
346
347
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

348

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
373
374
375
    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];
  }
376
377


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

Thomas Witkowski's avatar
Thomas Witkowski committed
395
      stdMpi.send(sendIt->first, dofs);
396
397
398
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
402
    stdMpi.startCommunication<double>(MPI_DOUBLE);
403
404

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

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

418

419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
  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);
      }
    }
  }


438
439
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
440
441
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

442
443
444
445
446
447
448
449
    // 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())
450
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
	}
	
	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);
    }    
468
469
470
471
472
473

    // 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();
474
475
476
477
  }


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
478
479
480
481
482
483
484
485
486
487
488
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


489
  void MeshDistributor::checkMeshChange()
490
  {
491
    FUNCNAME("MeshDistributor::checkMeshChange()");    
492

493
494
495
496
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "before_check_mesh");
#endif

497
498
499
500
501
    // === 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);
502

503
    if (recvAllValues == 0)
504
505
      return;

506
507
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
508

509
510
511
512
513
    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.
514
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
515
516

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
517
518
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
519
 	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
520
521

      for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
522
523
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
524
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
525

526
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
527
528
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
529
 	  allBound[it.getRank()].push_back(*it);	
530

531

532
      // === Check the boundaries and adapt mesh if necessary. ===
533
534
535
536
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

537
538
539
540
541
542
543
      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);
544
545
546
547

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

550
#if (DEBUG != 0)
551
    debug::writeMesh(feSpace, -1, "mesh");
552
553
554
555
556
557
#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. ===
558

559
    mesh->dofCompress();
560
    updateLocalGlobalNumbering();
561
562
563
564

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

    createPeriodicMap();
565
566
567
568
569


    // === The mesh has changed, so check if it is required to repartition the mesh. ===

    repartitionMesh();
570
571
572
  }

  
573
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
574
  {
575
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
576
577
578
579
580
581

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

583
584
585
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
586
	elCode.init(boundIt->rankObj);
587
588
589
590
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
591
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
592
    stdMpi.send(sendCodes);
593
    stdMpi.recv(allBound);
594
    stdMpi.startCommunication<uint64_t>(MPI_UNSIGNED_LONG);
595
 
596
    // === Compare received mesh structure codes. ===
597
    
598
599
    bool meshFitTogether = true;

600
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
601
     
602
603
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
604
      
605
606
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
607

608
609
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
610

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

614
615
616
617
618
619
	  bool b = startFitElementToMeshCode(recvCodes[i], 
					     boundIt->rankObj.el,
					     boundIt->rankObj.subObj,
					     boundIt->rankObj.ithObj, 
					     boundIt->rankObj.elType,
					     boundIt->rankObj.reverseMode);
620

621
	  if (b)
622
	    meshFitTogether = false;	  
623
 	}
624

625
	i++;
626
627
628
      }
    }

629
    return meshFitTogether;
630
  }
631
632


633
634
635
636
637
638
  bool MeshDistributor::startFitElementToMeshCode(MeshStructure &code, 
						  Element *el, 
						  GeoIndex subObj,
						  int ithObj, 
						  int elType,
						  bool reverseMode)
639
  {
640
    FUNCNAME("MeshDistributor::startFitElementToMeshCode()");
641

642
643
    TEST_EXIT_DBG(el)("No element given!\n");

644
645
    // If the code is empty, the element does not matter and the function can
    // return without chaning the mesh.
646
647
    if (code.empty())
      return false;
648

649
650
651
652
653
    // 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);
654

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

657
    bool meshChanged = false;
658
659
660
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

661
662
663
664
665
    // Test for reverse mode, in which the left and right children of elements
    // are flipped.
    if (reverseMode)
      traverseFlag |= Mesh::CALL_REVERSE_MODE;    

666

667
668
669
670
671
    // === 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.
672
      TraverseStack stack;
673
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
674
675
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
676

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

679
      meshChanged = fitElementToMeshCode(code, stack, subObj, ithObj, reverseMode);
680
681
      return meshChanged;
    }
682

683
684
685

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

686
    if (el->isLeaf()) {
687
688

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

692
      // Create traverse stack and traverse the mesh to the element.
693
      TraverseStack stack;
694
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
695
696
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
697
698
699

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

700
      // Code is not leaf, therefore refine the element.
701
      el->setMark(1);
702
703
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
704
      refineManager->refineFunction(elInfo);
705
      meshChanged = true;
706
    }
707

708
709
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
710
711
712
713
    if (reverseMode) {
      std::swap(s0, s1);
      std::swap(child0, child1);    
    }
714

715
716
717
    // === 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.                                          ===
718

719
720
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
721

722
723
724
    if (s0 != -1) {
      while (elInfo && elInfo->getElement() != child0)
	elInfo = stack.traverseNext(elInfo);     
725

726
727
728
729
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
    } else {
      while (elInfo && elInfo->getElement() != child1) 
	elInfo = stack.traverseNext(elInfo);      
730

731
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
732
733
    }

734

735
    return meshChanged;
736
737
  }

738

739
740
741
742
743
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
					     TraverseStack &stack,
					     GeoIndex subObj,
					     int ithObj, 
					     bool reverseMode)
744
  {
745
746
747
748
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");


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

750
751
    ElInfo *elInfo = stack.getElInfo();
    if (!elInfo)
752
      return false;
753

754
755
756
757

    // === 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.                                               ===
758
759
760
761
762
763
764

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

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

766
      return false;
767
    }
768

769
770
771
772
773
774

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


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

776
    if (el->isLeaf()) {
777
778
      TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");

779
      el->setMark(1);
780
781
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
782
      refineManager->refineFunction(elInfo);
783
      meshChanged = true;
784
785
    }

786
787
788
789
790
791

    // === 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());
792
793
794
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
795
      std::swap(s0, s1);
796
797
      std::swap(child0, child1);
    }
798

799
800
801
802
803
804
805
    
    // === Traverse left child. ===

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

806
      stack.traverseNext(elInfo);
807
      code.nextElement();
808
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
809
810
      elInfo = stack.getElInfo();
    } else {
811
812
813
814
      // 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.

815
816
      do {
	elInfo = stack.traverseNext(elInfo);
817
818
819
      } while (elInfo && elInfo->getElement() != child1); 

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

822
823
824
825
826
827
828
829
830
831
    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.
832
833

      code.nextElement();
834
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
835
    } else {
836
837
838
839
      // 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.

840
      int level = elInfo->getLevel();
841

842
843
844
845
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
846

847
848

    return meshChanged;
849
850
  }

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


863
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
864
				    std::map<int, const DegreeOfFreedom*> &dofMap)
865
  {
866
    FUNCNAME("MeshDistributor::deserialize()");
867
868

    int vecSize = 0;
869
    SerUtil::deserialize(in, vecSize);
870
    data.clear();