MeshDistributor.cc 65.4 KB
Newer Older
1
//
2
3
4
5
6
7
8
9
10
11
12
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include <algorithm>
14
15
#include <iostream>
#include <fstream>
16
17
#include <limits>
#include <stdint.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
18
#include <boost/lexical_cast.hpp>
19
20
#include <boost/filesystem.hpp>

21
#include "parallel/MeshDistributor.h"
22
#include "parallel/MeshManipulation.h"
23
#include "parallel/ParallelDebug.h"
24
#include "parallel/StdMpi.h"
25
#include "parallel/ParMetisPartitioner.h"
26
#include "parallel/MpiHelper.h"
27
28
29
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
#include "io/VtkWriter.h"
30
31
32
33
34
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
35
36
#include "DOFMatrix.h"
#include "DOFVector.h"
37
#include "SystemVector.h"
38
#include "ElementDofIterator.h"
39
40
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
41
#include "VertexVector.h"
42
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
43
44
#include "ProblemVec.h"
#include "ProblemInstat.h"
45
#include "Debug.h"
46

47
48
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
49
  using boost::lexical_cast;
50
  using namespace boost::filesystem;
51
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
52

53
54
55
56
57
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

58

59
  MeshDistributor::MeshDistributor(string str)
60
61
62
63
64
65
66
    : probStat(0),
      name(str),
      feSpace(NULL),
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
67
      nRankDofs(0),
68
      nOverallDofs(0),
69
      rstart(0),
70
      deserialized(false),
71
      writeSerializationFile(false),
72
      repartitioningAllowed(false),
73
      repartitionIthChange(20),
74
75
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
76
      debugOutputDir(""),
77
      lastMeshChangeIndex(0)
78
  {
79
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
80

81
82
83
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
84
85
86
87

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

89
90
    GET_PARAMETER(0, name + "->debug output dir", &debugOutputDir);
    GET_PARAMETER(0, name + "->repartition ith change", "%d", &repartitionIthChange);
91
92
93
94

    tmp = 0;
    GET_PARAMETER(0, name + "->log main rank", "%d", &tmp);
    Msg::outputMainRank = (tmp > 0);
95
96

    partitioner = new ParMetisPartitioner(&mpiComm);
97
98
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
99

100
  void MeshDistributor::initParallelization()
101
  {
102
    FUNCNAME("MeshDistributor::initParallelization()");
103
104
105
106

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

107
108
    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");
109

110
111
    elObjects.setFeSpace(feSpace);

112
113
114
    // 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).
115
    if (deserialized) {
116
117
      updateMacroElementInfo();

118
      setRankDofs();
119

120
      removePeriodicBoundaryConditions();
121

122
123
      macroElIndexMap.clear();
      macroElIndexTypeMap.clear();
124

125
      for (vector<MacroElement*>::iterator it = allMacroElements.begin();
126
	   it != allMacroElements.end(); ++it) {
127
128
129
130
	macroElIndexMap[(*it)->getIndex()] = (*it)->getElement();
	macroElIndexTypeMap[(*it)->getIndex()] = (*it)->getElType();
      }

131
      return;
132
    }
133

134

135
136
137
138
139
    // 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();

140
141
142
143
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

144
    // create an initial partitioning of the mesh
145
    partitioner->createInitialPartitioning();
146

147
    // set the element weights, which are 1 at the very first begin
148
    setInitialElementWeights();
149
150
151

    // and now partition the mesh    
    partitioner->fillCoarsePartitionVec(&oldPartitionVec);
152
153
154
155

    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");

156
157
    partitioner->fillCoarsePartitionVec(&partitionVec);

158

Thomas Witkowski's avatar
Thomas Witkowski committed
159
#if (DEBUG != 0)
160
161
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
162
163
164
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
165

166
      if (writePartMesh > 0) {
167
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex.vtu");
168
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
169
      }
170
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
171
#endif
172

173

174
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
175

176
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
177

178
#if (DEBUG != 0)
179
    ParallelDebug::printBoundaryInfo(*this);
180
181
#endif

182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
    for (deque<MacroElement*>::iterator it = mesh->firstMacroElement();
	 it != mesh->endOfMacroElements(); ++it) {
      for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
	if ((*it)->getNeighbour(i) && 
	    mesh->isPeriodicAssociation((*it)->getBoundary(i))) {
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
	}
      }
    }

    for (vector<MacroElement*>::iterator it = allMacroElements.begin();
	 it != allMacroElements.end(); ++it) {
      for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
	if ((*it)->getNeighbour(i) && 
	    mesh->isPeriodicAssociation((*it)->getBoundary(i))) {
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
	}
      }
    }
209

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
210
211
212
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
213

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
214

215
216
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
217

218
219
220
    // We have to remove the VertexVectors, which contain periodic assoiciations, 
    // because they are not valid anymore after some macro elements have been removed
    // and the corresponding DOFs were deleted.
221
    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
222
223
	 it != mesh->getPeriodicAssociations().end(); ++it)
      const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
224

225
    updateLocalGlobalNumbering();
226

227

228
229
230
231
232
233
    // === In 3D we have to make some test, if the resulting mesh is valid. If   ===
    // === it is not valid, there is no possiblity yet to fix this problem, just ===
    // === exit with an error message.                                           ===
    
    check3dValidMesh();

234
235
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
236
#if (DEBUG != 0)
237
    MSG("AMDiS runs in debug mode, so make some test ...\n");
238

239
    ParallelDebug::testAllElements(*this);
240
    debug::testSortedDofs(mesh, elMap);
241
242
    ParallelDebug::testInteriorBoundary(*this);
    ParallelDebug::testCommonDofs(*this, true);
243
    ParallelDebug::testGlobalIndexByCoords(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
244

245
    debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
246
247

    MSG("Debug mode tests finished!\n");
248
#endif
249

250

251
    // === Create periodic DOF mapping, if there are periodic boundaries. ===
252

253
    createPeriodicMap();
254

255

256
    // === Global refinements. ===
257
    
Thomas Witkowski's avatar
Thomas Witkowski committed
258
    int globalRefinement = 0;
259
    GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinement);
260
    
Thomas Witkowski's avatar
Thomas Witkowski committed
261
    if (globalRefinement > 0) {
262
      refineManager->globalRefine(mesh, globalRefinement);
263

264
      updateLocalGlobalNumbering();
265
266

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

269
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
270
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
271

272

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

275
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
276

277

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

280
    removePeriodicBoundaryConditions();
281
282
  }

283

284
285
286
287
288
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

    if (feSpace != NULL) {
289
      vector<FiniteElemSpace*> feSpaces = probVec->getFeSpaces();
290
291
292
293
294
295
296
297
298
299
300
      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");
301
      TEST_EXIT(mesh->getDofAdmin(0).getNumberOfPreDofs(VERTEX) == 0)
302
303
304
305
306
307
308
309
310
311
312
313
314
	("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());
      }

315
      partitioner->setMesh(mesh);
316
317
318
319
    }

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
320
321
322
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
323
      string filename = "";
324
325
326
327
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
328
329
330
331
332

      int tsModulo = -1;
      GET_PARAMETER(0, probVec->getName() + "->output->write every i-th timestep", 
		    "%d", &tsModulo);
      
333
      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
334
335
      writeSerializationFile = true;
    }    
336
337

    int readSerialization = 0;
338
339
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", 
		  &readSerialization);
340
    if (readSerialization) {
341
      string filename = "";
342
      GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
343
      filename += ".p" + lexical_cast<string>(mpiRank);
344
      MSG("Start deserialization with %s\n", filename.c_str());
345
      ifstream in(filename.c_str());
346
347
348
349

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

350
      probVec->deserialize(in);
351
      in.close();
352
353
      MSG("Deserialization from file: %s\n", filename.c_str());

354
355
356
357
358
359
      filename = "";
      GET_PARAMETER(0, name + "->input->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
360
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
361
362
363
364
365
366
367
368
369
      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;
370
371
372
373
374
375
    }

    probStat.push_back(probVec);
  }


376
  void MeshDistributor::exitParallelization()
377
  {}
378

379
  
380
  void MeshDistributor::testForMacroMesh()
381
  {
382
    FUNCNAME("MeshDistributor::testForMacroMesh()");
383
384
385
386
387
388
389

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
394
395
396
397
398
399
400
401
402
403
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

404

405
  void MeshDistributor::synchVector(DOFVector<double> &vec)
406
  {
407
    StdMpi<vector<double> > stdMpi(mpiComm);
408
409

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
410
	 sendIt != sendDofs.end(); ++sendIt) {
411
      vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
414
      
Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
417
418
419
420

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

Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
423
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
424

425
    stdMpi.startCommunication();
426

Thomas Witkowski's avatar
Thomas Witkowski committed
427
428
429
430
431
    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];
  }
432
433


434
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
435
  {
436
    int nComponents = vec.getSize();
437
    StdMpi<vector<double> > stdMpi(mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
438
439
440

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt) {
441
      vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
442
443
444
445
446
447
448
      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])]);
449
450
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
451
      stdMpi.send(sendIt->first, dofs);
452
453
454
    }

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

458
    stdMpi.startCommunication();
459
460

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
463
464

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
465
466
467
468
469
      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++];
470
471
472
473
      }
    }
  }

474

475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

    if (mesh->getDim() != 3)
      return;

    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
      for (int i = 0; i < 4; i++) {
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
	allEdges.insert(elObjects.getEdgeLocalMap(elData));
      }

      rankMacroEls.insert(elInfo->getElement()->getIndex());

      elInfo = stack.traverseNext(elInfo);
    }

    bool meshIsValid = true;

    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
      vector<ElementObjectData>& edgeEls = elObjects.getElements(*it);

      TEST_EXIT_DBG(edgeEls.size() > 0)
	("No edge %d/%d in elObjDB!\n", it->first, it->second);

      std::set<int> el0, el1;
      for (unsigned int i = 0; i < edgeEls.size(); i++)
	if (rankMacroEls.count(edgeEls[i].elIndex))
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
            
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

      if (el1.empty())
	continue;
      
      bool found = false;      
      do {
	found = false;
	for (std::set<int>::iterator elIt = el0.begin(); elIt != el0.end(); ++elIt) {
	  for (int i = 0; i < 4; i++) {
	    int neighEl = macroElementNeighbours[*elIt][i];
	    if (neighEl != -1 && el1.count(neighEl)) {
	      el0.insert(neighEl);
	      el1.erase(neighEl);
	      found = true;
	    }
	  }
	  
	  if (found)
	    break;
	}
      } while (found);
      
      if (!el1.empty()) {
	meshIsValid = false;
	break;
      }
    }

    if (!meshIsValid) {
      MSG("Error: mesh is not a valid 3D mesh on this rank!\n");
    }

    int foundError = !meshIsValid;
    mpi::globalAdd(foundError);
    TEST_EXIT(foundError == 0)("Error found on at least on rank!\n");
  }


553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
  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);
      }
    }
  }


572
573
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
574
575
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

576
577
578
579
580
581
582
583
    // 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())
584
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
	}
	
	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);
    }    
602
603
604

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
605
606
607
608
  }


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
609
610
611
612
613
614
615
616
617
618
619
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


620
  void MeshDistributor::checkMeshChange()
621
  {
622
    FUNCNAME("MeshDistributor::checkMeshChange()");
623

624
625
    double first = MPI::Wtime();

626

627
628
629
630
631
    // === 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);
632

633
    if (recvAllValues == 0)
634
635
      return;

636
637
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
638
   
639
    do {
640
641
      bool meshChanged = false;

642
643
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
644
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
645
646

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
647
648
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
649
 	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
650
651

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

656
657
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it) {
	if (it.getRank() == mpiRank) {
658
659
660
661
662
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE)) {
	    MeshStructure elCode;
	    elCode.init(it->rankObj);
	    
663
664
	    MeshManipulation mm(feSpace);
	    meshChanged |= !(mm.fitElementToMeshCode(elCode, it->neighObj));
665
	  }
666
667
668
669
670
671
	} else {
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
	    allBound[it.getRank()].push_back(*it);	
	}
      }
672

673

674
      // === Check the boundaries and adapt mesh if necessary. ===
675
676
677
678
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

679
      meshChanged |= checkAndAdaptBoundary(allBound);
680
681
682

      // === Check on all ranks if at least one rank's mesh has changed. ===

683
      int sendValue = static_cast<int>(meshChanged);
684
685
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
686
687

      MSG("Mesh changed on %d ranks!\n", recvAllValues);
688
    } while (recvAllValues != 0);
689

690
#if (DEBUG != 0)
691
    debug::writeMesh(feSpace, -1, debugOutputDir + "mesh");
692
693
694
#endif

    // === Because the mesh has been changed, update the DOF numbering and mappings. ===
695

696
697
    updateLocalGlobalNumbering();

698
    // === Update periodic mapping, if there are periodic boundaries. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
699
    createPeriodicMap();
700
701
702
703
704

    INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n", 
		  MPI::Wtime() - first);


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

707
    nMeshChangesAfterLastRepartitioning++;
708

709
710
711
712
    if (repartitioningAllowed && 
	nMeshChangesAfterLastRepartitioning >= repartitionIthChange) {
      repartitionMesh();
      nMeshChangesAfterLastRepartitioning = 0;
713
    }
714
715
716
  }

  
717
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
718
  {
719
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
720
721
722

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

727
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
728
729
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
730
	elCode.init(boundIt->rankObj);
731
732
733
734
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
735
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
736
    stdMpi.send(sendCodes);
737
738
739
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it)
      stdMpi.recv(it->first);
    stdMpi.startCommunication();
740
 
741
    // === Compare received mesh structure codes. ===
742
    
743
    bool meshChanged = false;
744

745
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
746
     
747
748
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
749
      
750
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
751
	   boundIt != it->second.end(); ++boundIt, i++) {
752

753
754
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
755

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

759
760
	  MeshManipulation mm(feSpace);
	  meshChanged |= mm.fitElementToMeshCode(recvCodes[i], boundIt->rankObj);
761
 	}
762
763
764
      }
    }

765
    return meshChanged;
766
  }
767
768


769
  void MeshDistributor::serialize(ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
770
  {    
771
    int vecSize = data.size();
772
    SerUtil::serialize(out, vecSize);
773
    for (int i = 0; i < vecSize; i++) {
774
      int dofIndex = *(data[i]);
775
      SerUtil::serialize(out, dofIndex);
776
777
778
779
    }
  }


780
781
  void MeshDistributor::deserialize(istream &in, DofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
782
  {
783
    FUNCNAME("MeshDistributor::deserialize()");
784
785

    int vecSize = 0;
786
    SerUtil::deserialize(in, vecSize);
787
    data.clear();
788
789
790
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
791
      SerUtil::deserialize(in, dofIndex);
792
793
794
795
796
797
798
799
800

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

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


801
  void MeshDistributor::serialize(ostream &out, RankToDofContainer &data)
802
803
  {
    int mapSize = data.size();
804
    SerUtil::serialize(out, mapSize);
805
806
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
807
      SerUtil::serialize(out, rank);
808
809
810
811
812
      serialize(out, it->second);
    }
  }

  
813
814
  void MeshDistributor::deserialize(istream &in, RankToDofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
815
  {
816
817
    data.clear();

818
    int mapSize = 0;
819
    SerUtil::deserialize(in, mapSize);
820
821
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
822
      SerUtil::deserialize(in, rank);
823
824
825
826
      deserialize(in, data[rank], dofMap);      
    }
  }

827

828
  void MeshDistributor::setInitialElementWeights() 
829
  {
830
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
831
832

    elemWeights.clear();
833
      
834
    string filename = "";
835
836
837
838
    GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

839
840
      ifstream infile;
      infile.open(filename.c_str(), ifstream::in);
841
842
843
844
845
846
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;
847

848
849
	elemWeights[elNum] = elWeight;
      }
850

851
      infile.close();
852
    } else {           
853
      TraverseStack stack;
854
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
855
      while (elInfo) {
856
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
857
	elInfo = stack.traverseNext(elInfo);
858
859
860
861
      }
    }
  }

862

863
864
865
866
867
868
869
  void MeshDistributor::repartitionMesh()
  {
    FUNCNAME("MeshDistributor::repartitionMesh()");

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

870
871
872
873
874

    // === First we check if the rank with the maximum number of DOFs has at  ===
    // === least 20% more DOFs than the rank with the minimum number of DOFs. ===
    // === In this case, the mesh will be repartition.                        ===

875
876
    int repartitioning = 0;

877
    vector<int> nDofsInRank(mpiSize);
878
    int nDofs = mesh->getDofAdmin(0).getUsedDofs();
879
880
881
882
    mpiComm.Gather(&nDofs, 1, MPI_INT, &(nDofsInRank[0]), 1, MPI_INT, 0);

    if (mpiRank == 0) {
      int nOverallDofs = 0;
883
884
      int minDofs = numeric_limits<int>::max();
      int maxDofs = numeric_limits<int>::min();
885
886
887
888
889
890
891
892
893
      for (int i = 0; i < mpiSize; i++) {
	nOverallDofs += nDofsInRank[i];
	minDofs = std::min(minDofs, nDofsInRank[i]);
	maxDofs = std::max(maxDofs, nDofsInRank[i]);
      }      
     
      MSG("Overall DOFs: %d    Min DOFs: %d    Max DOFs: %d\n", 
	  nOverallDofs, minDofs, maxDofs);

894
      if (static_cast<double>(maxDofs) / static_cast<double>(minDofs) > 1.2) 
895
896
897
898
899
900
901
902
903
904
905
	repartitioning = 1;

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


    if (repartitioning == 0)
      return;