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


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

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

50
51
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
52
  using boost::lexical_cast;
53
  using namespace boost::filesystem;
54
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
55

56
57
58
59
60
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

61

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

84
85
86
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
87
88
89
90

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

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

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

99
100
101
102
103
104
105
106
107
108
    string partStr = "parmetis";
    GET_PARAMETER(0, name + "->partitioner", &partStr);

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
111

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

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

119
120
    TEST_EXIT(feSpace)("No FE space has been defined for the mesh distributor!\n");
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
121

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

127
128
    elObjects.setFeSpace(feSpace);

129
130
131
    // If the problem has been already read from a file, we need only to set
    // isRankDofs to all matrices and rhs vector and to remove periodic 
    // boundary conditions (if there are some).
132
    if (deserialized) {
133
134
      updateMacroElementInfo();

135
      setRankDofs();
136

137
      removePeriodicBoundaryConditions();
138

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

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

148
      return;
149
    }
150

151

152
153
154
155
156
    // Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
    // only for macro meshes, so it will not work yet if the mesh is already refined
    // in some way.
    testForMacroMesh();

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

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

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

    // and now partition the mesh    
168
    oldPartitionMap = partitioner->getPartitionMap();
169
170
171
172

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

173
    partitionMap = partitioner->getPartitionMap();
174

175

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

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

190

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

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

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

199
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
    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;
	}
      }
    }
226

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

    removeMacroElements();
230

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
231

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

Thomas Witkowski's avatar
Thomas Witkowski committed
234

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

242
    updateLocalGlobalNumbering();
243

244

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

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

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

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

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

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

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

269
    createPeriodicMap();
270

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

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

283
      updateLocalGlobalNumbering();
284
285

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

288
      createPeriodicMap();
289
290
291
292

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

295

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

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

300

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

303
    removePeriodicBoundaryConditions();
304
305
  }

306

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

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

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

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

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

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

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

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

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

    probStat.push_back(probVec);
  }


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

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

427

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

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

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

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

448
    stdMpi.startCommunication();
449

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


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

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

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

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

481
    stdMpi.startCommunication();
482
483

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

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


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


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

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

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


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


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

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

649

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

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

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

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

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

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

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

696

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

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

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

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

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

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

717

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

720
721
    updateLocalGlobalNumbering();

722

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

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

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


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

734
    nMeshChangesAfterLastRepartitioning++;
735

736
737
738
739

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

740
741
742
743
    if (repartitioningAllowed && 
	nMeshChangesAfterLastRepartitioning >= repartitionIthChange) {
      repartitionMesh();
      nMeshChangesAfterLastRepartitioning = 0;
744
    }
745
746
747
  }

  
748
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
749
  {
750
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
751
752
753

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

758
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
759
760
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
761
	elCode.init(boundIt->rankObj);
762
763
764
765
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
766
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
767
    stdMpi.send(sendCodes);
768
769
770
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it)
      stdMpi.recv(it->first);
    stdMpi.startCommunication();
771
 
772
    // === Compare received mesh structure codes. ===
773
    
774
    bool meshChanged = false;
775

776
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
777
     
778
779
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
780
      
781
      for (vector<AtomicBoundary>::iterator boundIt = it->second.begin();
782
	   boundIt != it->second.end(); ++boundIt, i++) {
783

784
785
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
786

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

790
791
	  MeshManipulation mm(feSpace);
	  meshChanged |= mm.fitElementToMeshCode(recvCodes[i], boundIt->rankObj);
792
 	}
793
794
795
      }
    }

796
    return meshChanged;
797
  }
798
799


800
  void MeshDistributor::serialize(ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
801
  {    
802
    int vecSize = data.size();
803
    SerUtil::serialize(out, vecSize);
804
    for (int i = 0; i < vecSize; i++) {
805
      int dofIndex = *(data[i]);
806
      SerUtil::serialize(out, dofIndex);
807
808
809
810
    }
  }


811
812
  void MeshDistributor::deserialize(istream &in, DofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
813
  {
814
    FUNCNAME("MeshDistributor::deserialize()");
815
816

    int vecSize = 0;
817
    SerUtil::deserialize(in, vecSize);
818
    data.clear();
819
820
821
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
822
      SerUtil::deserialize(in, dofIndex);
823
824
825
826
827
828
829
830
831

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

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


832
  void MeshDistributor::serialize(ostream &out, RankToDofContainer &data)
833
834
  {
    int mapSize = data.size();
835
    SerUtil::serialize(out, mapSize);
836
837
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
838
      SerUtil::serialize(out, rank);
839
840
841
842
843
      serialize(out, it->second);
    }
  }

  
844
845
  void MeshDistributor::deserialize(istream &in, RankToDofContainer &data,
				    map<int, const DegreeOfFreedom*> &dofMap)
846
  {
847
848
    data.clear();

849
    int mapSize = 0;
850
    SerUtil::deserialize(in, mapSize);
851
852
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
853
      SerUtil::deserialize(in, rank);
854
855
856
857
      deserialize(in, data[rank], dofMap);      
    }
  }

858

859
  void MeshDistributor::setInitialElementWeights() 
860
  {
861
    FUNCNAME("MeshDistributor::setInitialElementWeights()");
862
863

    elemWeights.clear();
864
      
865
    string filename = "";
866
867
868
869
    GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

870
871
      ifstream infile;
      infile.open(filename.c_str(), ifstream::in);
872
873
874
875
876
877
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;
878

879
880
	elemWeights[elNum] = elWeight;
      }
881

882
      infile.close();
883
    } else {           
884
      TraverseStack stack;
885
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
886
      while (elInfo) {
887
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
888
	elInfo = stack.traverseNext(elInfo);