MeshDistributor.cc 69 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
      repartitionIthChange(20),
94
      repartitioningWaitAfterFail(20),
95
96
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
97
      repartitioningFailed(0),
98
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
99
      lastMeshChangeIndex(0),
100
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
101
      boundaryDofInfo(1),
102
      meshAdaptivity(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
103
      hasPeriodicBoundary(false),
104
105
      printTimings(false),
      printMemoryUsage(false)
106
  {
107
    FUNCNAME("MeshDistributor::MeshDistributor()");
Thomas Witkowski's avatar
Thomas Witkowski committed
108

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
118
    
119
120
    // === Create partitioner object. ===

121
    string partStr = "parmetis";
122
    Parameters::get(name + "->partitioner", partStr);
123
124

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

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

135
    if (partStr == "checker")
136
      partitioner = new CheckerPartitioner("parallel->partitioner", &mpiComm);
137

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

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

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

168
169
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
170
171
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
172

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


182
  void MeshDistributor::initParallelization()
183
  {
184
    FUNCNAME("MeshDistributor::initParallelization()");
185

186
187
188
    if (initialized)
      return;

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

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

199
200
201
202
203
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

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

216
    elObjDb.setFeSpace(feSpaces[0]);
217

Thomas Witkowski's avatar
Thomas Witkowski committed
218

219
220
221
    // 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).
222
    if (deserialized) {
223
224
      updateMacroElementInfo();

225
      removePeriodicBoundaryConditions();
226

227
      elObjDb.createMacroElementInfo(allMacroElements);
228

Thomas Witkowski's avatar
Thomas Witkowski committed
229
      updateLocalGlobalNumbering();
230

Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
      elObjDb.setData(partitionMap, levelData);

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

      initialized = true;
241
      return;
242
    }
243

244

Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
    // Create a very first pariitioning of the mesh.
    createInitialPartitioning();
247

248

249
#if (DEBUG != 0)    
250
251
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
252
253
#endif

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

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

268
    // Create interior boundary information.
269
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
270

271
272
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

273
274
275
276
277
    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))) {
278
279
280

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

281
282
283
284
285
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

    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))) {
296
297
298

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

299
300
301
302
303
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

306
307
308
	}
      }
    }
309

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

    removeMacroElements();
313

314
315
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
316

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

324

325
    updateLocalGlobalNumbering();
326

327

328
329
330
    // === In 3D we have to fix the mesh to allow local refinements. ===

    fix3dMeshRefinement();
331

332

333
334
    // === If in debug mode, make some tests. ===

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

338
    ParallelDebug::testAllElements(*this);
339
    debug::testSortedDofs(mesh, elMap);
340
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
341
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
342

343
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
344
345

    MSG("Debug mode tests finished!\n");
346
#endif
347

348
349
350
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

351
352
353
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
354

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

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

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

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

378

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

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

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

406
407
408
409
410
    initialPartitioner->createPartitionMap(partitionMap);

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


  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;

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


465
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
466
  {
467
    FUNCNAME("MeshDistributor::addProblemStat()");
468

469
470
471
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

472

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

506
507
508
    // === Check whether the stationary problem should be serialized. ===


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

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

528
529
530

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

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

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

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

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

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

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

577
    problemStat.push_back(probStat);
578

579
580
581

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

583
    if (initialized)
584
      removePeriodicBoundaryConditions(probStat);
585
586
587
  }


588
589
590
591
592
593
594
595
596
597
598
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


599
  void MeshDistributor::exitParallelization()
600
  {}
601

602
603
604
605
606
607
608
609
610

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

611
    dofMaps.push_back(&dofMap);   
612
613
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
614
615
616
617
618
619

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

620
621
    if (it != dofMaps.end())
      dofMaps.erase(it);
Thomas Witkowski's avatar
Thomas Witkowski committed
622
623
624
  }


625
  void MeshDistributor::testForMacroMesh()
626
  {
627
    FUNCNAME("MeshDistributor::testForMacroMesh()");
628
629
630
631
632
633
634

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

649

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

661
      StdMpi<vector<double> > stdMpi(mpiComm);
662

663
664
665
666
667
668
669
670
671
672
673
674
675
676
      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);
677
      }
678
679
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
      
      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++];
	  }
706
	}
707
708
709
710
      }
    }
  }

711

712
  void MeshDistributor::fix3dMeshRefinement()
713
  {
714
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
715
716
717
718

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

719
720
721
722
    MPI::COMM_WORLD.Barrier();
    double first = MPI::Wtime();


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

726
727
728
729
730
731
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

743

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

747
    FixRefinementPatch::connectedEdges.clear();
748

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

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

757
758

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

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

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

775
776
777
778

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

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

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

828

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

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

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

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


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

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


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

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

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


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


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


935
936
  void MeshDistributor::createMeshLevelStructure()
  {