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

34
35
namespace AMDiS {

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

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

44

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

64
65
66
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
67
68
69
70

    int tmp = 0;
    GET_PARAMETER(0, name + "->repartitioning", "%d", &tmp);
    repartitioningAllowed = (tmp > 0);
71
72
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
73

74
  void MeshDistributor::initParallelization()
75
  {
76
    FUNCNAME("MeshDistributor::initParallelization()");
77
78
79
80

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

81
82
    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");
83

84
85
86
    // 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).
87
    if (deserialized) {
88
      setRankDofs();
89
      removePeriodicBoundaryConditions();
90

91
      return;
92
    }
93

94
   
95
96
97
98
99
    // 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();

100
101
102
103
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

104
105
    // create an initial partitioning of the mesh
    partitioner->createPartitionData();
106

107
    // set the element weights, which are 1 at the very first begin
108
    setInitialElementWeights();
109
110
111
112
113
114

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

115

Thomas Witkowski's avatar
Thomas Witkowski committed
116
#if (DEBUG != 0)
117
118
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
119
120
121
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
122

123
124
      if (writePartMesh > 0) {
	debug::writeElementIndexMesh(mesh, "elementIndex.vtu");
125
	writePartitioningMesh("part.vtu");
126
      } else {
127
	MSG("Skip write part mesh!\n");
128
      }
129
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
130
#endif
131

132

133
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
134

135
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
136

137
#if (DEBUG != 0)
138
    ParallelDebug::printBoundaryInfo(*this);
139
140
#endif

141

142
143
144
145
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

146

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
147
148
149
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
150
151
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
152

153
154
155
156
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

157

158
159
    // === If in debug mode, make some tests. ===

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

169
    debug::writeMesh(feSpace, -1, "macro_mesh");   
170
#endif
171

172

173
174
    // === Create periodic dof mapping, if there are periodic boundaries. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
175
    createPeriodicMap();    
176

177
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
178

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

Thomas Witkowski's avatar
Thomas Witkowski committed
182
    if (globalRefinement > 0) {
183
      refineManager->globalRefine(mesh, globalRefinement);
184

185
#if (DEBUG != 0)
186
      debug::writeMesh(feSpace, -1, "gr_mesh");
187
188
#endif

189
      mesh->dofCompress();
190
      updateLocalGlobalNumbering();
191
192

     
193
      // === Update periodic mapping, if there are periodic boundaries. ===
194
      
195
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
196
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
197

198

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

201
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
202

203

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

206
    removePeriodicBoundaryConditions();
207
208
  }

209

210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
  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");
227
      TEST_EXIT(mesh->getDofAdmin(0).getNumberOfPreDOFs(0) == 0)
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
	("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;
246
247
248
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
249
250
251
252
253
      std::string filename = "";
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
254
255
256
257
258

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

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

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

276
      probVec->deserialize(in);
277
      in.close();
278
279
      MSG("Deserialization from file: %s\n", filename.c_str());

280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
      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;
296
297
298
299
300
301
    }

    probStat.push_back(probVec);
  }


302
  void MeshDistributor::exitParallelization()
303
  {}
304

305
  
306
  void MeshDistributor::updateDofAdmins()
307
  {
308
    FUNCNAME("MeshDistributor::updateDofAdmins()");
309
310

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
311
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDofAdmin(i));
312
313
314

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

323
324
325
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
326
327
328
    }
  }

329

330
  void MeshDistributor::testForMacroMesh()
331
  {
332
    FUNCNAME("MeshDistributor::testForMacroMesh()");
333
334
335
336
337
338
339

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
344
345
346
347
348
349
350
351
352
353
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

354

355
  void MeshDistributor::synchVector(DOFVector<double> &vec)
356
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
357
    StdMpi<std::vector<double> > stdMpi(mpiComm);
358
359

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
360
	 sendIt != sendDofs.end(); ++sendIt) {
361
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
364
      
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
367
368
369
370

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

Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
373
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
374

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

Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
379
380
381
    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];
  }
382
383


384
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
385
  {
386
    int nComponents = vec.getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
387
388
389
390
391
392
393
394
395
396
397
398
    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])]);
399
400
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
401
      stdMpi.send(sendIt->first, dofs);
402
403
404
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
408
    stdMpi.startCommunication<double>(MPI_DOUBLE);
409
410

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
411
412
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
413
414

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
419
      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++];
420
421
422
423
      }
    }
  }

424

425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
  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);
      }
    }
  }


444
445
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
446
447
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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

    // Remove periodic vertex associations
    for (std::map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
	 it != mesh->getPeriodicAssociations().end(); ++it)
478
      const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
479
    mesh->getPeriodicAssociations().clear();
480
481
482
483
  }


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
484
485
486
487
488
489
490
491
492
493
494
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


495
  void MeshDistributor::checkMeshChange()
496
  {
497
    FUNCNAME("MeshDistributor::checkMeshChange()");    
498

499
500
501
502
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "before_check_mesh");
#endif

503
504
505
506
507
    // === 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);
508

509
    if (recvAllValues == 0)
510
511
      return;

512
513
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
514

515
516
517
518
519
    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.
520
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
521
522

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

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

532
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
533
534
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
535
 	  allBound[it.getRank()].push_back(*it);	
536

537

538
      // === Check the boundaries and adapt mesh if necessary. ===
539
540
541
542
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

543
544
545
546
547
548
549
      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);
550
551
552
553

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

556
#if (DEBUG != 0)
557
    debug::writeMesh(feSpace, -1, "mesh");
558
559
560
561
562
563
#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. ===
564

565
    mesh->dofCompress();
566
    updateLocalGlobalNumbering();
567
568
569
570

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

    createPeriodicMap();
571
572
573
574


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

575
576
    if (repartitioningAllowed)
      repartitionMesh();
577
578
579
  }

  
580
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
581
  {
582
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
583
584
585
586
587
588

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

590
591
592
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
593
	elCode.init(boundIt->rankObj);
594
595
596
597
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
598
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
599
    stdMpi.send(sendCodes);
600
    stdMpi.recv(allBound);
601
    stdMpi.startCommunication<uint64_t>(MPI_UNSIGNED_LONG);
602
 
603
    // === Compare received mesh structure codes. ===
604
    
605
606
    bool meshFitTogether = true;

607
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
608
     
609
610
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
611
      
612
613
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
614

615
616
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
617

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

621
622
623
624
625
626
	  bool b = startFitElementToMeshCode(recvCodes[i], 
					     boundIt->rankObj.el,
					     boundIt->rankObj.subObj,
					     boundIt->rankObj.ithObj, 
					     boundIt->rankObj.elType,
					     boundIt->rankObj.reverseMode);
627

628
	  if (b)
629
	    meshFitTogether = false;	  
630
 	}
631

632
	i++;
633
634
635
      }
    }

636
    return meshFitTogether;
637
  }
638
639


640
641
642
643
644
645
  bool MeshDistributor::startFitElementToMeshCode(MeshStructure &code, 
						  Element *el, 
						  GeoIndex subObj,
						  int ithObj, 
						  int elType,
						  bool reverseMode)
646
  {
647
    FUNCNAME("MeshDistributor::startFitElementToMeshCode()");
648

649
650
    TEST_EXIT_DBG(el)("No element given!\n");

651
652
    // If the code is empty, the element does not matter and the function can
    // return without chaning the mesh.
653
654
    if (code.empty())
      return false;
655

656
657
658
659
660
    // 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);
661

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

664
    bool meshChanged = false;
665
666
667
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

668
669
670
671
672
    // Test for reverse mode, in which the left and right children of elements
    // are flipped.
    if (reverseMode)
      traverseFlag |= Mesh::CALL_REVERSE_MODE;    

673

674
675
676
677
678
    // === 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.
679
      TraverseStack stack;
680
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
681
682
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
683

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

686
      meshChanged = fitElementToMeshCode(code, stack, subObj, ithObj, reverseMode);
687
688
      return meshChanged;
    }
689

690
691
692

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

693
    if (el->isLeaf()) {
694
695

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

699
      // Create traverse stack and traverse the mesh to the element.
700
      TraverseStack stack;
701
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
702
703
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
704
705
706

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

707
      // Code is not leaf, therefore refine the element.
708
      el->setMark(1);
709
710
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
711
      refineManager->refineFunction(elInfo);
712
      meshChanged = true;
713
    }
714

715
716
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
717
718
719
720
    if (reverseMode) {
      std::swap(s0, s1);
      std::swap(child0, child1);    
    }
721

722
723
724
    // === 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.                                          ===
725

726
727
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
728

729
730
731
    if (s0 != -1) {
      while (elInfo && elInfo->getElement() != child0)
	elInfo = stack.traverseNext(elInfo);     
732

733
734
735
736
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
    } else {
      while (elInfo && elInfo->getElement() != child1) 
	elInfo = stack.traverseNext(elInfo);      
737

738
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
739
740
    }

741

742
    return meshChanged;
743
744
  }

745

746
747
748
749
750
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
					     TraverseStack &stack,
					     GeoIndex subObj,
					     int ithObj, 
					     bool reverseMode)
751
  {
752
753
754
755
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");


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

757
758
    ElInfo *elInfo = stack.getElInfo();
    if (!elInfo)
759
      return false;
760

761
762
763
764

    // === 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.                                               ===
765
766
767
768
769
770
771

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

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

773
      return false;
774
    }
775

776
777
778
779
780
781

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


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

783
    if (el->isLeaf()) {
784
785
      TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");

786
      el->setMark(1);
787
788
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
789
      refineManager->refineFunction(elInfo);
790
      meshChanged = true;
791
792
    }

793
794
795
796
797
798

    // === 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());
799
800
801
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
802
      std::swap(s0, s1);
803
804
      std::swap(child0, child1);
    }
805

806
807
808
809
810
811
812
    
    // === Traverse left child. ===

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

813
      stack.traverseNext(elInfo);
814
      code.nextElement();
815
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
816
817
      elInfo = stack.getElInfo();
    } else {
818
819
820
821
      // 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.

822
823
      do {
	elInfo = stack.traverseNext(elInfo);
824
825
826
      } while (elInfo && elInfo->getElement() != child1); 

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

829
830
831
832
833
834
835
836
837
838
    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.
839
840

      code.nextElement();
841
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
842
    } else {
843
844
845
846
      // 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.

847
      int level = elInfo->getLevel();
848

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

854
855

    return meshChanged;
856
857
  }

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


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

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

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

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


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

  
903
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
904
				    std::map<int, const DegreeOfFreedom*> &dofMap)
905
  {
906
907
    data.clear();

908
    int mapSize = 0;
909
    SerUtil::deserialize(in, mapSize);
910
911
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
912
      SerUtil::deserialize(in, rank);
913
914
915
916
      deserialize(in, data[rank], dofMap);      
    }
  }

917

918
  void MeshDistributor::setInitialElementWeights() 
919
  {
920
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
921
922

    elemWeights.clear();
923
924
925
926
927
928
929
930
931
932
933
934
935
936
      
    std::string filename = "";
    GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

      std::ifstream infile;
      infile.open(filename.c_str(), std::ifstream::in);
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;
937

938
939
	elemWeights[elNum] = elWeight;
      }
940

941
      infile.close();
942
    } else {           
943
      TraverseStack stack;
944
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
945
      while (elInfo) {
946
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
947
	elInfo = stack.traverseNext(elInfo);
948
949
950
951
      }
    }
  }

952

953
954
955
956
957
958
959
  void MeshDistributor::repartitionMesh()
  {
    FUNCNAME("MeshDistributor::repartitionMesh()");

    TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
      ("Only meshes with one DOFAdmin are supported!\n");

960
    int nDofs = mesh->getDofAdmin(0).getUsedDofs();
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
    int repartitioning = 0;

    std::vector<int> nDofsInRank(mpiSize);
    mpiComm.Gather(&nDofs, 1, MPI_INT, &(nDofsInRank[0]), 1, MPI_INT, 0);

    if (mpiRank == 0) {
      int nOverallDofs = 0;
      int minDofs = std::numeric_limits<int>::max();
      int maxDofs = std::numeric_limits<int>::min();
      for (int i = 0; i < mpiSize; i++) {
	nOverallDofs += nDofsInRank[i];
	minDofs = std::min(minDofs, nDofsInRank[i]);
	maxDofs = std::max(maxDofs, nDofsInRank[i]);
      }      
     
      MSG("Overall DOFs: %d    Min DOFs: %d    Max DOFs: %d\n", 
	  nOverallDofs, minDofs, maxDofs);

      if (static_cast<double>(maxDofs) / static_cast<double>(minDofs)) 
	repartitioning = 1;

      mpiComm.Bcast(&repartitioning, 1, MPI_INT, 0);
    } else {
      mpiComm.Bcast(&repartitioning, 1, MPI_INT, 0);
    }


    if (repartitioning == 0)
      return;

991
992
993
994
995
996
    DOFVector<double> tmpa(feSpace, "tmp");
    tmpa.set(mpiRank);
    VtkWriter::writeFile(tmpa, "before-repartition.vtu");

    MSG("USED-SIZE A: %d\n", mesh->getDofAdmin(0).getUsedDofs());

997
998
999
1000
    elemWeights.clear();

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);