MeshDistributor.cc 63 KB
Newer Older
1
//
2
3
4
5
6
7
8
9
10
11
12
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include <algorithm>
14
15
#include <iostream>
#include <fstream>
16
17
#include <limits>
#include <stdint.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
18
#include <boost/lexical_cast.hpp>
19
20
#include <boost/filesystem.hpp>

21
#include "parallel/MeshDistributor.h"
22
#include "parallel/MeshManipulation.h"
23
#include "parallel/ParallelDebug.h"
24
#include "parallel/StdMpi.h"
25
#include "parallel/ParMetisPartitioner.h"
26
#include "parallel/MpiHelper.h"
27
28
29
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
#include "io/VtkWriter.h"
30
31
32
33
34
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
35
36
#include "DOFMatrix.h"
#include "DOFVector.h"
37
#include "SystemVector.h"
38
#include "ElementDofIterator.h"
39
40
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
41
#include "VertexVector.h"
42
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
43
44
#include "ProblemVec.h"
#include "ProblemInstat.h"
45
#include "Debug.h"
46

47
48
namespace AMDiS {

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

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

58

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
97

98
  void MeshDistributor::initParallelization()
99
  {
100
    FUNCNAME("MeshDistributor::initParallelization()");
101
102
103
104

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

105
106
    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");
107

108
109
    elObjects.setFeSpace(feSpace);

110
111
112
    // 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).
113
    if (deserialized) {
114
115
      updateMacroElementInfo();

116
      setRankDofs();
117

118
      removePeriodicBoundaryConditions();
119

120
121
      macroElIndexMap.clear();
      macroElIndexTypeMap.clear();
122

123
124
      map<int, bool>& elementInRank = partitioner->getElementInRank();
      for (vector<MacroElement*>::iterator it = allMacroElements.begin();
125
126
	   it != allMacroElements.end(); ++it) {
	elementInRank[(*it)->getIndex()] = false;
127
128
129
130
	macroElIndexMap[(*it)->getIndex()] = (*it)->getElement();
	macroElIndexTypeMap[(*it)->getIndex()] = (*it)->getElType();
      }

131
      for (deque<MacroElement*>::iterator it = mesh->getMacroElements().begin();
132
133
134
	   it != mesh->getMacroElements().end(); ++it)
	elementInRank[(*it)->getIndex()] = true;      

135
      return;
136
    }
137

138

139
140
141
142
143
    // 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();

144
145
146
147
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

148
149
    // create an initial partitioning of the mesh
    partitioner->createPartitionData();
150

151
    // set the element weights, which are 1 at the very first begin
152
    setInitialElementWeights();
153
154
155

    // and now partition the mesh    
    partitioner->fillCoarsePartitionVec(&oldPartitionVec);
156
157
158
159

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

160
161
    partitioner->fillCoarsePartitionVec(&partitionVec);

162

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

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

177

178
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
179

180
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
181

182
#if (DEBUG != 0)
183
    ParallelDebug::printBoundaryInfo(*this);
184
185
#endif

186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
    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;
	}
      }
    }
213

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
214
215
216
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
217

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
218

219
220
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
221

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

229
    updateLocalGlobalNumbering();
230

231

232
233
234
235
236
237
    // === 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();

238
239
    // === If in debug mode, make some tests. ===

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

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

249
    debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
250
251

    MSG("Debug mode tests finished!\n");
252
#endif
253

254

255
    // === Create periodic DOF mapping, if there are periodic boundaries. ===
256

257
    createPeriodicMap();
258

259

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

268
      updateLocalGlobalNumbering();
269
270

     
271
272
      // === Update periodic mapping, if there are periodic boundaries. ===     

273
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
274
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
275

276

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

279
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
280

281

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

284
    removePeriodicBoundaryConditions();
285
286
  }

287

288
289
290
291
292
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

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

      partitioner = new ParMetisPartitioner(mesh, &mpiComm);
    }

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

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

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

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

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

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

    probStat.push_back(probVec);
  }


380
  void MeshDistributor::exitParallelization()
381
  {}
382

383
  
384
  void MeshDistributor::testForMacroMesh()
385
  {
386
    FUNCNAME("MeshDistributor::testForMacroMesh()");
387
388
389
390
391
392
393

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

408

409
  void MeshDistributor::synchVector(DOFVector<double> &vec)
410
  {
411
    StdMpi<vector<double> > stdMpi(mpiComm);
412
413

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

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

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

429
    stdMpi.startCommunication();
430

Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
433
434
435
    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];
  }
436
437


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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
455
      stdMpi.send(sendIt->first, dofs);
456
457
458
    }

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

462
    stdMpi.startCommunication();
463
464

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

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
469
470
471
472
473
      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++];
474
475
476
477
      }
    }
  }

478

479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
  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");
  }


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


576
577
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
578
579
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
609
610
611
612
  }


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


624
  void MeshDistributor::checkMeshChange()
625
  {
626
    FUNCNAME("MeshDistributor::checkMeshChange()");
627

628
629
    double first = MPI::Wtime();

630

631
632
633
634
635
    // === 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);
636

637
    if (recvAllValues == 0)
638
639
      return;

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

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

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

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

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

677

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

683
      meshChanged |= checkAndAdaptBoundary(allBound);
684
685
686

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

687
      int sendValue = static_cast<int>(meshChanged);
688
689
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
690
691

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

694
#if (DEBUG != 0)
695
    debug::writeMesh(feSpace, -1, debugOutputDir + "mesh");
696
697
698
#endif

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

700
701
    updateLocalGlobalNumbering();

702
703
704

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

Thomas Witkowski's avatar
Thomas Witkowski committed
705
    createPeriodicMap();
706
707
708
709
710

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


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

713
    nMeshChangesAfterLastRepartitioning++;
714

715
716
717
718
    if (repartitioningAllowed && 
	nMeshChangesAfterLastRepartitioning >= repartitionIthChange) {
      repartitionMesh();
      nMeshChangesAfterLastRepartitioning = 0;
719
    }
720
721
722
  }

  
723
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
724
  {
725
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
726
727
728

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

733
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
734
735
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
736
	elCode.init(boundIt->rankObj);
737
738
739
740
	sendCodes[it->first].push_back(elCode);
      }
    }

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

751
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
752
     
753
754
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
755
      
756
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
757
	   boundIt != it->second.end(); ++boundIt, i++) {
758

759
760
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
761

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

765
766
	  MeshManipulation mm(feSpace);
	  meshChanged |= mm.fitElementToMeshCode(recvCodes[i], boundIt->rankObj);
767
 	}
768
769
770
      }
    }

771
    return meshChanged;
772
  }
773
774


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


786
787
  void MeshDistributor::deserialize(istream &in, DofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
788
  {
789
    FUNCNAME("MeshDistributor::deserialize()");
790
791

    int vecSize = 0;
792
    SerUtil::deserialize(in, vecSize);
793
    data.clear();
794
795
796
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
797
      SerUtil::deserialize(in, dofIndex);
798
799
800
801
802
803
804
805
806

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

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


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

  
819
820
  void MeshDistributor::deserialize(istream &in, RankToDofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
821
  {
822
823
    data.clear();

824
    int mapSize = 0;
825
    SerUtil::deserialize(in, mapSize);
826
827
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
828
      SerUtil::deserialize(in, rank);
829
830
831
832
      deserialize(in, data[rank], dofMap);      
    }
  }

833

834
  void MeshDistributor::setInitialElementWeights() 
835
  {
836
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
837
838

    elemWeights.clear();
839
      
840
    string filename = "";
841
842
843
844
    GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

845
846
      ifstream infile;
      infile.open(filename.c_str(), ifstream::in);
847
848
849
850
851
852
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;
853

854
855
	elemWeights[elNum] = elWeight;
      }
856

857
      infile.close();
858
    } else {           
859
      TraverseStack stack;
860
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
861
      while (elInfo) {
862
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
863
	elInfo = stack.traverseNext(elInfo);
864
865
866
867
      }
    }
  }

868

869
870
871
872
873
874
875
  void MeshDistributor::repartitionMesh()
  {
    FUNCNAME("MeshDistributor::repartitionMesh()");

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

876
877
878
879
880

    // === First we check if the rank with the maximum number of DOFs has at  ===
    // === least 20% more DOFs than the rank with the minimum number of DOFs. ===
    // === In this case, the mesh will be repartition.                        ===

881
882
    int repartitioning = 0;

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

    if (mpiRank == 0) {
      int nOverallDofs = 0;
889
890
      int minDofs = numeric_limits<int>::max();
      int maxDofs = numeric_limits<int>::min();
891
892
893
894
895
896
897
898
899
      for (int i = 0; i < mpiSize; i++) {
	nOverallDofs += nDofsInRank[i];
	minDofs = std::min(minDofs, nDofsInRank[i]);
	maxDofs = std::max(maxDofs, nDofsInRank[i]);
      }      
     
      MSG("Overall DOFs: %d    Min DOFs: %d    Max DOFs: %d\n", 
	  nOverallDofs, minDofs, maxDofs);

900
      if (static_cast<double>(maxDofs) / static_cast<double>(minDofs) > 1.2) 
901
902
903
904
905
906
907
908
909
910
911
	repartitioning = 1;

      mpiComm.Bcast(&repartitioning, 1, MPI_INT, 0);
    } else {
      mpiComm.Bcast(&repartitioning, 1, MPI_INT, 0);
    }


    if (repartitioning == 0)
      return;

Thomas Witkowski's avatar