MeshDistributor.cc 65.1 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
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
    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;
	}
      }
    }
229

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
230
231
232
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
233

234
235
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
236

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

244
    updateLocalGlobalNumbering();
245

246
247
248
249
250
251
    // === 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();

252
253
    // === If in debug mode, make some tests. ===

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

257
    ParallelDebug::testAllElements(*this);
258
    debug::testSortedDofs(mesh, elMap);
259
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
260

261
    debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
262
263

    MSG("Debug mode tests finished!\n");
264
#endif
265

266
    // === Create periodic DOF mapping, if there are periodic boundaries. ===
267

268
    createPeriodicMap();
269

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

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

282
      updateLocalGlobalNumbering();
283
284

     
285
286
      // === Update periodic mapping, if there are periodic boundaries. ===     

287
      createPeriodicMap();
288
289
290
291

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

294

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

297
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
298

299

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

302
    removePeriodicBoundaryConditions();
303
304
  }

305

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

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

337
      partitioner->setMesh(mesh);
338
339
340
341
    }

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

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

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

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

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

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

    probStat.push_back(probVec);
  }


398
  void MeshDistributor::exitParallelization()
399
  {}
400

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

426

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

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

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

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

447
    stdMpi.startCommunication();
448

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


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

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

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

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

480
    stdMpi.startCommunication();
481
482

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

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
487
488
489
490
491
      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++];
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
573
574
  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");
  }


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


594
595
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
596
597
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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

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


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


642
  void MeshDistributor::checkMeshChange()
643
  {
644
    FUNCNAME("MeshDistributor::checkMeshChange()");
645

646
647
    double first = MPI::Wtime();

648

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

655
    if (recvAllValues == 0)
656
657
      return;

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

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

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

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

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

695

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

701
      meshChanged |= checkAndAdaptBoundary(allBound);
702
703
704

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

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

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

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

716

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

719
720
    updateLocalGlobalNumbering();

721

722
    // === Update periodic mapping, if there are periodic boundaries. ===
723

Thomas Witkowski's avatar
Thomas Witkowski committed
724
    createPeriodicMap();
725

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


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

733
    nMeshChangesAfterLastRepartitioning++;
734

735
736
737
738

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

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

    // === 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);
    }
763
764
765
  }

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

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

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

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

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

802
803
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
804

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

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

814
    return meshChanged;
815
  }
816
817


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


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

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

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

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


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

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

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

876

877
  void MeshDistributor::setInitialElementWeights() 
878
  {
879
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
880
881

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

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

897
898
	elemWeights[elNum] = elWeight;
      }
899

900
      infile.close();
901
    } else {           
902
      TraverseStack