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

51
52
namespace AMDiS {

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

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

62

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

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

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

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

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

100
101
102
103
104
105
106
107
108
    string partStr = "parmetis";
    GET_PARAMETER(0, name + "->partitioner", &partStr);

    if (partStr == "parmetis") 
      partitioner = new ParMetisPartitioner(&mpiComm);

    if (partStr == "zoltan")
      partitioner = new ZoltanPartitioner(&mpiComm);

Thomas Witkowski's avatar
Thomas Witkowski committed
109
110
111
    if (partStr == "simple")
      partitioner = new SimplePartitioner(&mpiComm);

Thomas Witkowski's avatar
Thomas Witkowski committed
112
113
114
115
    tmp = 0;
    GET_PARAMETER(0, name + "->box partitioning", "%d", &tmp);
    partitioner->setBoxPartitioning(static_cast<bool>(tmp));

116
    TEST_EXIT(partitioner)("Could not create partitioner \"%s\"!\n", partStr.c_str());
117
118
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
119

120
  void MeshDistributor::initParallelization()
121
  {
122
    FUNCNAME("MeshDistributor::initParallelization()");
123
124
125
126

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

127
128
    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");
129

130
131
132
    int a = 0;
    char *b = NULL;
    float zoltanVersion = 0.0;
133
134
    Zoltan_Initialize(a, &b, &zoltanVersion);

135
136
    elObjects.setFeSpace(feSpace);

137
138
139
    // 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).
140
    if (deserialized) {
141
142
      updateMacroElementInfo();

143
      setRankDofs();
144

145
      removePeriodicBoundaryConditions();
146

147
148
      macroElIndexMap.clear();
      macroElIndexTypeMap.clear();
149

150
      for (vector<MacroElement*>::iterator it = allMacroElements.begin();
151
	   it != allMacroElements.end(); ++it) {
152
153
154
155
	macroElIndexMap[(*it)->getIndex()] = (*it)->getElement();
	macroElIndexTypeMap[(*it)->getIndex()] = (*it)->getElType();
      }

156
      return;
157
    }
158

159

160
161
162
163
164
    // 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();

165
166
167
168
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

169
    // create an initial partitioning of the mesh
170
    partitioner->createInitialPartitioning();
171

172
    // set the element weights, which are 1 at the very first begin
173
    setInitialElementWeights();
174
175

    // and now partition the mesh    
176
177
    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
178
    partitioner->createPartitionMap(partitionMap);
179

180

Thomas Witkowski's avatar
Thomas Witkowski committed
181
#if (DEBUG != 0)
182
183
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
184
185
186
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
187

188
      if (writePartMesh > 0) {
189
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex.vtu");
190
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
191
      }
192
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
193
#endif
194

195

196
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
197

198
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
199

200
#if (DEBUG != 0)
201
    ParallelDebug::printBoundaryInfo(*this);
202
203
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
204
205
    MSG("HERE 1\n");

206
207
208
209
210
    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))) {
211
212
213

	  int neighIndex = (*it)->getNeighbour(i)->getIndex();

214
215
216
217
218
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
219
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
220
221
222
223
224
225
226
227
228
	}
      }
    }

    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))) {
229
230
231

	  int neighIndex = (*it)->getNeighbour(i)->getIndex();

232
233
234
235
236
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
237
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
238
239
240
	}
      }
    }
241

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
242
243
244
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
245

246
247
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
248

249
250
251
    // 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.
252
    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
253
254
	 it != mesh->getPeriodicAssociations().end(); ++it)
      const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
255

256
    updateLocalGlobalNumbering();
257

258
259
260
261
262
263
    // === 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();

264
265
    // === If in debug mode, make some tests. ===

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

269
    ParallelDebug::testAllElements(*this);
270
    debug::testSortedDofs(mesh, elMap);
271
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
272

273
    debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
274
275

    MSG("Debug mode tests finished!\n");
276
#endif
277

278
    // === Create periodic DOF mapping, if there are periodic boundaries. ===
279

280
    createPeriodicMap();
281

282
283
284
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
285

286
    // === Global refinements. ===
287
    
Thomas Witkowski's avatar
Thomas Witkowski committed
288
    int globalRefinement = 0;
289
    GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinement);
290
    
Thomas Witkowski's avatar
Thomas Witkowski committed
291
    if (globalRefinement > 0) {
292
      refineManager->globalRefine(mesh, globalRefinement);
293

294
      updateLocalGlobalNumbering();
295
296

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

299
      createPeriodicMap();
300
301
302
303

#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
304
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
305

306

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

309
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
310

311

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

314
    removePeriodicBoundaryConditions();
315
316
  }

317

318
319
320
321
322
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

    if (feSpace != NULL) {
323
      vector<FiniteElemSpace*> feSpaces = probVec->getFeSpaces();
324
325
326
327
328
329
330
331
332
333
334
      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");
335
      TEST_EXIT(mesh->getDofAdmin(0).getNumberOfPreDofs(VERTEX) == 0)
336
337
338
339
340
341
342
343
344
345
346
347
348
	("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());
      }

349
      partitioner->setMesh(mesh);
350
351
352
353
    }

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
354
355
356
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
357
      string filename = "";
358
359
360
361
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
362
363
364
365
366

      int tsModulo = -1;
      GET_PARAMETER(0, probVec->getName() + "->output->write every i-th timestep", 
		    "%d", &tsModulo);
      
367
      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
368
369
      writeSerializationFile = true;
    }    
370
371

    int readSerialization = 0;
372
373
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", 
		  &readSerialization);
374
    if (readSerialization) {
375
      string filename = "";
376
      GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
377
      filename += ".p" + lexical_cast<string>(mpiRank);
378
      MSG("Start deserialization with %s\n", filename.c_str());
379
      ifstream in(filename.c_str());
380
381
382
383

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

384
      probVec->deserialize(in);
385
      in.close();
386
387
      MSG("Deserialization from file: %s\n", filename.c_str());

388
389
390
391
392
393
      filename = "";
      GET_PARAMETER(0, name + "->input->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
394
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
395
396
397
398
399
400
401
402
403
      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;
404
405
406
407
408
409
    }

    probStat.push_back(probVec);
  }


410
  void MeshDistributor::exitParallelization()
411
  {}
412

413
  
414
  void MeshDistributor::testForMacroMesh()
415
  {
416
    FUNCNAME("MeshDistributor::testForMacroMesh()");
417
418
419
420
421
422
423

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
428
429
430
431
432
433
434
435
436
437
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

438

439
  void MeshDistributor::synchVector(DOFVector<double> &vec)
440
  {
441
    StdMpi<vector<double> > stdMpi(mpiComm);
442
443

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
444
	 sendIt != sendDofs.end(); ++sendIt) {
445
      vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
448
      
Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
451
452
453
454

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

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

459
    stdMpi.startCommunication();
460

Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
463
464
465
    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];
  }
466
467


468
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
469
  {
470
    int nComponents = vec.getSize();
471
    StdMpi<vector<double> > stdMpi(mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
472
473
474

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt) {
475
      vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
476
477
478
479
480
481
482
      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])]);
483
484
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
485
      stdMpi.send(sendIt->first, dofs);
486
487
488
    }

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

492
    stdMpi.startCommunication();
493
494

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
495
496
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
497
498

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
499
500
501
502
503
      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++];
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
562
563
564
565
566
567
568
569
570
571
  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()) {
572
	MSG("Unvalid 3D mesh with %d elements on this edge!\n", edgeEls.size());
573
574
575
576
577
578
579
580
581
582
583
584
585
586
	meshIsValid = false;
      }
    }

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


587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
  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);
      }
    }
  }


606
607
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
608
609
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

610
611
612
613
614
615
616
617
    // 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())
618
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
	}
	
	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);
    }    
636
637
638

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
639
640
641
642
  }


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
645
646
647
648
649
650
651
652
653
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


654
  void MeshDistributor::checkMeshChange()
655
  {
656
    FUNCNAME("MeshDistributor::checkMeshChange()");
657

658
659
    double first = MPI::Wtime();

660

661
662
663
664
665
    // === 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);
666

667
    if (recvAllValues == 0)
668
669
      return;

670
671
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
672
   
673
    do {
674
675
      bool meshChanged = false;

676
677
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
678
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
681
682
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
683
 	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
684
685

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

690
691
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it) {
	if (it.getRank() == mpiRank) {
692
693
694
695
696
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE)) {
	    MeshStructure elCode;
	    elCode.init(it->rankObj);
	    
697
698
	    MeshManipulation mm(feSpace);
	    meshChanged |= !(mm.fitElementToMeshCode(elCode, it->neighObj));
699
	  }
700
701
702
703
704
705
	} else {
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
	    allBound[it.getRank()].push_back(*it);	
	}
      }
706

707

708
      // === Check the boundaries and adapt mesh if necessary. ===
709
710
711
712
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

713
      meshChanged |= checkAndAdaptBoundary(allBound);
714
715
716

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

717
      int sendValue = static_cast<int>(meshChanged);
718
719
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
720
721

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

724
#if (DEBUG != 0)
725
    debug::writeMesh(feSpace, -1, debugOutputDir + "mesh");
726
727
#endif

728

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

731
732
    updateLocalGlobalNumbering();

733

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

Thomas Witkowski's avatar
Thomas Witkowski committed
736
    createPeriodicMap();
737

738
739
740
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
741
742


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

745
    nMeshChangesAfterLastRepartitioning++;
746

747
748
749
750

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

751
752
753
754
    if (repartitioningAllowed && 
	nMeshChangesAfterLastRepartitioning >= repartitionIthChange) {
      repartitionMesh();
      nMeshChangesAfterLastRepartitioning = 0;
755
    }
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774

    // === Print imbalance factor. ===

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

    if (mpiRank == 0) {
      int nOverallDofs = 0;
      int maxDofs = numeric_limits<int>::min();
      for (int i = 0; i < mpiSize; i++) {
	nOverallDofs += nDofsInRank[i];
	maxDofs = std::max(maxDofs, nDofsInRank[i]);
      }      
      int avrgDofs = nOverallDofs / mpiSize;
      double imbalance = (static_cast<double>(maxDofs - avrgDofs) /  avrgDofs) * 100.0;

      MSG("Imbalancing factor: %.1f\%\n", imbalance);
    }
775
776
777
  }

  
778
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
779
  {
780
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
781
782
783

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

788
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
789
790
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
791
	elCode.init(boundIt->rankObj);
792
793
794
795
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
796
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
797
    stdMpi.send(sendCodes);
798
799
800
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it)
      stdMpi.recv(it->first);
    stdMpi.startCommunication();
801
 
802
    // === Compare received mesh structure codes. ===
803
    
804
    bool meshChanged = false;
805

806
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
807
     
808
809
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
810
      
811
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
812
	   boundIt != it->second.end(); ++boundIt, i++) {
813

814
815
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
816

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

820
821
	  MeshManipulation mm(feSpace);
	  meshChanged |= mm.fitElementToMeshCode(recvCodes[i], boundIt->rankObj);
822
 	}
823
824
825
      }
    }

826
    return meshChanged;
827
  }
828
829


830
  void MeshDistributor::serialize(ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
831
  {    
832
    int vecSize = data.size();
833
    SerUtil::serialize(out, vecSize);
834
    for (int i = 0; i < vecSize; i++) {
835
      int dofIndex = *(data[i]);
836
      SerUtil::serialize(out, dofIndex);
837
838
839
840
    }
  }


841
842
  void MeshDistributor::deserialize(istream &in, DofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
843
  {
844
    FUNCNAME("MeshDistributor::deserialize()");
845
846

    int vecSize = 0;
847
    SerUtil::deserialize(in, vecSize);
848
    data.clear();
849
850
851
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
852
      SerUtil::deserialize(in, dofIndex);
853
854
855
856
857
858
859
860
861

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

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


862
  void MeshDistributor::serialize(ostream &out, RankToDofContainer &data)
863
864
  {
    int mapSize = data.size();
865
    SerUtil::serialize(out, mapSize);
866
867
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
868
      SerUtil::serialize(out, rank);
869
870
871
872
873
      serialize(out, it->second);
    }
  }

  
874
875
  void MeshDistributor::deserialize(istream &in, RankToDofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
876
  {
877
878
    data.clear();

879
    int mapSize = 0;
880
    SerUtil::deserialize(in, mapSize);
881
882
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
883
      SerUtil::deserialize(in, rank);
884
885
886
887
      deserialize(in, data[rank], dofMap);      
    }
  }

888

889
  void MeshDistributor::setInitialElementWeights() 
890
  {
891
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
892
893

    elemWeights.clear();