MeshDistributor.cc 68.9 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
48
49
50
51
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
52
53
#include "DOFMatrix.h"
#include "DOFVector.h"
54
#include "SystemVector.h"
55
#include "ElementDofIterator.h"
56
57
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
58
#include "VertexVector.h"
59
#include "MeshStructure.h"
60
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
61
#include "ProblemInstat.h"
62
#include "RefinementManager3d.h"
63
#include "Debug.h"
64
#include "Timer.h"
65

66
namespace AMDiS { namespace Parallel {
67

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

72
73
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

74
75
76
  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
77

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

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

111
    MPI::Intracomm &mpiComm = MPI::COMM_WORLD;
112
    mpiRank = mpiComm.Get_rank();
113

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
142
    if (partStr == "simple")
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
      partitioner = new SimplePartitioner("parallel->partitioner", &mpiComm);


    // === 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
163

164
    int tmp = 0;
165
    Parameters::get(name + "->box partitioning", tmp);
Thomas Witkowski's avatar
Thomas Witkowski committed
166
    partitioner->setBoxPartitioning(static_cast<bool>(tmp));
167
    initialPartitioner->setBoxPartitioning(static_cast<bool>(tmp));
Thomas Witkowski's avatar
Thomas Witkowski committed
168

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

172
173
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
174
175
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
176

Thomas Witkowski's avatar
Thomas Witkowski committed
177
178
  MeshDistributor::~MeshDistributor()
  {
179
    if (partitioner) {
Thomas Witkowski's avatar
Thomas Witkowski committed
180
      delete partitioner;
181
182
      partitioner = NULL;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
185
  }


186
  void MeshDistributor::initParallelization()
187
  {
188
    FUNCNAME("MeshDistributor::initParallelization()");
189

190
191
192
    if (initialized)
      return;

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

197
    TEST_EXIT(MPI::COMM_WORLD.Get_size() > 1)
198
      ("Parallelization does not work with only one process!\n");
199
    TEST_EXIT(feSpaces.size() > 0)
200
      ("No FE space has been defined for the mesh distributor!\n");
201
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
202

203
204
205
206
207
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
208
      doNext = false;
209
210
211
212
213
214
      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;
215
216
	  doNext = true;
	}
217
      }
218
219
    } while (doNext);

220
    elObjDb.setFeSpace(feSpaces[0]);
221

Thomas Witkowski's avatar
Thomas Witkowski committed
222

223
224
225
    // 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).
226
    if (deserialized) {
227
228
      updateMacroElementInfo();

229
      removePeriodicBoundaryConditions();
230

231
      elObjDb.createMacroElementInfo(allMacroElements);
232

Thomas Witkowski's avatar
Thomas Witkowski committed
233
      updateLocalGlobalNumbering();
234

Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
237
#if (DEBUG != 0)
238
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
239
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], 
240
				    *(dofMaps[0]), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
241
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
242
#endif    
243
244

      initialized = true;
245
      return;
246
    }
247

248

Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
    // Create a very first pariitioning of the mesh.
    createInitialPartitioning();
251

252

253
#if (DEBUG != 0)    
254
255
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
256
257
#endif

258
    if (mpiRank == 0) {
259
#if (DEBUG != 0)    
260
      int writePartMesh = 1;
261
262
263
#else
      int writePartMesh = 0;
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
264
      Parameters::get("parallel->debug->write mesh partitioning", writePartMesh);
265

266
      if (writePartMesh > 0) {
267
	debug::writeElementIndexMesh(mesh , debugOutputDir + "elementIndex");
268
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part");
269
      }
270
    }
271

272
    // Create interior boundary information.
273
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
274

275
276
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

277
278
279
280
281
    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))) {
282
283
284

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

285
286
287
288
289
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
290
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
291
292
293
294
295
296
297
298
299
	}
      }
    }

    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))) {
300
301
302

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

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

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

310
311
312
	}
      }
    }
313

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

    removeMacroElements();
317

318
319
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
320

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

328

329
    updateLocalGlobalNumbering();
330

331

332
333
334
    // === In 3D we have to fix the mesh to allow local refinements. ===

    fix3dMeshRefinement();
335

336

337
338
    // === If in debug mode, make some tests. ===

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

342
    ParallelDebug::testAllElements(*this);
343
    debug::testSortedDofs(mesh, elMap);
344
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
345
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
346

347
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
348
349

    MSG("Debug mode tests finished!\n");
350
#endif
351

352
353
354
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

355
356
357
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
358

359
    // === Global refinements. ===
360
    
Thomas Witkowski's avatar
Thomas Witkowski committed
361
    int globalRefinement = 0;
362
    Parameters::get(mesh->getName() + "->global refinements", globalRefinement);
363
    
Thomas Witkowski's avatar
Thomas Witkowski committed
364
    if (globalRefinement > 0) {
365
      refineManager->globalRefine(mesh, globalRefinement);
366
     
367
      updateLocalGlobalNumbering();
368

369
370
371
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
372
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
373

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
374
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
375
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
376
377
      elObjDb.clear();

378
    initialized = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
379
    MSG("Init parallelization needed %.5f seconds\n", MPI::Wtime() - first);
380
381
  }

382

Thomas Witkowski's avatar
Thomas Witkowski committed
383
384
385
386
387
388
389
390
391
392
393
394
395
396
  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
397
398
    bool useInitialPartitioning = 
      initialPartitioner->createInitialPartitioning();
Thomas Witkowski's avatar
Thomas Witkowski committed
399
400
401
402

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

403
    if (!useInitialPartitioning) {   
Thomas Witkowski's avatar
Thomas Witkowski committed
404
      // and now partition the mesh    
405
406
      bool partitioningSucceed = 
	initialPartitioner->partition(elemWeights, INITIAL);
Thomas Witkowski's avatar
Thomas Witkowski committed
407
408
      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
    }
409

410
411
412
413
414
    initialPartitioner->createPartitionMap(partitionMap);

    if (initialPartitioner != partitioner) {
      *partitioner = *initialPartitioner;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
  }


  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;

452
      /*int nProc = */ArhReader::readMetaData(filename, arhElInRank, arhElCodeSize);
Thomas Witkowski's avatar
Thomas Witkowski committed
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
      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);    
    }
  }


469
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
470
  {
471
    FUNCNAME("MeshDistributor::addProblemStat()");
472

473
474
475
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

476

477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
    // === 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.                                 ===
    
    if (mesh != NULL) {
      TEST_EXIT(mesh == probStat->getMesh())
	("Does not yet support for different meshes!\n");
    } else {
      mesh = probStat->getMesh();
    }
494
495
496
497
498
499
500
501
502
503
    
    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());
504
    }
505
506
    
    partitioner->setMesh(mesh);
507
    initialPartitioner->setMesh(mesh);
508
    
509

510
511
512
    // === Check whether the stationary problem should be serialized. ===


513
514
    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
515
    Parameters::get(probStat->getName() + "->output->write serialization",  
516
		    writeSerialization);
517
    if (writeSerialization && !writeSerializationFile) {
518
      string filename = "";
519
      Parameters::get(name + "->output->serialization filename", filename);
520
521
522
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
523
524

      int tsModulo = -1;
525
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
526
		      tsModulo);
527
      
528
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
529
530
      writeSerializationFile = true;
    }    
531

532
533
534

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

535
    int readSerialization = 0;
536
    Parameters::get(probStat->getName() + "->input->read serialization", 
537
		    readSerialization);
538
    if (readSerialization) {
539
      string filename = "";
540
      Parameters::get(probStat->getName() + "->input->serialization filename", 
541
		      filename);
542
      filename += ".p" + lexical_cast<string>(mpiRank);
543
      MSG("Start deserialization with %s\n", filename.c_str());
544
      ifstream in(filename.c_str());
545
546
547
548

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

549
550
551
552
553
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

554
      probStat->deserialize(in);
555
      in.close();
556
557
      MSG("Deserialization from file: %s\n", filename.c_str());

558
      filename = "";
559
      Parameters::get(name + "->input->serialization filename", filename);
560
561
562
563
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
564
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
565
566
567
568
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
569
570
571
572
573

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
574
575
576
577
578
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
579
580
    }

581
    problemStat.push_back(probStat);
582

583
584
585

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

587
    if (initialized)
588
      removePeriodicBoundaryConditions(probStat);
589
590
591
  }


592
593
594
595
596
597
598
599
600
601
602
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

    if (globalMeshDistributor == NULL)
      globalMeshDistributor = new MeshDistributor();

    globalMeshDistributor->addProblemStat(probStat);
  }


603
  void MeshDistributor::exitParallelization()
604
  {}
605

606
607
608
609
610
611
612
613
614

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

615
    dofMaps.push_back(&dofMap);   
616
617
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
618
619
620
621
622
623

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

624
625
    if (it != dofMaps.end())
      dofMaps.erase(it);
Thomas Witkowski's avatar
Thomas Witkowski committed
626
627
628
  }


629
  void MeshDistributor::testForMacroMesh()
630
  {
631
    FUNCNAME("MeshDistributor::testForMacroMesh()");
632
633
634
635
636
637
638

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
643
644
645
646
647
648
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

649
    TEST_EXIT(nMacroElements >= MPI::COMM_WORLD.Get_size())
650
651
652
      ("The mesh has less macro elements than number of mpi processes!\n");
  }

653

654
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
655
  {
656
    FUNCNAME("MeshDistributor::synchVector()");
657
658
659
660
661
662
663
   
    int nLevels = levelData.getNumberOfLevels();
    for (int level = nLevels - 1; level >= 0; level--) {
      MPI::Intracomm &mpiComm = levelData.getMpiComm(level);
      
      if (mpiComm == MPI::COMM_SELF)
	continue;
664

665
      StdMpi<vector<double> > stdMpi(mpiComm);
666

667
668
669
670
671
672
673
674
675
676
677
678
679
680
      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);
681
      }
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
      
      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++];
	  }
710
	}
711
712
713
714
      }
    }
  }

715

716
  void MeshDistributor::fix3dMeshRefinement()
717
  {
718
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
719
720
721
722

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

723
    Timer t;
724

725
726
727
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

728
729
730
731
732
733
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
734
      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
735
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
736
737
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
738
739
740
741
742
743
744
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

745

746
747
    // === Reset fixed edges, and start new search for edges that must  ===
    // === be fixed.                                                    ===
748

749
    FixRefinementPatch::connectedEdges.clear();
750

751
    // Traverse all edges
752
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
753
      // For each edge get all macro elements that contain this edge.
754
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
755
756
757
758

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

759
760

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

763
764
765
      // All elements on this edge.
      std::set<int> allEls;
      // Maps from each element index to the local edge index of the common edge.
766
      map<int, int> edgeNoInEl;
767
768
769

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
770
	  allEls.insert(edgeEls[i].elIndex);
771
772
773
	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
774
           
775
      TEST_EXIT_DBG(!allEls.empty())("Should not happen!\n");
776

777
778
779
780

      // If there is only one element in rank on this edge, there is nothing 
      // to do.
      if (allEls.size() == 1)
781
782
	continue;
      
783
784
785
786
787
      // 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);   
788
      
789
790
791
792
793
794
795
796
797
798
799
800
      // 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) {
801
	    pair<Element*, int> edge1 = 
802
	      make_pair(elObjDb.getElementPtr(*elIt), edgeNoInEl[*elIt]);
803

804
#if (DEBUG != 0)
805
806
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);
807
	    
808
	    WorldVector<double> c0, c1;
809
810
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
811
812
813
814
815
816
817
818
819
820
821
822
	    
	    // 	      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
823

824
825
826
827
828
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}	  
      }
    }
829

830

831
    MSG("Fix 3D mesh refinement needed %.5f seconds\n", t.elapsed());
832
  }
833

834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853

  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;
854
855
	  }
	}
856
857
858
	
	if (found)
	  break;
859
      }
860
    } while (found);
861

862
863
864
    disconnectedEls.push_back(firstElem);
    if (otherElems.size())
      helpToFix(otherElems, edgeNoInEl, disconnectedEls);
865
866
867
  }


868
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
869
					   int level,
870
					   DofContainer& dofs)
871
872
  {
    DofContainerSet dofSet;
873
    for (DofComm::Iterator it(dofComm[level].getSendDofs(), feSpace); 
874
	 !it.end(); it.nextRank())
875
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
876
    for (DofComm::Iterator it(dofComm[level].getRecvDofs(), feSpace); 
877
	 !it.end(); it.nextRank())
878
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
879
880
881
882
883
884

    dofs.clear();
    dofs.insert(dofs.begin(), dofSet.begin(), dofSet.end());
  }


885
886
887
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
    // Remove periodic boundaries in boundary manager on matrices and vectors.
888
889
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
890
891
892

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
893
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
894
895
896
897
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
898
899
900

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
901
902
903
  }


904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
  void MeshDistributor::removePeriodicBoundaryConditions(ProblemStatSeq *probStat)
  {
    int nComponents = probStat->getNumComponents();

    for (int j = 0; j < nComponents; j++) {
      for (int k = 0; k < nComponents; k++) {
	DOFMatrix* mat = probStat->getSystemMatrix(j, k);
	if (mat && mat->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
      }
      
      if (probStat->getSolution()->getDOFVector(j)->getBoundaryManager())
	removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
      
      if (probStat->getRhs()->getDOFVector(j)->getBoundaryManager())
	removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat->getRhs()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
    }
  }


924
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
925
926
927
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {