MeshDistributor.cc 65.3 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(mesh)("No mesh has been defined for mesh distribution!\n");

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

83
      return;
84
    }
85
    
86

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

92

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

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

116
    ParallelDebug::testAllElements(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
117
#endif
118

119

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

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

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

128

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

    createLocalGlobalNumbering();

133

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

    removeMacroElements();
137
138
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
139

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

    updateDofAdmins();

144

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

Thomas Witkowski's avatar
Thomas Witkowski committed
147
#if (DEBUG != 0)
148
    MSG("AMDiS runs in debug mode, so make some test ...\n");
149
    debug::testSortedDofs(mesh, elMap);
150
151
    ParallelDebug::testInteriorBoundary(*this);
    ParallelDebug::testCommonDofs(*this, true);
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
287
    }

    probStat.push_back(probVec);
  }


  void MeshDistributor::exitParallelization(AdaptInfo *adaptInfo)
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
329
330
331
332
333
334
335
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

336

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

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

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

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

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

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


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

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

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

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

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

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

406

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


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

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

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


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


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

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

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

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

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

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

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

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

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

516

517
      // === Check the boundaries and adapt mesh if necessary. ===
518

519
520
521
522
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

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

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

536
537

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

546
    mesh->dofCompress();
547
    updateLocalGlobalNumbering();
548
549
550
551

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

    createPeriodicMap();
552
553
554
  }

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

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

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

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

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

590
591
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
592

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

596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
	  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);

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

614
615
	  if (b)
	    meshFitTogether = false;
616
 	}
617

618
	i++;
619
620
621
      }
    }

622
    return meshFitTogether;
623
  }
624
625


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

635
636
    TEST_EXIT_DBG(el)("No element given!\n");

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

642
643
644
645
646
    // 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);
647

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

650
    bool meshChanged = false;
651
652
653
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

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

659

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

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

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

676
677
678

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

679
    if (el->isLeaf()) {
680
681

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

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

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

693
      // Code is not leaf, therefore refine the element.
694
      el->setMark(1);
695
696
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
697
      refineManager->refineFunction(elInfo);
698
699
700
701

      MSG("A-NEW ELs %d %d\n", el->getChild(0)->getIndex(), el->getChild(1)->getIndex());

      meshChanged = true;
702
    }
703

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

711
712
713
    // === 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.                                          ===
714

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

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

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

727
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
728
729
    }

730

731
    return meshChanged;
732
733
  }

734

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


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

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

750
751
752
753

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

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

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

762
      return false;
763
    }
764

765
766
767
768
769
770

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


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

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

775
      el->setMark(1);
776
777
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
778
      refineManager->refineFunction(elInfo);
779
      meshChanged = true;
780
781
    }

782
783
784
785
786
787

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

795
796
797
798
799
800
801
    
    // === Traverse left child. ===

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

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

811
812
      do {
	elInfo = stack.traverseNext(elInfo);
813
814
815
      } while (elInfo && elInfo->getElement() != child1); 

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

818
819
820
821
822
823
824
825
826
827
    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.
828
829

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

836
      int level = elInfo->getLevel();
837

838
839
840
841
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
842

843
844

    return meshChanged;
845
846
  }

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


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

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

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

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


880
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)
881
882
  {
    int mapSize = data.size();
883
    SerUtil::serialize(out, mapSize);
884
885
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
886
      SerUtil::serialize(out, rank);
887
888
889
890
891