MeshDistributor.cc 63.7 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/MeshPartitioner.h"
26
#include "parallel/ParMetisPartitioner.h"
27
#include "parallel/MpiHelper.h"
28
29
30
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
#include "io/VtkWriter.h"
31
32
33
34
35
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
36
37
#include "DOFMatrix.h"
#include "DOFVector.h"
38
#include "SystemVector.h"
39
#include "ElementDofIterator.h"
40
41
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
42
#include "VertexVector.h"
43
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
44
45
#include "ProblemVec.h"
#include "ProblemInstat.h"
46
#include "Debug.h"
47

48
49
namespace AMDiS {

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

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

59

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
100

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

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

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

111
112
    elObjects.setFeSpace(feSpace);

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

119
      setRankDofs();
120

121
      removePeriodicBoundaryConditions();
122

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

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

132
      return;
133
    }
134

135

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

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

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

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

    // and now partition the mesh    
152
    partitioner->getPartitionMap(oldPartitionMap);
153
154
155
156

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

157
    partitioner->getPartitionMap(partitionMap);
158

159

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

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

174

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

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

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

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
209
    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;
	}
      }
    }
210

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

    removeMacroElements();
214

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
215

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

Thomas Witkowski's avatar
Thomas Witkowski committed
218

219
220
221
    // 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.
222
    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
223
224
	 it != mesh->getPeriodicAssociations().end(); ++it)
      const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
225

226
    updateLocalGlobalNumbering();
227

228

229
230
231
232
233
234
    // === 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();

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

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

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

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

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

251

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

254
    createPeriodicMap();
255

256

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

265
      updateLocalGlobalNumbering();
266
267

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

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

273

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

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

278

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

281
    removePeriodicBoundaryConditions();
282
283
  }

284

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

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

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

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

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

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

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

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

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

    probStat.push_back(probVec);
  }


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

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

405

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

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

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

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

426
    stdMpi.startCommunication();
427

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


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

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

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

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

459
    stdMpi.startCommunication();
460
461

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

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
466
467
468
469
470
      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++];
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
553
  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");
  }


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


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

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

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


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


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

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

627

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

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

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

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

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

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

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

674

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

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

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

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

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

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

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

697
698
    updateLocalGlobalNumbering();

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

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


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

708
    nMeshChangesAfterLastRepartitioning++;
709

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

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

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

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

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

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

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

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

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

766
    return meshChanged;
767
  }
768
769


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


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

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

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

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


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

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

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

828

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

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

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

849
850
	elemWeights[elNum] = elWeight;
      }
851

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

863

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

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

871
872
873
874
875

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

876
877
    int repartitioning = 0;

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

    if (mpiRank == 0) {
      int nOverallDofs = 0;
884
885
      int minDofs = numeric_limits<int>::max();
      int maxDofs = numeric_limits<int>::min();
886
887
888
889
890
891
892
893
894
      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);

895
      if (static_cast<double>(maxDofs) / static_cast<double>(minDofs) > 1.2) 
Thomas Witkowski's avatar