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

200
201
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
    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;
	}
      }
    }
227

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

    removeMacroElements();
231

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
232

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

Thomas Witkowski's avatar
Thomas Witkowski committed
235

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

243
    updateLocalGlobalNumbering();
244

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
260
    ParallelDebug::testInteriorBoundary(*this);
    ParallelDebug::testCommonDofs(*this, true);
261
    ParallelDebug::testGlobalIndexByCoords(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
262

263
    debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
264
265

    MSG("Debug mode tests finished!\n");
266
#endif
267

268
    // === Create periodic DOF mapping, if there are periodic boundaries. ===
269

270
    createPeriodicMap();
271

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

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

284
      updateLocalGlobalNumbering();
285
286

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

289
      createPeriodicMap();
290
291
292
293

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

296

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

299
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
300

301

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

304
    removePeriodicBoundaryConditions();
305
306
  }

307

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

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

339
      partitioner->setMesh(mesh);
340
341
342
343
    }

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

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

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

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

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

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

    probStat.push_back(probVec);
  }


400
  void MeshDistributor::exitParallelization()
401
  {}
402

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

428

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

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

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

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

449
    stdMpi.startCommunication();
450

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


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

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

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

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

482
    stdMpi.startCommunication();
483
484

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

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
489
490
491
492
493
      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++];
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
575
576
  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");
  }


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


596
597
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
598
599
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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

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


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


644
  void MeshDistributor::checkMeshChange()
645
  {
646
    FUNCNAME("MeshDistributor::checkMeshChange()");
647

648
649
    double first = MPI::Wtime();

650

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

657
    if (recvAllValues == 0)
658
659
      return;

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

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

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

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

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

697

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

703
      meshChanged |= checkAndAdaptBoundary(allBound);
704
705
706

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

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

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

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

718

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

721
722
    updateLocalGlobalNumbering();

723

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

Thomas Witkowski's avatar
Thomas Witkowski committed
726
    createPeriodicMap();
727

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


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

735
    nMeshChangesAfterLastRepartitioning++;
736

737
738
739
740

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

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

    // === 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);
    }
765
766
767
  }

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

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

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

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

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

804
805
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
806

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

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

816
    return meshChanged;
817
  }
818
819


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


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

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

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

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


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

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

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

878

879
  void MeshDistributor::setInitialElementWeights() 
880
  {
881
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
882
883

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

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

899
900
	elemWeights[elNum] = elWeight;
      }
901

902