MeshDistributor.cc 65.3 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
101
102
103
104
105
106
107
    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
108
109
110
111
    tmp = 0;
    GET_PARAMETER(0, name + "->box partitioning", "%d", &tmp);
    partitioner->setBoxPartitioning(static_cast<bool>(tmp));

112
    TEST_EXIT(partitioner)("Could not create partitioner \"%s\"!\n", partStr.c_str());
113
114
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
115

116
  void MeshDistributor::initParallelization()
117
  {
118
    FUNCNAME("MeshDistributor::initParallelization()");
119
120
121
122

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

123
124
    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");
125

126
127
128
    int a = 0;
    char *b = NULL;
    float zoltanVersion = 0.0;
129
130
    Zoltan_Initialize(a, &b, &zoltanVersion);

131
132
    elObjects.setFeSpace(feSpace);

133
134
135
    // 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).
136
    if (deserialized) {
137
138
      updateMacroElementInfo();

139
      setRankDofs();
140

141
      removePeriodicBoundaryConditions();
142

143
144
      macroElIndexMap.clear();
      macroElIndexTypeMap.clear();
145

146
      for (vector<MacroElement*>::iterator it = allMacroElements.begin();
147
	   it != allMacroElements.end(); ++it) {
148
149
150
151
	macroElIndexMap[(*it)->getIndex()] = (*it)->getElement();
	macroElIndexTypeMap[(*it)->getIndex()] = (*it)->getElType();
      }

152
      return;
153
    }
154

155

156
157
158
159
160
    // 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();

161
162
163
164
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

165
    // create an initial partitioning of the mesh
166
    partitioner->createInitialPartitioning();
167

168
    // set the element weights, which are 1 at the very first begin
169
    setInitialElementWeights();
170
171

    // and now partition the mesh    
172
173
    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
174
    partitioner->createPartitionMap(partitionMap);
175

176

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

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

191

192
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
193

194
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
195

196
#if (DEBUG != 0)
197
    ParallelDebug::printBoundaryInfo(*this);
198
199
#endif

Thomas Witkowski's avatar
Thomas Witkowski committed
200
201
    MSG("HERE 1\n");

202
203
204
205
206
    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))) {
207
208
209

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

210
211
212
213
214
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
215
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
216
217
218
219
220
221
222
223
224
	}
      }
    }

    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))) {
225
226
227

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

228
229
230
231
232
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
233
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
234
235
236
	}
      }
    }
237

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
238
239
240
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
241

242
243
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
244

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

252
    updateLocalGlobalNumbering();
253

254
255
256
257
258
259
    // === 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();

260
261
    // === If in debug mode, make some tests. ===

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

265
    ParallelDebug::testAllElements(*this);
266
    debug::testSortedDofs(mesh, elMap);
267
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
268

269
    debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
270
271

    MSG("Debug mode tests finished!\n");
272
#endif
273

274
    // === Create periodic DOF mapping, if there are periodic boundaries. ===
275

276
    createPeriodicMap();
277

278
279
280
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
281

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

290
      updateLocalGlobalNumbering();
291
292

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

295
      createPeriodicMap();
296
297
298
299

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

302

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

305
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
306

307

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

310
    removePeriodicBoundaryConditions();
311
312
  }

313

314
315
316
317
318
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

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

345
      partitioner->setMesh(mesh);
346
347
348
349
    }

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

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

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

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

380
      probVec->deserialize(in);
381
      in.close();
382
383
      MSG("Deserialization from file: %s\n", filename.c_str());

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

    probStat.push_back(probVec);
  }


406
  void MeshDistributor::exitParallelization()
407
  {}
408

409
  
410
  void MeshDistributor::testForMacroMesh()
411
  {
412
    FUNCNAME("MeshDistributor::testForMacroMesh()");
413
414
415
416
417
418
419

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
424
425
426
427
428
429
430
431
432
433
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

434

435
  void MeshDistributor::synchVector(DOFVector<double> &vec)
436
  {
437
    StdMpi<vector<double> > stdMpi(mpiComm);
438
439

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
453
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
454

455
    stdMpi.startCommunication();
456

Thomas Witkowski's avatar
Thomas Witkowski committed
457
458
459
460
461
    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];
  }
462
463


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
481
      stdMpi.send(sendIt->first, dofs);
482
483
484
    }

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

488
    stdMpi.startCommunication();
489
490

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
491
492
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
493
494

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
495
496
497
498
499
      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++];
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
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
  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");
  }


583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
  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);
      }
    }
  }


602
603
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
604
605
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
635
636
637
638
  }


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


650
  void MeshDistributor::checkMeshChange()
651
  {
652
    FUNCNAME("MeshDistributor::checkMeshChange()");
653

654
655
    double first = MPI::Wtime();

656

657
658
659
660
661
    // === 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);
662

663
    if (recvAllValues == 0)
664
665
      return;

666
667
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
668
   
669
    do {
670
671
      bool meshChanged = false;

672
673
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
674
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
675
676

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

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

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

703

704
      // === Check the boundaries and adapt mesh if necessary. ===
705
706
707
708
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

709
      meshChanged |= checkAndAdaptBoundary(allBound);
710
711
712

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

713
      int sendValue = static_cast<int>(meshChanged);
714
715
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
716
717

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

720
#if (DEBUG != 0)
721
    debug::writeMesh(feSpace, -1, debugOutputDir + "mesh");
722
723
#endif

724

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

727
728
    updateLocalGlobalNumbering();

729

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

Thomas Witkowski's avatar
Thomas Witkowski committed
732
    createPeriodicMap();
733

734
735
736
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
737
738


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

741
    nMeshChangesAfterLastRepartitioning++;
742

743
744
745
746

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

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

    // === 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);
    }
771
772
773
  }

  
774
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
775
  {
776
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
777
778
779

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

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

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

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

810
811
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
812

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

816
817
	  MeshManipulation mm(feSpace);
	  meshChanged |= mm.fitElementToMeshCode(recvCodes[i], boundIt->rankObj);
818
 	}
819
820
821
      }
    }

822
    return meshChanged;
823
  }
824
825


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


837
838
  void MeshDistributor::deserialize(istream &in, DofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
839
  {
840
    FUNCNAME("MeshDistributor::deserialize()");
841
842

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

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

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


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

  
870
871
  void MeshDistributor::deserialize(istream &in, RankToDofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
872
  {
873
874
    data.clear();

875
    int mapSize = 0;
876
    SerUtil::deserialize(in, mapSize);
877
878
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
879
      SerUtil::deserialize(in, rank);
880
881
882
883
      deserialize(in, data[rank], dofMap);      
    }
  }

884

885
  void MeshDistributor::setInitialElementWeights() 
886
  {
887
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
888
889

    elemWeights.clear();
890
      
891
    string filename = "";
892
893
894
895
    GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());