MeshDistributor.cc 66.6 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
#include <boost/filesystem.hpp>
6
#include <limits>
7

8
9
#include "parallel/MeshDistributor.h"
#include "parallel/ParallelDebug.h"
10
#include "parallel/StdMpi.h"
11
#include "parallel/ParMetisPartitioner.h"
12
13
14
15
16
17
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.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

32
33
namespace AMDiS {

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

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

42

43
44
45
46
47
48
49
50
  MeshDistributor::MeshDistributor(std::string str)
    : probStat(0),
      name(str),
      feSpace(NULL),
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
51
      initialPartitionMesh(true),
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
    // create an initial partitioning of the mesh
    partitioner->createPartitionData();
    // set the element weights, which are 1 at the very first begin
97
    setInitialElementWeights();
98
    // and now partition the mesh
99
    partitionMesh();   
100

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

108
109
      if (writePartMesh > 0) {
	debug::writeElementIndexMesh(mesh, "elementIndex.vtu");
110
	writePartitioningMesh("part.vtu");
111
      } else {
112
	MSG("Skip write part mesh!\n");
113
      }
114
    }
115
    ParallelDebug::testAllElements(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
116
#endif
117

118

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

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

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

127

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

    createLocalGlobalNumbering();

132

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

    removeMacroElements();
136
137
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
138

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

    updateDofAdmins();

143

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

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

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

157

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

Thomas Witkowski's avatar
Thomas Witkowski committed
160
    createPeriodicMap();    
161

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

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

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

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

174
      mesh->dofCompress();
175
      updateLocalGlobalNumbering();
176
177

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

183

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

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

188

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

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

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

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

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

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

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

    probStat.push_back(probVec);
  }


287
  void MeshDistributor::exitParallelization()
288
  {}
289

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

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

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

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

314

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

339

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

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

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

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

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

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


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

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

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

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

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

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

409

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


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

433
434
435
436
437
438
439
440
    // 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())
441
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
	}
	
	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);
    }    
459
460
461
462
463
464

    // 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();
465
466
467
468
  }


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


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

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

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

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

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

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

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

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

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

522

523
      // === Check the boundaries and adapt mesh if necessary. ===
524
525
526
527
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

528
529
530
531
532
533
534
      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);
535
536
537
538

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

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

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

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

    createPeriodicMap();
556
557
558
559
560


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

    repartitionMesh();
561
562
563
  }

  
564
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
565
  {
566
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
567
568
569
570
571
572

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

574
575
576
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
577
	elCode.init(boundIt->rankObj);
578
579
580
581
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
582
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
583
    stdMpi.send(sendCodes);
584
585
586
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
587
    // === Compare received mesh structure codes. ===
588
    
589
590
    bool meshFitTogether = true;

591
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
592
     
593
594
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
595
      
596
597
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
598

599
600
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
601

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

605
606
607
608
609
610
	  bool b = startFitElementToMeshCode(recvCodes[i], 
					     boundIt->rankObj.el,
					     boundIt->rankObj.subObj,
					     boundIt->rankObj.ithObj, 
					     boundIt->rankObj.elType,
					     boundIt->rankObj.reverseMode);
611

612
613
	  if (b)
	    meshFitTogether = false;
614
 	}
615

616
	i++;
617
618
619
      }
    }

620
    return meshFitTogether;
621
  }
622
623


624
625
626
627
628
629
  bool MeshDistributor::startFitElementToMeshCode(MeshStructure &code, 
						  Element *el, 
						  GeoIndex subObj,
						  int ithObj, 
						  int elType,
						  bool reverseMode)
630
  {
631
    FUNCNAME("MeshDistributor::startFitElementToMeshCode()");
632

633
634
    TEST_EXIT_DBG(el)("No element given!\n");

635
636
    // If the code is empty, the element does not matter and the function can
    // return without chaning the mesh.
637
638
    if (code.empty())
      return false;
639

640
641
642
643
644
    // 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);
645

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

648
    bool meshChanged = false;
649
650
651
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

652
653
654
655
656
    // Test for reverse mode, in which the left and right children of elements
    // are flipped.
    if (reverseMode)
      traverseFlag |= Mesh::CALL_REVERSE_MODE;    

657

658
659
660
661
662
    // === 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.
663
      TraverseStack stack;
664
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
665
666
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
667

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

670
      meshChanged = fitElementToMeshCode(code, stack, subObj, ithObj, reverseMode);
671
672
      return meshChanged;
    }
673

674
675
676

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

677
    if (el->isLeaf()) {
678
679

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

683
      // Create traverse stack and traverse the mesh to the element.
684
      TraverseStack stack;
685
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
686
687
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
688
689
690

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

691
      // Code is not leaf, therefore refine the element.
692
      el->setMark(1);
693
694
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
695
      refineManager->refineFunction(elInfo);
696
      meshChanged = true;
697
    }
698

699
700
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
701
702
703
704
    if (reverseMode) {
      std::swap(s0, s1);
      std::swap(child0, child1);    
    }
705

706
707
708
    // === 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.                                          ===
709

710
711
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
712

713
714
715
    if (s0 != -1) {
      while (elInfo && elInfo->getElement() != child0)
	elInfo = stack.traverseNext(elInfo);     
716

717
718
719
720
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
    } else {
      while (elInfo && elInfo->getElement() != child1) 
	elInfo = stack.traverseNext(elInfo);      
721

722
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
723
724
    }

725

726
    return meshChanged;
727
728
  }

729

730
731
732
733
734
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
					     TraverseStack &stack,
					     GeoIndex subObj,
					     int ithObj, 
					     bool reverseMode)
735
  {
736
737
738
739
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");


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

741
742
    ElInfo *elInfo = stack.getElInfo();
    if (!elInfo)
743
      return false;
744

745
746
747
748

    // === 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.                                               ===
749
750
751
752
753
754
755

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

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

757
      return false;
758
    }
759

760
761
762
763
764
765

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


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

767
    if (el->isLeaf()) {
768
769
      TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");

770
      el->setMark(1);
771
772
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
773
      refineManager->refineFunction(elInfo);
774
      meshChanged = true;
775
776
    }

777
778
779
780
781
782

    // === 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());
783
784
785
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
786
      std::swap(s0, s1);
787
788
      std::swap(child0, child1);
    }
789

790
791
792
793
794
795
796
    
    // === Traverse left child. ===

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

797
      stack.traverseNext(elInfo);
798
      code.nextElement();
799
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
800
801
      elInfo = stack.getElInfo();
    } else {
802
803
804
805
      // 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.

806
807
      do {
	elInfo = stack.traverseNext(elInfo);
808
809
810
      } while (elInfo && elInfo->getElement() != child1); 

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

813
814
815
816
817
818
819
820
821
822
    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.
823
824

      code.nextElement();
825
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
826
    } else {
827
828
829
830
      // 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.

831
      int level = elInfo->getLevel();
832

833
834
835
836
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
837

838
839

    return meshChanged;
840
841
  }

842
  
843
  void MeshDistributor::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
844
  {    
845
    int vecSize = data.size();
846
    SerUtil::serialize(out, vecSize);
847
848
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
849
      SerUtil::serialize(out, dofIndex);
850
851
852
853
    }
  }


854
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
855
				    std::map<int, const DegreeOfFreedom*> &dofMap)
856
  {
857
    FUNCNAME("MeshDistributor::deserialize()");
858
859

    int vecSize = 0;
860
    SerUtil::deserialize(in, vecSize);
861
    data.clear();
862
863
864
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
865
      SerUtil::deserialize(in, dofIndex);
866
867
868
869
870
871
872
873
874

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

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