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

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

31
32
namespace AMDiS {

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

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

41

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
66

67
  void MeshDistributor::initParallelization(AdaptInfo *adaptInfo)
68
  {
69
    FUNCNAME("MeshDistributor::initParallelization()");
70
71
72
73

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

74
75
    TEST_EXIT(feSpace)("No FE space has been defined for the mesh distributor!\n");
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
76

77
78
79
    // If the problem has been already read from a file, we need only to set
    // isRankDofs to all matrices and rhs vector and to remove periodic 
    // boundary conditions (if there are some).
80
    if (deserialized) {
81
      setRankDofs();
82
      removePeriodicBoundaryConditions();
83

84
      return;
85
    }
86

87
   
88
89
90
91
92
    // Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
    // only for macro meshes, so it will not work yet if the mesh is already refined
    // in some way.
    testForMacroMesh();

93
94
95
96
97
98
99
    // create an initial partitioning of the mesh
    partitioner->createPartitionData();
    // set the element weights, which are 1 at the very first begin
    setElemWeights(adaptInfo);
    // and now partition the mesh
    partitionMesh(adaptInfo);   

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

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

117

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

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

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

126

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

    createLocalGlobalNumbering();

131

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

    removeMacroElements();
135
136
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
137

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

    updateDofAdmins();

142

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

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

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

155

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

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

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

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

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

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

172
      mesh->dofCompress();
173
      updateLocalGlobalNumbering();
174
175

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

181

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

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

186

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

189
    removePeriodicBoundaryConditions();
190
191
  }

192

193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

    if (feSpace != NULL) {
      std::vector<FiniteElemSpace*> feSpaces = probVec->getFeSpaces();
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	TEST_EXIT(feSpace == feSpaces[i])
	  ("Parallelizaton is not supported for multiple FE spaces!\n");
      }
    } else {
      feSpace = probVec->getFeSpace(0);
      mesh = feSpace->getMesh();
      info = probVec->getInfo();
      
      TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
	("Only meshes with one DOFAdmin are supported!\n");
      TEST_EXIT(mesh->getDOFAdmin(0).getNumberOfPreDOFs(0) == 0)
	("Wrong pre dof number for DOFAdmin!\n");
      
      switch (mesh->getDim()) {
      case 2:
	refineManager = new RefinementManager2d();
	break;
      case 3:
	refineManager = new RefinementManager3d();
	break;
      default:
	ERROR_EXIT("This should not happen for dim = %d!\n", mesh->getDim());
      }

      partitioner = new ParMetisPartitioner(mesh, &mpiComm);
    }

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
229
230
231
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
232
233
234
235
236
      std::string filename = "";
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
237
238
239
240
241

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

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

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

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

263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
      filename = "";
      GET_PARAMETER(0, name + "->input->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
      std::string rankFilename = filename + ".p" + lexical_cast<std::string>(mpiRank);
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
279
280
281
282
283
284
285
    }

    probStat.push_back(probVec);
  }


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

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

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

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

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

312

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

337

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

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

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

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

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

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


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

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

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

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

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

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

407

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


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

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

    // 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();
463
464
465
466
  }


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


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

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

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

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

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

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

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

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

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

517

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

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

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

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

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

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

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

    createPeriodicMap();
554
555
556
  }

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

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

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

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

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

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

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

598
#if 0
599
600
601
602
603
604
605
606
607
608
	  MSG("CHECK EL %d %d %d %d WITH El %d %d %d %d in RANK %d\n",
	      boundIt->rankObj.elIndex,
	      boundIt->rankObj.subObj,
	      boundIt->rankObj.ithObj, 
	      boundIt->rankObj.elType,	      
	      boundIt->neighObj.elIndex,
	      boundIt->neighObj.subObj,
	      boundIt->neighObj.ithObj, 
	      boundIt->neighObj.elType,	      
	      it->first);
609
#endif
610
611
612
613
614
615
616

	  bool b = startFitElementToMeshCode(recvCodes[i], 
					     boundIt->rankObj.el,
					     boundIt->rankObj.subObj,
					     boundIt->rankObj.ithObj, 
					     boundIt->rankObj.elType,
					     boundIt->rankObj.reverseMode);
617

618
619
	  if (b)
	    meshFitTogether = false;
620
 	}
621

622
	i++;
623
624
625
      }
    }

626
    return meshFitTogether;
627
  }
628
629


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

639
640
    TEST_EXIT_DBG(el)("No element given!\n");

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

646
647
648
649
650
    // 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);
651

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

654
    bool meshChanged = false;
655
656
657
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

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

663

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

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

676
      meshChanged = fitElementToMeshCode(code, stack, subObj, ithObj, reverseMode);
677
678
      return meshChanged;
    }
679

680
681
682

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

683
    if (el->isLeaf()) {
684
685

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

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

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

697
      // Code is not leaf, therefore refine the element.
698
      el->setMark(1);
699
700
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
701
      refineManager->refineFunction(elInfo);
702
      MSG("REFINE EL %d %d %d\n", elInfo->getElement()->getIndex(), el->getChild(0)->getIndex(), el->getChild(1)->getIndex());
703
      meshChanged = true;
704
    }
705

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

713
714
715
    // === 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.                                          ===
716

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

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

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

729
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
730
731
    }

732

733
    return meshChanged;
734
735
  }

736

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


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

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

752
753
754
755

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

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

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

764
      return false;
765
    }
766

767
768
769
770
771
772

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


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

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

777
      el->setMark(1);
778
779
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
780
      refineManager->refineFunction(elInfo);
781
      MSG("REFINE EL %d %d %d\n", elInfo->getElement()->getIndex(), el->getChild(0)->getIndex(), el->getChild(1)->getIndex());
782
      meshChanged = true;
783
784
    }

785
786
787
788
789
790

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

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

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

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

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

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

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

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

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

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

846
847

    return meshChanged;
848
849
  }

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


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

    int vecSize = 0;
868
    SerUtil::deserialize(in, vecSize);
869
    data.clear();
870
871
872
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
873
      SerUtil::deserialize(in, dofIndex);
874
875
876
877
878
879
880
881
882

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

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


883
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)