MeshDistributor.cc 69.3 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

65
namespace AMDiS { namespace Parallel {
66

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

71
72
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
141
    if (partStr == "simple")
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
      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
162

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
175

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


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

189
190
191
    if (initialized)
      return;

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
221

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

228
      removePeriodicBoundaryConditions();
229

230
      elObjDb.createMacroElementInfo(allMacroElements);
231

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

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

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

      initialized = true;
244
      return;
245
    }
246

247

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

251

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

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

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

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

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

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

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

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

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

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

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

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

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

309
310
311
	}
      }
    }
312

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

    removeMacroElements();
316

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

Thomas Witkowski's avatar
Thomas Witkowski committed
319

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

327

328
    updateLocalGlobalNumbering();
329

330

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

    fix3dMeshRefinement();
334

335

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

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

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

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

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

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

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

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

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

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

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

381

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

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

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

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

    if (initialPartitioner != partitioner) {
      *partitioner = *initialPartitioner;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
414
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
  }


  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;

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


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

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

475

476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
    // === 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();
    }
493
494
495
496
497
498
499
500
501
502
    
    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());
503
    }
504
505
    
    partitioner->setMesh(mesh);
506
    initialPartitioner->setMesh(mesh);
507
    
508

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


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

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

531
532
533

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

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

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

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

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

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

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

580
    problemStat.push_back(probStat);
581

582
583
584

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

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


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

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

    globalMeshDistributor->addProblemStat(probStat);
  }


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

605
606
607
608
609
610
611
612
613

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

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

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

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

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


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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

652

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

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

666
667
668
669
670
671
672
673
674
675
676
677
678
679
      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);
680
      }
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
      
      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++];
	  }
709
	}
710
711
712
713
      }
    }
  }

714

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

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

722
723
724
725
    MPI::COMM_WORLD.Barrier();
    double first = MPI::Wtime();


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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

746

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

750
    FixRefinementPatch::connectedEdges.clear();
751

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

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

760
761

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

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

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

778
779
780
781

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

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

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

831

832
833
834
    MPI::COMM_WORLD.Barrier();
    MSG("Fix 3D mesh refinement needed %.5f seconds\n", MPI::Wtime() - first);
  }
835

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

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

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


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

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


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

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

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


906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
  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()));
    }
  }


926
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
927
928
929
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
930
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
931
	boundaryMap.erase(it++);
Thomas Witkowski's avatar