MeshDistributor.cc 88.9 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20
21


Thomas Witkowski's avatar
Thomas Witkowski committed
22
#include <algorithm>
23
24
#include <iostream>
#include <fstream>
25
26
#include <limits>
#include <stdint.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
27
#include <boost/lexical_cast.hpp>
28
29
#include <boost/filesystem.hpp>

30
#include "parallel/MeshDistributor.h"
31
#include "parallel/MeshManipulation.h"
32
33
#ifndef NDEBUG
  #include "parallel/ParallelDebug.h"
34
#endif
35
#include "parallel/StdMpi.h"
36
#include "parallel/MeshPartitioner.h"
37
#include "parallel/ParMetisPartitioner.h"
38
#include "parallel/ZoltanPartitioner.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
39
#include "parallel/SimplePartitioner.h"
40
#include "parallel/CheckerPartitioner.h"
41
#include "parallel/MpiHelper.h"
42
#include "parallel/DofComm.h"
43
#include "parallel/ParallelProblemStat.h"
44
45
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
46
#include "io/MacroWriter.h"
47
#include "io/VtkWriter.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
48
#include "io/ArhReader.h"
49
#include "io/Arh2Reader.h"
50
51
52
53
54
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
55
56
#include "DOFMatrix.h"
#include "DOFVector.h"
57
#include "SystemVector.h"
58
#include "ElementDofIterator.h"
59
60
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
61
#include "VertexVector.h"
62
#include "MeshStructure.h"
63
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
64
#include "ProblemInstat.h"
65
#include "RefinementManager3d.h"
66
67
#ifndef NDEBUG
  #include "Debug.h"
68
#endif
69
#include "Timer.h"
70
#include "io/MacroReader.h"
71

72
namespace AMDiS { namespace Parallel {
73

Thomas Witkowski's avatar
Thomas Witkowski committed
74
  using boost::lexical_cast;
75
  using namespace boost::filesystem;
76
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
77

78
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;
79

80
81
82
  const Flag MeshDistributor::BOUNDARY_SUBOBJ_SORTED              = 0X01L;
  const Flag MeshDistributor::BOUNDARY_FILL_INFO_SEND_DOFS        = 0X02L;
  const Flag MeshDistributor::BOUNDARY_FILL_INFO_RECV_DOFS        = 0X04L;
Thomas Witkowski's avatar
Thomas Witkowski committed
83

84
85
86
87
88
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

89
  MeshDistributor::MeshDistributor()
90
    : problemStat(0),
91
      initialized(false),
92
      name("parallel"),
93
      macroMesh(NULL),
94
95
96
      refineManager(NULL),
      partitioner(NULL),
      initialPartitioner(NULL),
97
      deserialized(false),
98
      writeSerializationFile(false),
99
      repartitioningAllowed(false),
100
      repartitionOnlyOnce(false),
101
      repartitionIthChange(20),
102
      repartitioningWaitAfterFail(20),
103
104
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
105
      repartitioningFailed(0),
106
      debugOutputDir(""),
107
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
108
      boundaryDofInfo(1),
109
      meshAdaptivity(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
110
      hasPeriodicBoundary(false),
111
112
      printTimings(false),
      printMemoryUsage(false)
113
  {
114
    FUNCNAME("MeshDistributor::MeshDistributor()");
Thomas Witkowski's avatar
Thomas Witkowski committed
115

116
    MPI::Intracomm &mpiComm = MPI::COMM_WORLD;
117
    mpiRank = mpiComm.Get_rank();
118

119
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
120
    Parameters::get(name + "->debug output dir", debugOutputDir);
121
    Parameters::get(name + "->repartition only once", repartitionOnlyOnce);
122
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
123
    Parameters::get(name + "->repartition wait after fail", repartitioningWaitAfterFail);
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
124
    Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
125

126
    nMeshChangesAfterLastRepartitioning = repartitionIthChange - 1;
127

128
129
    // === Create partitioner object. ===

130
    string partStr = "parmetis";
131
    Parameters::get(name + "->partitioner", partStr);
132

133
    if (partStr == "parmetis")
134
      partitioner = new ParMetisPartitioner("parallel->partitioner", &mpiComm);
135

Thomas Witkowski's avatar
Thomas Witkowski committed
136
137
    if (partStr == "zoltan") {
#ifdef HAVE_ZOLTAN
138
      partitioner = new ZoltanPartitioner("parallel->partitioner", &mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
139
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
140
      ERROR_EXIT("AMDiS was compiled without Zoltan support. Therefore you cannot make use of it!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
#endif
    }
143

144
    if (partStr == "checker")
145
      partitioner = new CheckerPartitioner("parallel->partitioner", &mpiComm);
146

Thomas Witkowski's avatar
Thomas Witkowski committed
147
    if (partStr == "simple")
148
149
      partitioner = new SimplePartitioner("parallel->partitioner", &mpiComm);

150
151
152
    if (!partitioner) {
      ERROR_EXIT("Unknown partitioner or no partitioner specified!\n");
    }
153
154
155
156
157
158
159

    // === Create initial partitioner object. ===

    partStr = "";
    Parameters::get(name + "->initial partitioner", partStr);
    if (partStr == "") {
      initialPartitioner = partitioner;
160
    } else {
161
      if (partStr == "checker") {
162
	initialPartitioner =
163
164
165
166
167
168
169
170
	  new CheckerPartitioner("parallel->initial partitioner", &mpiComm);
      } else {
	ERROR_EXIT("Not yet supported, but very easy to implement!\n");
      }
    }


    // === And read some more parameters. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
171

172
    int tmp = 0;
173
    Parameters::get(name + "->box partitioning", tmp);
Thomas Witkowski's avatar
Thomas Witkowski committed
174
    partitioner->setBoxPartitioning(static_cast<bool>(tmp));
175
    initialPartitioner->setBoxPartitioning(static_cast<bool>(tmp));
Thomas Witkowski's avatar
Thomas Witkowski committed
176

Thomas Witkowski's avatar
Thomas Witkowski committed
177
    Parameters::get(name + "->print timings", printTimings);
178
    Parameters::get(name + "->print memory usage", printMemoryUsage);
Thomas Witkowski's avatar
Thomas Witkowski committed
179

180
181
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
182
183
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
184

Thomas Witkowski's avatar
Thomas Witkowski committed
185
186
  MeshDistributor::~MeshDistributor()
  {
187
    if (partitioner) {
Thomas Witkowski's avatar
Thomas Witkowski committed
188
      delete partitioner;
189
      partitioner = NULL;
190
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
191
  }
192

193
194
195
196
197
  void MeshDistributor::addInterchangeVector(SystemVector *vec)
  {
    for (int i = 0; i < vec->getSize(); i++)
      interchangeVectors.push_back(vec->getDOFVector(i));
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
198

199
200
201
202
203
  void MeshDistributor::removeInterchangeVector(SystemVector* vec)
  {
    for (int i = 0; i < vec->getSize(); i++)
      removeInterchangeVector(vec->getDOFVector(i));
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
204

205
  void MeshDistributor::initParallelization()
206
  {
207
    FUNCNAME("MeshDistributor::initParallelization()");
208

209
210
211
    if (initialized)
      return;

Thomas Witkowski's avatar
Thomas Witkowski committed
212
    double first = MPI::Wtime();
213
    MSG("Initialization phase 1 needed %.5f seconds\n",
214
	first - ParallelProblemStat::initTimeStamp);
Thomas Witkowski's avatar
Thomas Witkowski committed
215

216
    TEST_EXIT(MPI::COMM_WORLD.Get_size() > 1)
217
      ("Parallelization does not work with only one process!\n");
218
    TEST_EXIT(feSpaces.size() > 0)
219
      ("No FE space has been defined for the mesh distributor!\n");
220
    TEST_EXIT(meshes.size() > 0)("No mesh has been defined for the mesh distributor!\n");
221

222

223
224
225
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

226
    //TODO Maybe unnecessary anymore
227
228
    bool doNext = false;
    do {
229
      doNext = false;
230
231
232
233
234
235
      for (unsigned int i = 0; i < feSpaces.size() - 1; i++) {
	if (feSpaces[i]->getBasisFcts()->getDegree() >
	    feSpaces[i + 1]->getBasisFcts()->getDegree()) {
	  const FiniteElemSpace *tmp = feSpaces[i + 1];
	  feSpaces[i + 1] = feSpaces[i];
	  feSpaces[i] = tmp;
236
237
	  doNext = true;
	}
238
      }
239
240
    } while (doNext);

241
242
243
244
    // Sort FE spaces seperately for each mesh.
    map<Mesh*, vector<const FiniteElemSpace*> >::iterator iter = meshToFeSpaces.begin();
    while(iter != meshToFeSpaces.end()) {
      vector<const FiniteElemSpace*>& spaces = iter->second;
245

246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
      bool doNext = false;
      do {
	doNext = false;
	for (size_t i = 0; i < spaces.size() - 1; i++) {
	  if (spaces[i]->getBasisFcts()->getDegree() >
	      spaces[i + 1]->getBasisFcts()->getDegree()) {
	    const FiniteElemSpace *tmp = spaces[i + 1];
	    spaces[i + 1] = spaces[i];
	    spaces[i] = tmp;
	    doNext = true;
	  }
	}
      } while (doNext);
      ++iter;
    }
261

262
263
264
265
    // 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();
266

267
#ifndef NDEBUG
268
    // Check whether meshes come from the same macro mesh. The way is to compare
269
270
271
    // the node coords of each macro element in the meshes.
    debug::ElementIdxToCoords macroCoords;
    debug::createNodeCoords(macroMesh, macroCoords);
272

273
    for (size_t i = 1; i < meshes.size(); i++)
274
      debug::testNodeCoords(meshes[i], macroCoords);
275

276
#endif
277

278
    // Initialize dof communicators which have been created in addProblemStat.
279
    for (size_t i = 0; i < meshes.size(); i++)
280
      dofComms[meshes[i]].init(levelData, meshToFeSpaces[meshes[i]]);
Thomas Witkowski's avatar
Thomas Witkowski committed
281

282

283
284
    // Initialize element object DB
    elObjDb.init(meshes, meshToFeSpaces);
285

286
    // If the problem has been already read from a file, we need only to set
287
    // isRankDofs to all matrices and rhs vector and to remove periodic
288
    // boundary conditions (if there are some).
289
    if (deserialized) {
290
291
      updateMacroElementInfo();

292
      removePeriodicBoundaryConditions();
293

294
      elObjDb.createMacroElementInfo(allMacroElements);
295

296
      updateDofRelatedStruct();
297

Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
      elObjDb.setData(partitionMap, levelData);

300
#ifndef NDEBUG
301
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
302
303
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1],
				    *(dofMaps[0]),
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
304
				    debugOutputDir + "mpi-dbg", "dat");
305
#endif
306
307

      initialized = true;
308
      return;
309
    }
310

311

Thomas Witkowski's avatar
Thomas Witkowski committed
312
313
    // Create a very first pariitioning of the mesh.
    createInitialPartitioning();
314

315

316
#ifndef NDEBUG
317
318
319
320
    std::vector<debug::ElementIdxToDofs> elMap(meshes.size());
    for (size_t i = 0; i < meshes.size(); i++) {
      debug::createSortedDofs(meshes[i], elMap[i]);
    }
321

322
323
    if (mpiRank == 0) {
      int writePartMesh = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
324
      Parameters::get("parallel->debug->write mesh partitioning", writePartMesh);
325

326
      if (writePartMesh > 0) {
327
	debug::writeElementIndexMesh(macroMesh , debugOutputDir + "elementIndex");
328
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part");
329
      }
330
    }
331
#endif
332

333
    // Create interior boundary information.
334
    createInteriorBoundary(true);
335

336
    // === Remove neighbourhood relations due to periodic bounday conditions. ===
337

338
339
340
341
342
343
    for (size_t i = 0; i < meshes.size(); i++) {
      for (deque<MacroElement*>::iterator it = meshes[i]->firstMacroElement();
	 it != meshes[i]->endOfMacroElements(); ++it) {
	for (int j = 0; j < macroMesh->getGeo(NEIGH); j++) {
	  if ((*it)->getNeighbour(j) &&
	    meshes[i]->isPeriodicAssociation((*it)->getBoundary(j))) {
344

345
	    int neighIndex = (*it)->getNeighbour(j)->getIndex();
346

347
348
349
	    (*it)->getNeighbour(j)->setNeighbour((*it)->getOppVertex(j), NULL);
	    (*it)->setNeighbour(j, NULL);
	    (*it)->setBoundary(j, 0);
350

351
352
353
354
355
	    if (i == 0) {
	      macroElementNeighbours[(*it)->getIndex()][j] = -1;
	      macroElementNeighbours[neighIndex][(*it)->getOppVertex(j)] = -1;
	    }
	  }
356
357
358
	}
      }
    }
359

360
361
362
363
    for (map<Mesh*, std::vector<MacroElement*> >::iterator iter = allMacroElements.begin();
	iter != allMacroElements.end(); iter++) {
      Mesh* mesh = iter->first;
      vector<MacroElement*>& allMacros = iter->second;
364

365
366
367
368
      for (vector<MacroElement*>::iterator it = allMacros.begin();
	 it != allMacros.end(); ++it) {
	for (int i = 0; i < macroMesh->getGeo(NEIGH); i++) {
	  if ((*it)->getNeighbour(i) &&
369
	    mesh->isPeriodicAssociation((*it)->getBoundary(i))) {
370

371
	    int neighIndex = (*it)->getNeighbour(i)->getIndex();
372

373
374
375
	    (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	    (*it)->setNeighbour(i, NULL);
	    (*it)->setBoundary(i, 0);
376

377
378
379
380
381
	    if (i == 0) {
	      macroElementNeighbours[(*it)->getIndex()][i] = -1;
	      macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
	    }
	  }
382
383
384
	}
      }
    }
385

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

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
388
    removeMacroElements();
389

390
    // === Create new global and local DOF numbering. ===
391
392

    // We have to remove the VertexVectors, which contain periodic assoiciations,
393
394
    // because they are not valid anymore after some macro elements have been
    // removed and the corresponding DOFs were deleted.
395
396
397
398
    for (size_t i = 0; i < meshes.size(); i++)
      for (map<BoundaryType, VertexVector*>::iterator it = meshes[i]->getPeriodicAssociations().begin();
	  it != meshes[i]->getPeriodicAssociations().end(); ++it)
	const_cast<DOFAdmin&>(meshes[i]->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
399

400

401

402
    updateDofRelatedStruct();
403

404
405
406
    // === In 3D we have to fix the mesh to allow local refinements. ===

    fix3dMeshRefinement();
407

408
409
    // === If in debug mode, make some tests. ===

410
#ifndef NDEBUG
411
    MSG("AMDiS runs in debug mode, so make some test ...\n");
412

413
    ParallelDebug::testAllElements(*this);
414
415
    for (size_t i = 0; i < meshes.size(); i++)
      debug::testSortedDofs(meshes[i], elMap[i]);
416
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
417
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
418

419
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");
420
421

    MSG("Debug mode tests finished!\n");
422
#endif
423

424
    // Remove periodic boundary conditions in sequential problem definition.
425
426
    removePeriodicBoundaryConditions();

427
#ifndef NDEBUG
428
429
    ParallelDebug::testPeriodicBoundary(*this);
#endif
430

431
    // === Global refinements. ===
432
433
    bool oneMeshRefined = false;
    for (size_t i = 0; i < meshes.size(); i++) {
434

435
436
      int globalRefinement = 0;
      Parameters::get(meshes[i]->getName() + "->global refinements", globalRefinement);
437

438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
      if (globalRefinement > 0) {
        oneMeshRefined = true;
	bool doRefineInter = true;
	Parameters::get(meshes[i]->getName() + "->refinement interpol", doRefineInter);
	std::map<DOFVector<double>*, RefineCoarsenOperation> rememberOp;
	if (!doRefineInter) {
	  // No refinement during initial global refinement
	  for (int iadmin = 0; iadmin < meshes[i]->getNumberOfDOFAdmin(); iadmin++) {
	    std::list<DOFIndexedBase*>::iterator it;
	    DOFAdmin* admin = const_cast<DOFAdmin*>(&meshes[i]->getDofAdmin(iadmin));
	    std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
	    for (it = admin->beginDOFIndexed(); it != end; it++) {
	      DOFVector<double>* vec = dynamic_cast<DOFVector<double>*>(*it);
	      if (vec) {
		rememberOp[vec] = vec->getRefineOperation();
		vec->setRefineOperation(NO_OPERATION);
	      }
455
456
457
	    }
	  }
	}
458
459
460
461
462
463
464
465
466
467
468
469
	refineManager->globalRefine(meshes[i], globalRefinement);
	if (!doRefineInter) {
	  // No refinement during initial global refinement
	  for (int iadmin = 0; iadmin < meshes[i]->getNumberOfDOFAdmin(); iadmin++) {
	    std::list<DOFIndexedBase*>::iterator it;
	    DOFAdmin* admin = const_cast<DOFAdmin*>(&meshes[i]->getDofAdmin(iadmin));
	    std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
	    for (it = admin->beginDOFIndexed(); it != end; it++) {
	      DOFVector<double>* vec = dynamic_cast<DOFVector<double>*>(*it);
	      if (vec)
		vec->setRefineOperation(rememberOp[vec]);
	    }
470
471
	  }
	}
472
	updateDofRelatedStruct(meshes[i]);
473

474
#ifndef NDEBUG
475
      ParallelDebug::testPeriodicBoundary(*this);
476
#endif
477
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
478
    }
479
480
    if (oneMeshRefined)
      updateLocalGlobalNumbering();
481

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
482
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
483
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
484
485
      elObjDb.clear();

486
    initialized = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
487
    MSG("Init parallelization needed %.5f seconds\n", MPI::Wtime() - first);
488
489
  }

490

Thomas Witkowski's avatar
Thomas Witkowski committed
491
492
493
494
495
496
497
498
  void MeshDistributor::createInitialPartitioning()
  {
    FUNCNAME("MeshDistributor::createInitialPartitioning()");

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

499
    // Create an initial partitioning of the mesh
500
    bool useInitialPartitioning =
501
      initialPartitioner->createInitialPartitioning();
Thomas Witkowski's avatar
Thomas Witkowski committed
502

503
    // Set the element weights, which are 1 at the very first begin
Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
    setInitialElementWeights();

506
507
508
    if (!useInitialPartitioning) {
      // And now partition the mesh
      bool partitioningSucceed =
509
	initialPartitioner->partition(elemWeights, INITIAL);
Thomas Witkowski's avatar
Thomas Witkowski committed
510
511
      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
    }
512

513
514
515
516
517
    initialPartitioner->createPartitionMap(partitionMap);

    if (initialPartitioner != partitioner) {
      *partitioner = *initialPartitioner;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
518
519
520
  }


521
  void MeshDistributor::setInitialElementWeights()
Thomas Witkowski's avatar
Thomas Witkowski committed
522
523
524
525
  {
    FUNCNAME("MeshDistributor::setInitialElementWeights()");

    elemWeights.clear();
526

Thomas Witkowski's avatar
Thomas Witkowski committed
527
    string filename = "";
528
    Parameters::get(macroMesh->getName() + "->macro weights", filename);
Thomas Witkowski's avatar
Thomas Witkowski committed
529
530
531
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

532
533
      std::ifstream infile;
      infile.open(filename.c_str(), std::ifstream::in);
Thomas Witkowski's avatar
Thomas Witkowski committed
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;

	elemWeights[elNum] = elWeight;
      }

      infile.close();

      return;
    }

    filename = "";
    Parameters::get("parallel->partitioner->read meta arh", filename);
    if (filename != "") {
      map<int, int> arhElInRank;
      map<int, int> arhElCodeSize;

555
      /*int nProc = */ io::Arh2Reader::readMetaData(filename, arhElInRank, arhElCodeSize);
556
      for (map<int, int>::iterator it = arhElCodeSize.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
557
558
559
560
561
562
563
	   it != arhElCodeSize.end(); ++it)
	elemWeights[it->first] = it->second;

      return;
    }

    TraverseStack stack;
564
    ElInfo *elInfo = stack.traverseFirst(macroMesh, -1, Mesh::CALL_LEAF_EL);
Thomas Witkowski's avatar
Thomas Witkowski committed
565
566
    while (elInfo) {
      elemWeights[elInfo->getElement()->getIndex()] = 1.0;
567
      elInfo = stack.traverseNext(elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
568
569
570
571
    }
  }


572
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
573
  {
574
    FUNCNAME("MeshDistributor::addProblemStat()");
575

576
577
578
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

579

580
    // === Add all FE spaces from stationary problem. ===
581

582
    vector<const FiniteElemSpace*> newFeSpaces = probStat->getFeSpaces();
583

584
    for (size_t i = 0; i < newFeSpaces.size(); i++)
585
      if (find(feSpaces.begin(), feSpaces.end(), newFeSpaces[i]) ==
586
	  feSpaces.end())
587
588
589
590
      {
	// Find new fespace, then add mesh if it's new
	Mesh* mesh = newFeSpaces[i]->getMesh();
	if (find(meshes.begin(), meshes.end(), mesh) == meshes.end()) {
591
592
593
	  TEST_EXIT(meshes.size() < 2)
	    ("Currently max two meshes supported on parallel mode.\n");

594
595
596
597
	  meshes.push_back(mesh);
	  dofComms[mesh] = MultiLevelDofComm();
	  lastMeshChangeIndexs[mesh] = 0;
	}
598

599
	feSpaces.push_back(newFeSpaces[i]);
600
601
	meshToFeSpaces[mesh].push_back(newFeSpaces[i]);
      }
602

603
604
    // === Set the first mesh to be the macro mesh and pass it to partitioner.    ===
    // === Create a corresponding refinement manager object.                      ===
605

606
607
    if (problemStat.empty()) {
      macroMesh = meshes[0];
608

609
610
611
612
613
614
615
616
617
618
      switch (macroMesh->getDim()) {
      case 2:
	refineManager = new RefinementManager2d();
	break;
      case 3:
	refineManager = new RefinementManager3d();
	break;
      default:
	ERROR_EXIT("This should not happen for dim = %d!\n", macroMesh->getDim());
      }
619

620
621
      partitioner->setMesh(macroMesh);
      initialPartitioner->setMesh(macroMesh);
622
    }
623

624

625
626
627
    // === Check whether the stationary problem should be serialized. ===


628
629
    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
630
    Parameters::get(probStat->getName() + "->output->write serialization",
631
		    writeSerialization);
632
    if (writeSerialization && !writeSerializationFile) {
633
      string filename = "";
634
      Parameters::get(name + "->output->serialization filename", filename);
635

636
637
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
638
639

      int tsModulo = -1;
640
      Parameters::get(probStat->getName() + "->output->write every i-th timestep",
641
		      tsModulo);
642

643
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
644
      writeSerializationFile = true;
645
    }
646

647
648
649

    // === Check whether the stationary problem should be deserialized. ===

650
    int readSerialization = 0;
651
    Parameters::get(probStat->getName() + "->input->read serialization",
652
		    readSerialization);
653
    if (readSerialization) {
654
      string filename = "";
655
      Parameters::get(probStat->getName() + "->input->serialization filename",
656
		      filename);
657
      filename += ".p" + lexical_cast<string>(mpiRank);
658
      MSG("Start deserialization with %s\n", filename.c_str());
659
      std::ifstream in(filename.c_str());
660
661
662
663

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

664
      // Read the revision number of the AMDiS version which was used to create
665
666
667
668
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

669
      probStat->deserialize(in);
670
      in.close();
671
672
      MSG("Deserialization from file: %s\n", filename.c_str());

673
      filename = "";
674
      Parameters::get(name + "->input->serialization filename", filename);
675

676
677
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
678

679
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
680
      in.open(rankFilename.c_str());
681

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

685
      // Read the revision number of the AMDiS version which was used to create
686
687
688
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
689

690
691
692
693
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
694
695
    }

696
    problemStat.push_back(probStat);
697

698
699
700

    // === If the mesh distributor is already initialized, don't forget to  ===
    // === remove the periodic boundary conditions on these objects.        ===
701

702
    if (initialized)
703
      removePeriodicBoundaryConditions(probStat);
704
705
706
  }


707
708
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
709
    if (globalMeshDistributor == NULL)
710
711
712
713
714
715
      globalMeshDistributor = new MeshDistributor();

    globalMeshDistributor->addProblemStat(probStat);
  }


716
  void MeshDistributor::exitParallelization()
717
  {}
718

719
720
721
722
723
724
725
726
727

  void MeshDistributor::registerDofMap(ParallelDofMapping &dofMap)
  {
    FUNCNAME("MeshDistributor::registerDofMap()");

    TEST_EXIT(find(dofMaps.begin(), dofMaps.end(), &dofMap) ==
	      dofMaps.end())
      ("Parallel DOF mapping already registerd in mesh distributor object!\n");

728
    dofMaps.push_back(&dofMap);
729
730
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
731
732
733

  void MeshDistributor::removeDofMap(ParallelDofMapping &dofMap)
  {
734
    vector<ParallelDofMapping*>::iterator it =
Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
      find(dofMaps.begin(), dofMaps.end(), &dofMap);

737
738
    if (it != dofMaps.end())
      dofMaps.erase(it);
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
741
  }


742
  void MeshDistributor::testForMacroMesh()
743
  {
744
    FUNCNAME("MeshDistributor::testForMacroMesh()");
745

746
    for (size_t i = 0; i < meshes.size(); i++) {
747

748
      int nMacroElements = 0;
749

750
751
752
753
754
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(meshes[i], -1, Mesh::CALL_LEAF_EL);
      while (elInfo) {
	TEST_EXIT(elInfo->getLevel() == 0)
	  ("Mesh is already refined! This does not work with parallelization!\n");
755

756
757
	TEST_EXIT(elInfo->getType() == 0)
	  ("Only macro elements with level 0 are supported!\n");
758

759
	nMacroElements++;
760

761
762
	elInfo = stack.traverseNext(elInfo);
      }
763

764
765
766
      TEST_EXIT(nMacroElements >= MPI::COMM_WORLD.Get_size())
	("The mesh has less macro elements than number of mpi processes!\n");
    }
767
768
  }

769

770
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
771
  {
772
    FUNCNAME_DBG("MeshDistributor::synchVector()");
773

774
775
776
    int nLevels = levelData.getNumberOfLevels();
    for (int level = nLevels - 1; level >= 0; level--) {
      MPI::Intracomm &mpiComm = levelData.getMpiComm(level);
777

778
779
      if (mpiComm == MPI::COMM_SELF)
	continue;
780

781
782
783
784
785
786
787
788
789
790
      vector<Mesh*> sysMeshes;

      // std::set of mesh pointer cannot be used for parallel
      for (size_t i = 0; i < meshes.size(); i++)
        for (int j = 0; j < vec.getSize(); j++)
          if (meshes[i] == vec.getFeSpace(j)->getMesh()) {
            sysMeshes.push_back(meshes[i]);
            break;
          }

791
      TEST_EXIT_DBG(sysMeshes.size() <= 2)
792
        ("This really should not happen.\n");
793

794
795
796
797
      vector<Mesh*>::iterator meshIt = sysMeshes.begin();
      while(meshIt != sysMeshes.end())
      {
	StdMpi<vector<double> > stdMpi(mpiComm);
798