MeshDistributor.cc 81.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * 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.
 * 
 ******************************************************************************/
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
#include "parallel/ParallelDebug.h"
33
#include "parallel/StdMpi.h"
34
#include "parallel/MeshPartitioner.h"
35
#include "parallel/ParMetisPartitioner.h"
36
#include "parallel/ZoltanPartitioner.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
37
#include "parallel/SimplePartitioner.h"
38
#include "parallel/CheckerPartitioner.h"
39
#include "parallel/MpiHelper.h"
40
#include "parallel/DofComm.h"
41
#include "parallel/ParallelProblemStat.h"
42
43
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
44
#include "io/MacroWriter.h"
45
#include "io/VtkWriter.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
46
#include "io/ArhReader.h"
47
#include "io/Arh2Reader.h"
48
49
50
51
52
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
53
54
#include "DOFMatrix.h"
#include "DOFVector.h"
55
#include "SystemVector.h"
56
#include "ElementDofIterator.h"
57
58
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
59
#include "VertexVector.h"
60
#include "MeshStructure.h"
61
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
62
#include "ProblemInstat.h"
63
#include "RefinementManager3d.h"
64
#include "Debug.h"
65
#include "Timer.h"
66
#include "io/MacroReader.h"
67

68
namespace AMDiS { namespace Parallel {
69

Thomas Witkowski's avatar
Thomas Witkowski committed
70
  using boost::lexical_cast;
71
  using namespace boost::filesystem;
72
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
73

74
  MeshDistributor* MeshDistributor::globalMeshDistributor = nullptr;
75

76
77
78
  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
79

80
81
82
83
84
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

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

113
    MPI::Intracomm &mpiComm = MPI::COMM_WORLD;
114
    mpiRank = mpiComm.Get_rank();
115

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

123
    nMeshChangesAfterLastRepartitioning = repartitionIthChange - 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
124
    
125
126
    // === Create partitioner object. ===

127
    string partStr = "parmetis";
128
    Parameters::get(name + "->partitioner", partStr);
129
130

    if (partStr == "parmetis") 
131
      partitioner = new ParMetisPartitioner("parallel->partitioner", &mpiComm);
132

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

141
    if (partStr == "checker")
142
      partitioner = new CheckerPartitioner("parallel->partitioner", &mpiComm);
143

Thomas Witkowski's avatar
Thomas Witkowski committed
144
    if (partStr == "simple")
145
146
      partitioner = new SimplePartitioner("parallel->partitioner", &mpiComm);

147
148
149
    if (!partitioner) {
      ERROR_EXIT("Unknown partitioner or no partitioner specified!\n");
    }
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167

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

    partStr = "";
    Parameters::get(name + "->initial partitioner", partStr);
    if (partStr == "") {
      initialPartitioner = partitioner;
    } else {    
      if (partStr == "checker") {
	initialPartitioner = 
	  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
168

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

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

177
178
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
179
180
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
181

Thomas Witkowski's avatar
Thomas Witkowski committed
182
183
  MeshDistributor::~MeshDistributor()
  {
184
    if (partitioner) {
Thomas Witkowski's avatar
Thomas Witkowski committed
185
      delete partitioner;
186
      partitioner = nullptr;
187
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
188
  }
189
190
191
192
193
194
  
  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
195

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

202
  void MeshDistributor::initParallelization()
203
  {
204
    FUNCNAME("MeshDistributor::initParallelization()");
205

206
207
208
    if (initialized)
      return;

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

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

219

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

    bool doNext = false;
    do {
225
      doNext = false;
226
227
228
229
230
231
      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;
232
233
	  doNext = true;
	}
234
      }
235
236
    } while (doNext);

237
    elObjDb.setFeSpace(feSpaces[0]);
238

Thomas Witkowski's avatar
Thomas Witkowski committed
239

240
241
242
    // 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).
243
    if (deserialized) {
244
245
      updateMacroElementInfo();

246
      removePeriodicBoundaryConditions();
247

248
      elObjDb.createMacroElementInfo(allMacroElements);
249

Thomas Witkowski's avatar
Thomas Witkowski committed
250
      updateLocalGlobalNumbering();
251

Thomas Witkowski's avatar
Thomas Witkowski committed
252
253
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
254
#if (DEBUG != 0)
255
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
256
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], 
257
				    *(dofMaps[0]), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
258
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
259
#endif    
260
261

      initialized = true;
262
      return;
263
    }
264

265

Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
    // Create a very first pariitioning of the mesh.
    createInitialPartitioning();
268

269

270
#if (DEBUG != 0)    
271
272
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
273
274
#endif

275
    if (mpiRank == 0) {
276
#if (DEBUG != 0)    
277
      int writePartMesh = 1;
278
279
280
#else
      int writePartMesh = 0;
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
281
      Parameters::get("parallel->debug->write mesh partitioning", writePartMesh);
282

283
      if (writePartMesh > 0) {
284
	debug::writeElementIndexMesh(mesh , debugOutputDir + "elementIndex");
285
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part");
286
      }
287
    }
288

289
    // Create interior boundary information.
290
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
291

292

293
294
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

295
296
297
298
299
    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))) {
300
301
302

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

303
304
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), nullptr);
	  (*it)->setNeighbour(i, nullptr);
305
306
307
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
308
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
309
310
311
312
313
314
315
316
317
	}
      }
    }

    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))) {
318
319
320

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

321
322
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), nullptr);
	  (*it)->setNeighbour(i, nullptr);
323
324
325
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
326
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
327

328
329
330
	}
      }
    }
331

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

    removeMacroElements();
335

336
337
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
338

339
    // We have to remove the VertexVectors, which contain periodic assoiciations, 
340
341
    // because they are not valid anymore after some macro elements have been
    // removed and the corresponding DOFs were deleted.
342
    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
343
344
	 it != mesh->getPeriodicAssociations().end(); ++it)
      const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
345

346

347
    updateLocalGlobalNumbering();
348

349

350
351
352
    // === In 3D we have to fix the mesh to allow local refinements. ===

    fix3dMeshRefinement();
353

354

355
356
    // === If in debug mode, make some tests. ===

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

360
    ParallelDebug::testAllElements(*this);
361
    debug::testSortedDofs(mesh, elMap);
362
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
363
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
364

365
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
366
367

    MSG("Debug mode tests finished!\n");
368
#endif
369

370
371
372
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

373
374
375
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
376

377
    // === Global refinements. ===
378

Thomas Witkowski's avatar
Thomas Witkowski committed
379
    int globalRefinement = 0;
380
    Parameters::get(mesh->getName() + "->global refinements", globalRefinement);
381
    
Thomas Witkowski's avatar
Thomas Witkowski committed
382
    if (globalRefinement > 0) {
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
      bool doRefineInter = true;
      Parameters::get(mesh->getName() + "->refinement interpol", doRefineInter);
      std::map<DOFVector<double>*, RefineCoarsenOperation> rememberOp;
      if (!doRefineInter) {
	// no refinement during initial global refinement
	for (int iadmin = 0; iadmin < mesh->getNumberOfDOFAdmin(); iadmin++) {
	  std::list<DOFIndexedBase*>::iterator it;
	  DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->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);
	    }
	  }
	}
      }
401
      refineManager->globalRefine(mesh, globalRefinement);
402
403
404
405
406
407
408
409
410
411
412
413
414
      if (!doRefineInter) {
	// no refinement during initial global refinement
	for (int iadmin = 0; iadmin < mesh->getNumberOfDOFAdmin(); iadmin++) {
	  std::list<DOFIndexedBase*>::iterator it;
	  DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->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]);
	  }
	}
      }
415
      updateLocalGlobalNumbering();
416

417
418
419
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
420
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
421

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
422
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
423
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
424
425
      elObjDb.clear();

426
    initialized = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
427
    MSG("Init parallelization needed %.5f seconds\n", MPI::Wtime() - first);
428
429
  }

430

Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
433
434
435
436
437
438
439
440
441
442
443
444
  void MeshDistributor::createInitialPartitioning()
  {
    FUNCNAME("MeshDistributor::createInitialPartitioning()");

    // 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();

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

    // create an initial partitioning of the mesh
445
446
    bool useInitialPartitioning = 
      initialPartitioner->createInitialPartitioning();
Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
449
450

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

451
    if (!useInitialPartitioning) {   
Thomas Witkowski's avatar
Thomas Witkowski committed
452
      // and now partition the mesh    
453
454
      bool partitioningSucceed = 
	initialPartitioner->partition(elemWeights, INITIAL);
Thomas Witkowski's avatar
Thomas Witkowski committed
455
456
      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
    }
457

458
459
460
461
462
    initialPartitioner->createPartitionMap(partitionMap);

    if (initialPartitioner != partitioner) {
      *partitioner = *initialPartitioner;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
  }


  void MeshDistributor::setInitialElementWeights() 
  {
    FUNCNAME("MeshDistributor::setInitialElementWeights()");

    elemWeights.clear();
      
    string filename = "";
    Parameters::get(mesh->getName() + "->macro weights", filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

      ifstream infile;
      infile.open(filename.c_str(), ifstream::in);
      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;

500
      /*int nProc = */ io::Arh2Reader::readMetaData(filename, arhElInRank, arhElCodeSize);
Thomas Witkowski's avatar
Thomas Witkowski committed
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
      for (map<int, int>::iterator it = arhElCodeSize.begin(); 
	   it != arhElCodeSize.end(); ++it)
	elemWeights[it->first] = it->second;

      return;
    }

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      elemWeights[elInfo->getElement()->getIndex()] = 1.0;
      elInfo = stack.traverseNext(elInfo);    
    }
  }


517
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
518
  {
519
    FUNCNAME("MeshDistributor::addProblemStat()");
520

521
522
523
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

524

525
526
527
528
529
530
531
532
533
534
535
    // === Add all FE spaces from stationary problem. ===

    vector<const FiniteElemSpace*> newFeSpaces = probStat->getFeSpaces();
    for (int i = 0; i < static_cast<int>(newFeSpaces.size()); i++)
      if (find(feSpaces.begin(), feSpaces.end(), newFeSpaces[i]) == 
	  feSpaces.end())
	feSpaces.push_back(newFeSpaces[i]);

    // === Add mesh of stationary problem and create a corresponding  ===
    // === refinement manager object.                                 ===
    
536
    if (mesh != nullptr) {
537
      TEST_EXIT(mesh == probStat->getMesh())
538
	("Does not yet support for different meshes! %s\n", probStat->getMesh()->getName().c_str());
539
540
541
    } else {
      mesh = probStat->getMesh();
    }
542
543
544
545
546
547
548
549
550
551
    
    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());
552
    }
553
554
    
    partitioner->setMesh(mesh);
555
    initialPartitioner->setMesh(mesh);
556
    
557

558
559
560
    // === Check whether the stationary problem should be serialized. ===


561
562
    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
563
    Parameters::get(probStat->getName() + "->output->write serialization",  
564
		    writeSerialization);
565
    if (writeSerialization && !writeSerializationFile) {
566
      string filename = "";
567
      Parameters::get(name + "->output->serialization filename", filename);
568
569
570
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
571
572

      int tsModulo = -1;
573
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
574
		      tsModulo);
575
      
576
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
577
578
      writeSerializationFile = true;
    }    
579

580
581
582

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

583
    int readSerialization = 0;
584
    Parameters::get(probStat->getName() + "->input->read serialization", 
585
		    readSerialization);
586
    if (readSerialization) {
587
      string filename = "";
588
      Parameters::get(probStat->getName() + "->input->serialization filename", 
589
		      filename);
590
      filename += ".p" + lexical_cast<string>(mpiRank);
591
      MSG("Start deserialization with %s\n", filename.c_str());
592
      ifstream in(filename.c_str());
593
594
595
596

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

597
598
599
600
601
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

602
      probStat->deserialize(in);
603
      in.close();
604
605
      MSG("Deserialization from file: %s\n", filename.c_str());

606
      filename = "";
607
      Parameters::get(name + "->input->serialization filename", filename);
608
609
610
611
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
612
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
613
614
615
616
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
617
618
619
620
621

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
622
623
624
625
626
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
627
628
    }

629
    problemStat.push_back(probStat);
630

631
632
633

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

635
    if (initialized)
636
      removePeriodicBoundaryConditions(probStat);
637
638
639
  }


640
641
642
643
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

644
    if (globalMeshDistributor == nullptr)
645
646
647
648
649
650
      globalMeshDistributor = new MeshDistributor();

    globalMeshDistributor->addProblemStat(probStat);
  }


651
  void MeshDistributor::exitParallelization()
652
  {}
653

654
655
656
657
658
659
660
661
662

  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");

663
    dofMaps.push_back(&dofMap);   
664
665
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
666
667
668
669
670
671

  void MeshDistributor::removeDofMap(ParallelDofMapping &dofMap)
  {
    vector<ParallelDofMapping*>::iterator it = 
      find(dofMaps.begin(), dofMaps.end(), &dofMap);

672
673
    if (it != dofMaps.end())
      dofMaps.erase(it);
Thomas Witkowski's avatar
Thomas Witkowski committed
674
675
676
  }


677
  void MeshDistributor::testForMacroMesh()
678
  {
679
    FUNCNAME("MeshDistributor::testForMacroMesh()");
680
681
682
683
684
685
686

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
691
692
693
694
695
696
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

697
    TEST_EXIT(nMacroElements >= MPI::COMM_WORLD.Get_size())
698
699
700
      ("The mesh has less macro elements than number of mpi processes!\n");
  }

701

702
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
703
  {
704
    FUNCNAME("MeshDistributor::synchVector()");
705
706
707
708
709
710
711
   
    int nLevels = levelData.getNumberOfLevels();
    for (int level = nLevels - 1; level >= 0; level--) {
      MPI::Intracomm &mpiComm = levelData.getMpiComm(level);
      
      if (mpiComm == MPI::COMM_SELF)
	continue;
712

713
      StdMpi<vector<double> > stdMpi(mpiComm);
714

715
716
717
718
719
720
721
722
723
724
725
726
727
728
      for (DofComm::Iterator it(dofComm[level].getSendDofs()); 
	   !it.end(); it.nextRank()) {
	vector<double> dofs;
	
	for (int i = 0; i < vec.getSize(); i++) {
	  DOFVector<double> &dofVec = *(vec.getDOFVector(i));
	  for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
	    dofs.push_back(dofVec[it.getDofIndex()]);
	}
	
	int rank = it.getRank();
	// if (level > 0)
	//   rank = levelData.mapRank(rank, 0, level);
	stdMpi.send(rank, dofs);
729
      }
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
      
      for (DofComm::Iterator it(dofComm[level].getRecvDofs()); 
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	// if (level > 0)
	//   rank = levelData.mapRank(rank, 0, level);
	stdMpi.recv(rank);
      }
      
      stdMpi.startCommunication();
      
      for (DofComm::Iterator it(dofComm[level].getRecvDofs()); 
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	// if (level > 0)
	//   rank = levelData.mapRank(rank, 0, level);
	
	int counter = 0;
	vector<double> &recvData =  stdMpi.getRecvData(rank);
	
	for (int i = 0; i < vec.getSize(); i++) {
	  DOFVector<double> &dofVec = *(vec.getDOFVector(i));
	  
	  for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof()) {
	    TEST_EXIT_DBG(counter < recvData.size())
	      ("Recv data from rank %d has only %d entries!\n", rank, recvData.size());
	    dofVec[it.getDofIndex()] = recvData[counter++];
	  }
758
	}
759
760
761
762
      }
    }
  }

763

764
  void MeshDistributor::fix3dMeshRefinement()
765
  {
766
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
767
768
769
770

    if (mesh->getDim() != 3)
      return;

771
    Timer t;
772

773
774
775
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

776
777
778
779
780
781
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
782
      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
783
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
784
785
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
786
787
788
789
790
791
792
      }

      rankMacroEls.insert(elInfo->getElement()->getIndex());

      elInfo = stack.traverseNext(elInfo);
    }

793

794
795
    // === Reset fixed edges, and start new search for edges that must  ===
    // === be fixed.                                                    ===
796

797
    FixRefinementPatch::connectedEdges.clear();
798

799
    // Traverse all edges
800
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
801
      // For each edge get all macro elements that contain this edge.
802
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
803
804
805
806

      TEST_EXIT_DBG(edgeEls.size() > 0)
	("No edge %d/%d in elObjDB!\n", it->first, it->second);

807
808

      // We have now to check, if the current edge connects two macro elements,
809
      // which are otherwise not connected. 
810

811
812
813
      // All elements on this edge.
      std::set<int> allEls;
      // Maps from each element index to the local edge index of the common edge.
814
      map<int, int> edgeNoInEl;
815
816
817

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
818
	  allEls.insert(edgeEls[i].elIndex);
819
820
821
	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
822
           
823
      TEST_EXIT_DBG(!allEls.empty())("Should not happen!\n");
824

825
826
827
828

      // If there is only one element in rank on this edge, there is nothing 
      // to do.
      if (allEls.size() == 1)
829
830
	continue;
      
831
832
833
834
835
      // Create a set of element index sets. All element within one set are 
      // neighbours of each other, but can not reach any other element of the
      // other index sets by neighbourhood relations.
      vector<std::set<int> > disconnectedEls;
      helpToFix(allEls, edgeNoInEl, disconnectedEls);   
836
      
837
838
839
840
841
842
843
844
845
846
847
848
      // Fix the edges
      for (int dc0 = 0; dc0 < static_cast<int>(disconnectedEls.size()); dc0++) {
	for (int dc1 = 0; dc1 < static_cast<int>(disconnectedEls.size()); dc1++) {
	  if (dc0 == dc1)
	    continue;   
	  
	  int elIdx = *(disconnectedEls[dc1].begin());
	  pair<Element*, int> edge0 = 
	    make_pair(elObjDb.getElementPtr(elIdx), edgeNoInEl[elIdx]);
	  
	  for (std::set<int>::iterator elIt = disconnectedEls[dc0].begin();
	       elIt != disconnectedEls[dc0].end(); ++elIt) {
849
	    pair<Element*, int> edge1 = 
850
	      make_pair(elObjDb.getElementPtr(*elIt), edgeNoInEl[*elIt]);
851

852
#if (DEBUG != 0)
853
854
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);
855
	    
856
	    WorldVector<double> c0, c1;
857
858
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
859
860
861
862
863
864
865
866
867
868
869
870
	    
	    // 	      MSG("FOUND EDGE %d/%d <-> %d/%d\n", 
	    // 		  edge0.first->getIndex(), edge0.second,
	    // 		  edge1.first->getIndex(), edge1.second);
	    
	    // 	      MSG("FOUND EDGE: %d %d with coords %f %f %f   %f %f %f\n",
	    // 		  dofEdge0.first, dofEdge0.second, 
	    // 		  c0[0], c0[1], c0[2],
	    // 		  c1[0], c1[1], c1[2]);
	    
	    TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n");
#endif
871

872
873
874
875
876
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}	  
      }
    }
877

878

879
    MSG("Fix 3D mesh refinement needed %.5f seconds\n", t.elapsed());
880
  }
881

882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901

  void MeshDistributor::helpToFix(std::set<int> &elems, 
				  map<int, int> &edgeNoInEl,
				  vector<std::set<int> > &disconnectedEls)
  {
    std::set<int> firstElem;
    firstElem.insert(*(elems.begin()));

    std::set<int> otherElems(++(elems.begin()), elems.end());
    
    bool found = false;      
    do {
      found = false;
      for (std::set<int>::iterator elIt = firstElem.begin(); elIt != firstElem.end(); ++elIt) {
	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
	  int neighEl = macroElementNeighbours[*elIt][i];
	  if (neighEl != -1 && otherElems.count(neighEl)) {
	    firstElem.insert(neighEl);
	    otherElems.erase(neighEl);
	    found = true;
902
903
	  }
	}
904
905
906
	
	if (found)
	  break;
907
      }
908
    } while (found);
909

910
911
912
    disconnectedEls.push_back(firstElem);
    if (otherElems.size())
      helpToFix(otherElems, edgeNoInEl, disconnectedEls);
913
914
915
  }


916
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,