MeshDistributor.cc 63.9 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
#include <boost/filesystem.hpp>
20
#include <zoltan_cpp.h>
21

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

50
51
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
52
  using boost::lexical_cast;
53
  using namespace boost::filesystem;
54
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
55

56
57
58
59
60
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

61

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

84
85
86
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
87
88
89
90

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

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

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

99
100
    //    partitioner = new ParMetisPartitioner(&mpiComm);
    partitioner = new ZoltanPartitioner(&mpiComm);
101
102
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
103

104
  void MeshDistributor::initParallelization()
105
  {
106
    FUNCNAME("MeshDistributor::initParallelization()");
107
108
109
110

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

111
112
    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");
113

114
115
116
117
118
    int a;
    char *b;
    float zoltanVersion;
    Zoltan_Initialize(a, &b, &zoltanVersion);

119
120
    elObjects.setFeSpace(feSpace);

121
122
123
    // 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).
124
    if (deserialized) {
125
126
      updateMacroElementInfo();

127
      setRankDofs();
128

129
      removePeriodicBoundaryConditions();
130

131
132
      macroElIndexMap.clear();
      macroElIndexTypeMap.clear();
133

134
      for (vector<MacroElement*>::iterator it = allMacroElements.begin();
135
	   it != allMacroElements.end(); ++it) {
136
137
138
139
	macroElIndexMap[(*it)->getIndex()] = (*it)->getElement();
	macroElIndexTypeMap[(*it)->getIndex()] = (*it)->getElType();
      }

140
      return;
141
    }
142

143

144
145
146
147
148
    // 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();

149
150
151
152
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

153
    // create an initial partitioning of the mesh
154
    partitioner->createInitialPartitioning();
155

156
    // set the element weights, which are 1 at the very first begin
157
    setInitialElementWeights();
158
159

    // and now partition the mesh    
160
    partitioner->getPartitionMap(oldPartitionMap);
161
162
163
164

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

165
    partitioner->getPartitionMap(partitionMap);
166

167

Thomas Witkowski's avatar
Thomas Witkowski committed
168
#if (DEBUG != 0)
169
170
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
171
172
173
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
174

175
      if (writePartMesh > 0) {
176
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex.vtu");
177
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
178
      }
179
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
180
#endif
181

182

183
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
184

185
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
186

187
#if (DEBUG != 0)
188
    ParallelDebug::printBoundaryInfo(*this);
189
190
#endif

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
    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;
	}
      }
    }
218

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
219
220
221
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
222

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
223

224
225
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
226

227
228
229
    // 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.
230
    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
231
232
	 it != mesh->getPeriodicAssociations().end(); ++it)
      const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
233

234
    updateLocalGlobalNumbering();
235

236

237
238
239
240
241
242
    // === 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();

243
244
    // === If in debug mode, make some tests. ===

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

248
    ParallelDebug::testAllElements(*this);
249
    debug::testSortedDofs(mesh, elMap);
250
251
    ParallelDebug::testInteriorBoundary(*this);
    ParallelDebug::testCommonDofs(*this, true);
252
    ParallelDebug::testGlobalIndexByCoords(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
253

254
    debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
255
256

    MSG("Debug mode tests finished!\n");
257
#endif
258

259

260
    // === Create periodic DOF mapping, if there are periodic boundaries. ===
261

262
    createPeriodicMap();
263

264

265
    // === Global refinements. ===
266
    
Thomas Witkowski's avatar
Thomas Witkowski committed
267
    int globalRefinement = 0;
268
    GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinement);
269
    
Thomas Witkowski's avatar
Thomas Witkowski committed
270
    if (globalRefinement > 0) {
271
      refineManager->globalRefine(mesh, globalRefinement);
272

273
      updateLocalGlobalNumbering();
274
275

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

278
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
280

281

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

284
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
285

286

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

289
    removePeriodicBoundaryConditions();
290
291
  }

292

293
294
295
296
297
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

    if (feSpace != NULL) {
298
      vector<FiniteElemSpace*> feSpaces = probVec->getFeSpaces();
299
300
301
302
303
304
305
306
307
308
309
      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");
310
      TEST_EXIT(mesh->getDofAdmin(0).getNumberOfPreDofs(VERTEX) == 0)
311
312
313
314
315
316
317
318
319
320
321
322
323
	("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());
      }

324
      partitioner->setMesh(mesh);
325
326
327
328
    }

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
329
330
331
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
332
      string filename = "";
333
334
335
336
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
337
338
339
340
341

      int tsModulo = -1;
      GET_PARAMETER(0, probVec->getName() + "->output->write every i-th timestep", 
		    "%d", &tsModulo);
      
342
      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
343
344
      writeSerializationFile = true;
    }    
345
346

    int readSerialization = 0;
347
348
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", 
		  &readSerialization);
349
    if (readSerialization) {
350
      string filename = "";
351
      GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
352
      filename += ".p" + lexical_cast<string>(mpiRank);
353
      MSG("Start deserialization with %s\n", filename.c_str());
354
      ifstream in(filename.c_str());
355
356
357
358

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

359
      probVec->deserialize(in);
360
      in.close();
361
362
      MSG("Deserialization from file: %s\n", filename.c_str());

363
364
365
366
367
368
      filename = "";
      GET_PARAMETER(0, name + "->input->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
369
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
370
371
372
373
374
375
376
377
378
      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;
379
380
381
382
383
384
    }

    probStat.push_back(probVec);
  }


385
  void MeshDistributor::exitParallelization()
386
  {}
387

388
  
389
  void MeshDistributor::testForMacroMesh()
390
  {
391
    FUNCNAME("MeshDistributor::testForMacroMesh()");
392
393
394
395
396
397
398

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
403
404
405
406
407
408
409
410
411
412
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

413

414
  void MeshDistributor::synchVector(DOFVector<double> &vec)
415
  {
416
    StdMpi<vector<double> > stdMpi(mpiComm);
417
418

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
419
	 sendIt != sendDofs.end(); ++sendIt) {
420
      vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
423
      
Thomas Witkowski's avatar
Thomas Witkowski committed
424
425
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
426
427
428
429

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

Thomas Witkowski's avatar
Thomas Witkowski committed
430
431
432
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
433

434
    stdMpi.startCommunication();
435

Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438
439
440
    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];
  }
441
442


443
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
444
  {
445
    int nComponents = vec.getSize();
446
    StdMpi<vector<double> > stdMpi(mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
449

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt) {
450
      vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
453
454
455
456
457
      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])]);
458
459
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
460
      stdMpi.send(sendIt->first, dofs);
461
462
463
    }

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

467
    stdMpi.startCommunication();
468
469

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
472
473

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
474
475
476
477
478
      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++];
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
553
554
555
556
557
558
559
560
561
  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");
  }


562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
  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);
      }
    }
  }


581
582
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
583
584
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

585
586
587
588
589
590
591
592
    // 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())
593
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
	}
	
	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);
    }    
611
612
613

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
614
615
616
617
  }


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
618
619
620
621
622
623
624
625
626
627
628
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


629
  void MeshDistributor::checkMeshChange()
630
  {
631
    FUNCNAME("MeshDistributor::checkMeshChange()");
632

633
634
    double first = MPI::Wtime();

635

636
637
638
639
640
    // === 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);
641

642
    if (recvAllValues == 0)
643
644
      return;

645
646
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
647
   
648
    do {
649
650
      bool meshChanged = false;

651
652
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
653
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
654
655

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
656
657
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
658
 	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
659
660

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

665
666
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it) {
	if (it.getRank() == mpiRank) {
667
668
669
670
671
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE)) {
	    MeshStructure elCode;
	    elCode.init(it->rankObj);
	    
672
673
	    MeshManipulation mm(feSpace);
	    meshChanged |= !(mm.fitElementToMeshCode(elCode, it->neighObj));
674
	  }
675
676
677
678
679
680
	} else {
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
	    allBound[it.getRank()].push_back(*it);	
	}
      }
681

682

683
      // === Check the boundaries and adapt mesh if necessary. ===
684
685
686
687
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

688
      meshChanged |= checkAndAdaptBoundary(allBound);
689
690
691

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

692
      int sendValue = static_cast<int>(meshChanged);
693
694
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
695
696

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

699
#if (DEBUG != 0)
700
    debug::writeMesh(feSpace, -1, debugOutputDir + "mesh");
701
702
703
#endif

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

705
706
    updateLocalGlobalNumbering();

707
    // === Update periodic mapping, if there are periodic boundaries. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
708
    createPeriodicMap();
709
710
711
712
713

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


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

716
    nMeshChangesAfterLastRepartitioning++;
717

718
719
720
721
    if (repartitioningAllowed && 
	nMeshChangesAfterLastRepartitioning >= repartitionIthChange) {
      repartitionMesh();
      nMeshChangesAfterLastRepartitioning = 0;
722
    }
723
724
725
  }

  
726
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
727
  {
728
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
729
730
731

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

736
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
737
738
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
739
	elCode.init(boundIt->rankObj);
740
741
742
743
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
744
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
745
    stdMpi.send(sendCodes);
746
747
748
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it)
      stdMpi.recv(it->first);
    stdMpi.startCommunication();
749
 
750
    // === Compare received mesh structure codes. ===
751
    
752
    bool meshChanged = false;
753

754
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
755
     
756
757
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
758
      
759
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
760
	   boundIt != it->second.end(); ++boundIt, i++) {
761

762
763
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
764

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

768
769
	  MeshManipulation mm(feSpace);
	  meshChanged |= mm.fitElementToMeshCode(recvCodes[i], boundIt->rankObj);
770
 	}
771
772
773
      }
    }

774
    return meshChanged;
775
  }
776
777


778
  void MeshDistributor::serialize(ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
779
  {    
780
    int vecSize = data.size();
781
    SerUtil::serialize(out, vecSize);
782
    for (int i = 0; i < vecSize; i++) {
783
      int dofIndex = *(data[i]);
784
      SerUtil::serialize(out, dofIndex);
785
786
787
788
    }
  }


789
790
  void MeshDistributor::deserialize(istream &in, DofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
791
  {
792
    FUNCNAME("MeshDistributor::deserialize()");
793
794

    int vecSize = 0;
795
    SerUtil::deserialize(in, vecSize);
796
    data.clear();
797
798
799
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
800
      SerUtil::deserialize(in, dofIndex);
801
802
803
804
805
806
807
808
809

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

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


810
  void MeshDistributor::serialize(ostream &out, RankToDofContainer &data)
811
812
  {
    int mapSize = data.size();
813
    SerUtil::serialize(out, mapSize);
814
815
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
816
      SerUtil::serialize(out, rank);
817
818
819
820
821
      serialize(out, it->second);
    }
  }

  
822
823
  void MeshDistributor::deserialize(istream &in, RankToDofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
824
  {
825
826
    data.clear();

827
    int mapSize = 0;
828
    SerUtil::deserialize(in, mapSize);
829
830
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
831
      SerUtil::deserialize(in, rank);
832
833
834
835
      deserialize(in, data[rank], dofMap);      
    }
  }

836

837
  void MeshDistributor::setInitialElementWeights() 
838
  {
839
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
840
841

    elemWeights.clear();
842
      
843
    string filename = "";
844
845
846
847
    GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

848
849
      ifstream infile;
      infile.open(filename.c_str(), ifstream::in);
850
851
852
853
854
855
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;
856

857
858
	elemWeights[elNum] = elWeight;
      }
859

860
      infile.close();
861
    } else {           
862
      TraverseStack stack;
863
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
864
      while (elInfo) {
865
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
866
	elInfo = stack.traverseNext(elInfo);
867
868
869
870
      }
    }
  }

871

872
873
874
875
876
877
878
  void MeshDistributor::repartitionMesh()
  {
    FUNCNAME("MeshDistributor::repartitionMesh()");

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

879
880
881
882
883

    // === 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.                        ===

884
885
    int repartitioning = 0;

886
    vector<int> nDofsInRank(mpiSize);
887
    int nDofs = mesh->getDofAdmin(0).getUsedDofs();
888
889
890
891
    mpiComm.Gather(&nDofs, 1, MPI_INT, &(nDofsInRank[0]), 1, MPI_INT, 0);

    if (mpiRank == 0) {
      int nOverallDofs = 0;