MeshDistributor.cc 65 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
108
    string partStr = "parmetis";
    GET_PARAMETER(0, name + "->partitioner", &partStr);

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

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

    TEST_EXIT(partitioner)("Could not create partitioner \"%s\"!\n", partStr.c_str());
109
110
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
111

112
  void MeshDistributor::initParallelization()
113
  {
114
    FUNCNAME("MeshDistributor::initParallelization()");
115
116
117
118

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

119
120
    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");
121

122
123
124
    int a = 0;
    char *b = NULL;
    float zoltanVersion = 0.0;
125
126
    Zoltan_Initialize(a, &b, &zoltanVersion);

127
128
    elObjects.setFeSpace(feSpace);

129
130
131
    // 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).
132
    if (deserialized) {
133
134
      updateMacroElementInfo();

135
      setRankDofs();
136

137
      removePeriodicBoundaryConditions();
138

139
140
      macroElIndexMap.clear();
      macroElIndexTypeMap.clear();
141

142
      for (vector<MacroElement*>::iterator it = allMacroElements.begin();
143
	   it != allMacroElements.end(); ++it) {
144
145
146
147
	macroElIndexMap[(*it)->getIndex()] = (*it)->getElement();
	macroElIndexTypeMap[(*it)->getIndex()] = (*it)->getElType();
      }

148
      return;
149
    }
150

151

152
153
154
155
156
    // 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();

157
158
159
160
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

161
    // create an initial partitioning of the mesh
162
    partitioner->createInitialPartitioning();
163

164
    // set the element weights, which are 1 at the very first begin
165
    setInitialElementWeights();
166
167

    // and now partition the mesh    
168
169
    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
170
    partitioner->createPartitionMap(partitionMap);
171

172

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

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

187

188
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
189

190
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
191

192
#if (DEBUG != 0)
193
    ParallelDebug::printBoundaryInfo(*this);
194
195
#endif

196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
    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;
	}
      }
    }
223

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
224
225
226
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
227

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
228

229
230
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
231

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

239
    updateLocalGlobalNumbering();
240

241

242
243
244
245
246
247
    // === 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();

248
249
    // === If in debug mode, make some tests. ===

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

253
    ParallelDebug::testAllElements(*this);
254
    debug::testSortedDofs(mesh, elMap);
255
256
    ParallelDebug::testInteriorBoundary(*this);
    ParallelDebug::testCommonDofs(*this, true);
257
    ParallelDebug::testGlobalIndexByCoords(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
258

259
    debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
260
261

    MSG("Debug mode tests finished!\n");
262
#endif
263

264
    // === Create periodic DOF mapping, if there are periodic boundaries. ===
265

266
    createPeriodicMap();
267

268
269
270
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
271

272
    // === Global refinements. ===
273
    
Thomas Witkowski's avatar
Thomas Witkowski committed
274
    int globalRefinement = 0;
275
    GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinement);
276
    
Thomas Witkowski's avatar
Thomas Witkowski committed
277
    if (globalRefinement > 0) {
278
      refineManager->globalRefine(mesh, globalRefinement);
279

280
      updateLocalGlobalNumbering();
281
282

     
283
284
      // === Update periodic mapping, if there are periodic boundaries. ===     

285
      createPeriodicMap();
286
287
288
289

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

292

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

295
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
296

297

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

300
    removePeriodicBoundaryConditions();
301
302
  }

303

304
305
306
307
308
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

    if (feSpace != NULL) {
309
      vector<FiniteElemSpace*> feSpaces = probVec->getFeSpaces();
310
311
312
313
314
315
316
317
318
319
320
      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");
321
      TEST_EXIT(mesh->getDofAdmin(0).getNumberOfPreDofs(VERTEX) == 0)
322
323
324
325
326
327
328
329
330
331
332
333
334
	("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());
      }

335
      partitioner->setMesh(mesh);
336
337
338
339
    }

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
340
341
342
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
343
      string filename = "";
344
345
346
347
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
348
349
350
351
352

      int tsModulo = -1;
      GET_PARAMETER(0, probVec->getName() + "->output->write every i-th timestep", 
		    "%d", &tsModulo);
      
353
      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
354
355
      writeSerializationFile = true;
    }    
356
357

    int readSerialization = 0;
358
359
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", 
		  &readSerialization);
360
    if (readSerialization) {
361
      string filename = "";
362
      GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
363
      filename += ".p" + lexical_cast<string>(mpiRank);
364
      MSG("Start deserialization with %s\n", filename.c_str());
365
      ifstream in(filename.c_str());
366
367
368
369

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

370
      probVec->deserialize(in);
371
      in.close();
372
373
      MSG("Deserialization from file: %s\n", filename.c_str());

374
375
376
377
378
379
      filename = "";
      GET_PARAMETER(0, name + "->input->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
380
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
381
382
383
384
385
386
387
388
389
      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;
390
391
392
393
394
395
    }

    probStat.push_back(probVec);
  }


396
  void MeshDistributor::exitParallelization()
397
  {}
398

399
  
400
  void MeshDistributor::testForMacroMesh()
401
  {
402
    FUNCNAME("MeshDistributor::testForMacroMesh()");
403
404
405
406
407
408
409

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
414
415
416
417
418
419
420
421
422
423
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

424

425
  void MeshDistributor::synchVector(DOFVector<double> &vec)
426
  {
427
    StdMpi<vector<double> > stdMpi(mpiComm);
428
429

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
430
	 sendIt != sendDofs.end(); ++sendIt) {
431
      vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
434
      
Thomas Witkowski's avatar
Thomas Witkowski committed
435
436
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
437
438
439
440

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

Thomas Witkowski's avatar
Thomas Witkowski committed
441
442
443
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
444

445
    stdMpi.startCommunication();
446

Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
449
450
451
    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];
  }
452
453


454
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
455
  {
456
    int nComponents = vec.getSize();
457
    StdMpi<vector<double> > stdMpi(mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
460

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt) {
461
      vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
462
463
464
465
466
467
468
      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])]);
469
470
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
471
      stdMpi.send(sendIt->first, dofs);
472
473
474
    }

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

478
    stdMpi.startCommunication();
479
480

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
481
482
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
483
484

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
485
486
487
488
489
      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++];
490
491
492
493
      }
    }
  }

494

495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
  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");
  }


573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
  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);
      }
    }
  }


592
593
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
594
595
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

596
597
598
599
600
601
602
603
    // 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())
604
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
	}
	
	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);
    }    
622
623
624

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
625
626
627
628
  }


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
631
632
633
634
635
636
637
638
639
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


640
  void MeshDistributor::checkMeshChange()
641
  {
642
    FUNCNAME("MeshDistributor::checkMeshChange()");
643

644
645
    double first = MPI::Wtime();

646

647
648
649
650
651
    // === 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);
652

653
    if (recvAllValues == 0)
654
655
      return;

656
657
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
658
   
659
    do {
660
661
      bool meshChanged = false;

662
663
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
664
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
665
666

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
667
668
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
669
 	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
670
671

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

676
677
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it) {
	if (it.getRank() == mpiRank) {
678
679
680
681
682
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE)) {
	    MeshStructure elCode;
	    elCode.init(it->rankObj);
	    
683
684
	    MeshManipulation mm(feSpace);
	    meshChanged |= !(mm.fitElementToMeshCode(elCode, it->neighObj));
685
	  }
686
687
688
689
690
691
	} else {
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
	    allBound[it.getRank()].push_back(*it);	
	}
      }
692

693

694
      // === Check the boundaries and adapt mesh if necessary. ===
695
696
697
698
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

699
      meshChanged |= checkAndAdaptBoundary(allBound);
700
701
702

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

703
      int sendValue = static_cast<int>(meshChanged);
704
705
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
706
707

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

710
#if (DEBUG != 0)
711
    debug::writeMesh(feSpace, -1, debugOutputDir + "mesh");
712
713
#endif

714

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

717
718
    updateLocalGlobalNumbering();

719

720
    // === Update periodic mapping, if there are periodic boundaries. ===
721

Thomas Witkowski's avatar
Thomas Witkowski committed
722
    createPeriodicMap();
723

724
725
726
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
727
728


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

731
    nMeshChangesAfterLastRepartitioning++;
732

733
734
735
736

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

737
738
739
740
    if (repartitioningAllowed && 
	nMeshChangesAfterLastRepartitioning >= repartitionIthChange) {
      repartitionMesh();
      nMeshChangesAfterLastRepartitioning = 0;
741
    }
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760

    // === 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);
    }
761
762
763
  }

  
764
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
765
  {
766
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
767
768
769

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

774
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
775
776
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
777
	elCode.init(boundIt->rankObj);
778
779
780
781
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
782
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
783
    stdMpi.send(sendCodes);
784
785
786
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it)
      stdMpi.recv(it->first);
    stdMpi.startCommunication();
787
 
788
    // === Compare received mesh structure codes. ===
789
    
790
    bool meshChanged = false;
791

792
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
793
     
794
795
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
796
      
797
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
798
	   boundIt != it->second.end(); ++boundIt, i++) {
799

800
801
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
802

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

806
807
	  MeshManipulation mm(feSpace);
	  meshChanged |= mm.fitElementToMeshCode(recvCodes[i], boundIt->rankObj);
808
 	}
809
810
811
      }
    }

812
    return meshChanged;
813
  }
814
815


816
  void MeshDistributor::serialize(ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
817
  {    
818
    int vecSize = data.size();
819
    SerUtil::serialize(out, vecSize);
820
    for (int i = 0; i < vecSize; i++) {
821
      int dofIndex = *(data[i]);
822
      SerUtil::serialize(out, dofIndex);
823
824
825
826
    }
  }


827
828
  void MeshDistributor::deserialize(istream &in, DofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
829
  {
830
    FUNCNAME("MeshDistributor::deserialize()");
831
832

    int vecSize = 0;
833
    SerUtil::deserialize(in, vecSize);
834
    data.clear();
835
836
837
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
838
      SerUtil::deserialize(in, dofIndex);
839
840
841
842
843
844
845
846
847

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

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


848
  void MeshDistributor::serialize(ostream &out, RankToDofContainer &data)
849
850
  {
    int mapSize = data.size();
851
    SerUtil::serialize(out, mapSize);
852
853
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
854
      SerUtil::serialize(out, rank);
855
856
857
858
859
      serialize(out, it->second);
    }
  }

  
860
861
  void MeshDistributor::deserialize(istream &in, RankToDofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
862
  {
863
864
    data.clear();

865
    int mapSize = 0;
866
    SerUtil::deserialize(in, mapSize);
867
868
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
869
      SerUtil::deserialize(in, rank);
870
871
872
873
      deserialize(in, data[rank], dofMap);      
    }
  }

874

875
  void MeshDistributor::setInitialElementWeights() 
876
  {
877
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
878
879

    elemWeights.clear();
880
      
881
    string filename = "";
882
883
884
885
    GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

886
887
      ifstream infile;
      infile.open(filename.c_str(), ifstream::in);
888
889
890
891
892
893
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;
894

895
896
	elemWeights[elNum] = elWeight;
      }
897

898
      infile.close();
899
    } else {           
900
      TraverseStack stack;