MeshDistributor.cc 67.7 KB
Newer Older
1
//
2
3
4
5
6
7
8
9
10
11
12
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include <algorithm>
14
15
#include <iostream>
#include <fstream>
16
17
#include <limits>
#include <stdint.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
18
#include <boost/lexical_cast.hpp>
19
20
#include <boost/filesystem.hpp>

21
#include "parallel/MeshDistributor.h"
22
#include "parallel/MeshManipulation.h"
23
#include "parallel/ParallelDebug.h"
24
#include "parallel/StdMpi.h"
25
#include "parallel/MeshPartitioner.h"
26
#include "parallel/ParMetisPartitioner.h"
27
#include "parallel/ZoltanPartitioner.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28
#include "parallel/SimplePartitioner.h"
29
#include "parallel/CheckerPartitioner.h"
30
#include "parallel/MpiHelper.h"
31
#include "parallel/DofComm.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
32
#include "parallel/ParallelProblemStatBase.h"
33
34
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
35
#include "io/MacroWriter.h"
36
#include "io/VtkWriter.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
37
#include "io/ArhReader.h"
38
39
40
41
42
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
43
44
#include "DOFMatrix.h"
#include "DOFVector.h"
45
#include "SystemVector.h"
46
#include "ElementDofIterator.h"
47
48
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
49
#include "VertexVector.h"
50
#include "MeshStructure.h"
51
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
52
#include "ProblemInstat.h"
53
#include "RefinementManager3d.h"
54
#include "Debug.h"
55

56
57
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
58
  using boost::lexical_cast;
59
  using namespace boost::filesystem;
60
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
61

62
63
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

64
65
66
  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
67

68
69
70
71
72
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

73
  MeshDistributor::MeshDistributor()
74
    : problemStat(0),
75
      initialized(false),
76
      name("parallel"),
77
78
79
      mesh(NULL),
      refineManager(NULL),
      partitioner(NULL),
80
      deserialized(false),
81
      writeSerializationFile(false),
82
      repartitioningAllowed(false),
83
      repartitionIthChange(20),
84
      repartitioningWaitAfterFail(20),
85
86
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
87
      repartitioningFailed(0),
88

89
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
90
      lastMeshChangeIndex(0),
91
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
92
      boundaryDofInfo(1),
93
      meshAdaptivity(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
94
95
      hasPeriodicBoundary(false),
      printTimings(false)
96
  {
97
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
98

99
    mpiComm = MPI::COMM_WORLD;
100
101
    mpiRank = mpiComm.Get_rank();
    mpiSize = mpiComm.Get_size();
102

103
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
104
105
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
106
    Parameters::get(name + "->repartition wait after fail", repartitioningWaitAfterFail);
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
107
    Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
Thomas Witkowski's avatar
Thomas Witkowski committed
108
    
109
    string partStr = "parmetis";
110
    Parameters::get(name + "->partitioner", partStr);
111
112
113
114

    if (partStr == "parmetis") 
      partitioner = new ParMetisPartitioner(&mpiComm);

Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
    if (partStr == "zoltan") {
#ifdef HAVE_ZOLTAN
117
      partitioner = new ZoltanPartitioner(&mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
118
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
119
      ERROR_EXIT("AMDiS was compiled without Zoltan support. Therefore you cannot make use of it!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
#endif
    }
122

123
124
125
    if (partStr == "checker")
      partitioner = new CheckerPartitioner(&mpiComm);

Thomas Witkowski's avatar
Thomas Witkowski committed
126
127
128
    if (partStr == "simple")
      partitioner = new SimplePartitioner(&mpiComm);

129
    int tmp = 0;
130
    Parameters::get(name + "->box partitioning", tmp);
Thomas Witkowski's avatar
Thomas Witkowski committed
131
132
    partitioner->setBoxPartitioning(static_cast<bool>(tmp));

Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
    Parameters::get(name + "->print timings", printTimings);

135
    TEST_EXIT(partitioner)("Could not create partitioner \"%s\"!\n", partStr.c_str());
136
137
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
138

Thomas Witkowski's avatar
Thomas Witkowski committed
139
140
141
142
143
144
145
  MeshDistributor::~MeshDistributor()
  {
    if (partitioner)
      delete partitioner;
  }


146
  void MeshDistributor::initParallelization()
147
  {
148
    FUNCNAME("MeshDistributor::initParallelization()");
149

150
151
152
    if (initialized)
      return;

Thomas Witkowski's avatar
Thomas Witkowski committed
153
154
155
156
    double first = MPI::Wtime();
    MSG("Initialization phase 1 needed %.5f seconds\n", 
	first - ParallelProblemStatBase::initTimeStamp);

157
158
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");
159
    TEST_EXIT(feSpaces.size() > 0)
160
      ("No FE space has been defined for the mesh distributor!\n");
161
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
162

163
164
165
166
167
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
168
      doNext = false;
169
170
171
172
173
174
      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;
175
176
	  doNext = true;
	}
177
      }
178
179
    } while (doNext);

180
    elObjDb.setFeSpace(feSpaces[0]);
181

Thomas Witkowski's avatar
Thomas Witkowski committed
182
183
184
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

185
186
187
    // 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).
188
    if (deserialized) {
Thomas Witkowski's avatar
Thomas Witkowski committed
189
190
      createMeshLevelStructure();
      
191
192
      updateMacroElementInfo();

193
      removePeriodicBoundaryConditions();
194

195
      elObjDb.createMacroElementInfo(allMacroElements);
196

Thomas Witkowski's avatar
Thomas Witkowski committed
197
      updateLocalGlobalNumbering();
198

Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
201
#if (DEBUG != 0)
202
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
203
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], 
204
				    *(dofMaps[0]), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
205
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
206
#endif    
207
208

      initialized = true;
209
      return;
210
    }
211

212

Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
    // Create a very first pariitioning of the mesh.
    createInitialPartitioning();
215

216

217
#if (DEBUG != 0)    
218
219
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
220
221
#endif

222
    if (mpiRank == 0) {
223
#if (DEBUG != 0)    
224
      int writePartMesh = 1;
225
226
227
#else
      int writePartMesh = 0;
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
228
      Parameters::get("parallel->debug->write mesh partitioning", writePartMesh);
229

230
      if (writePartMesh > 0) {
231
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex");
232
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
233
      }
234
    }
235

236
    // Create interior boundary information.
237
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
238

239
240
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

241
242
243
244
245
    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))) {
246
247
248

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

249
250
251
252
253
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
254
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
255
256
257
258
259
260
261
262
263
	}
      }
    }

    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))) {
264
265
266

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

267
268
269
270
271
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

274
275
276
	}
      }
    }
277

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

    removeMacroElements();
281

282
283
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
284

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

292
    updateLocalGlobalNumbering();
293

294

295
296
297
    // === In 3D we have to fix the mesh to allow local refinements. ===

    fix3dMeshRefinement();
298

299

300
301
    // === If in debug mode, make some tests. ===

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

305
    ParallelDebug::testAllElements(*this);
306
    debug::testSortedDofs(mesh, elMap);
307
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
308
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
309

310
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
311
312

    MSG("Debug mode tests finished!\n");
313
#endif
314

315
316
317
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

318
319
320
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
321

322
    // === Global refinements. ===
323
    
Thomas Witkowski's avatar
Thomas Witkowski committed
324
    int globalRefinement = 0;
325
    Parameters::get(mesh->getName() + "->global refinements", globalRefinement);
326
    
Thomas Witkowski's avatar
Thomas Witkowski committed
327
    if (globalRefinement > 0) {
328
      refineManager->globalRefine(mesh, globalRefinement);
329
     
330
      updateLocalGlobalNumbering();
331

332
333
334
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
335
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
336

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
337
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
338
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
339
340
      elObjDb.clear();

341
    initialized = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
342
    MSG("Init parallelization needed %.5f seconds\n", MPI::Wtime() - first);
343
344
  }

345

Thomas Witkowski's avatar
Thomas Witkowski committed
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
  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
    partitioner->createInitialPartitioning();

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

    string filename = "";
    Parameters::get("parallel->partitioner->read meta arh", filename);
    if (filename == "" || ArhReader::readMetaData(filename) != mpiSize) {
      // and now partition the mesh    
      bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
    }
372
373

    partitioner->createPartitionMap(partitionMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
  }


  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;

      int nProc = ArhReader::readMetaData(filename, arhElInRank, arhElCodeSize);
      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);    
    }
  }


428
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
429
  {
430
    FUNCNAME("MeshDistributor::addProblemStat()");
431

432
433
434
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

435

436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
    // === 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();
    }
454
455
456
457
458
459
460
461
462
463
    
    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());
464
    }
465
466
467
    
    partitioner->setMesh(mesh);
    
468

469
470
471
    // === Check whether the stationary problem should be serialized. ===


472
473
    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
474
    Parameters::get(probStat->getName() + "->output->write serialization",  
475
		    writeSerialization);
476
    if (writeSerialization && !writeSerializationFile) {
477
      string filename = "";
478
      Parameters::get(name + "->output->serialization filename", filename);
479
480
481
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
482
483

      int tsModulo = -1;
484
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
485
		      tsModulo);
486
      
487
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
488
489
      writeSerializationFile = true;
    }    
490

491
492
493

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

494
    int readSerialization = 0;
495
    Parameters::get(probStat->getName() + "->input->read serialization", 
496
		    readSerialization);
497
    if (readSerialization) {
498
      string filename = "";
499
      Parameters::get(probStat->getName() + "->input->serialization filename", 
500
		      filename);
501
      filename += ".p" + lexical_cast<string>(mpiRank);
502
      MSG("Start deserialization with %s\n", filename.c_str());
503
      ifstream in(filename.c_str());
504
505
506
507

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

508
509
510
511
512
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

513
      probStat->deserialize(in);
514
      in.close();
515
516
      MSG("Deserialization from file: %s\n", filename.c_str());

517
      filename = "";
518
      Parameters::get(name + "->input->serialization filename", filename);
519
520
521
522
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
523
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
524
525
526
527
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
528
529
530
531
532

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
533
534
535
536
537
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
538
539
    }

540
    problemStat.push_back(probStat);
541

542
543
544

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

546
    if (initialized)
547
      removePeriodicBoundaryConditions(probStat);
548
549
550
  }


551
552
553
554
555
556
557
558
559
560
561
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


562
  void MeshDistributor::exitParallelization()
563
  {}
564

565
566
567
568
569
570
571
572
573
574
575
576

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

    dofMaps.push_back(&dofMap);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591

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

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

    TEST_EXIT(it != dofMaps.end())
      ("Cannot find Parallel DOF mapping object which should be removed!\n");

    dofMaps.erase(it);
  }


592
  void MeshDistributor::testForMacroMesh()
593
  {
594
    FUNCNAME("MeshDistributor::testForMacroMesh()");
595
596
597
598
599
600
601

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
606
607
608
609
610
611
612
613
614
615
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

    TEST_EXIT(nMacroElements >= mpiSize)
      ("The mesh has less macro elements than number of mpi processes!\n");
  }

616

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
617
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
618
  {
619
620
    FUNCNAME("MeshDistributor::synchVector()");

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
621
    TEST_EXIT_DBG(level >= 0 && level <= 1)("Wrong level number!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
622

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
623
624
625
626
627
628
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

    StdMpi<vector<double> > stdMpi(levelComm);

    for (DofComm::Iterator it(dc.getSendDofs()); 
629
	 !it.end(); it.nextRank()) {
630
      vector<double> dofs;
631

632
633
      for (int i = 0; i < vec.getSize(); i++) {
	DOFVector<double> &dofVec = *(vec.getDOFVector(i));
Thomas Witkowski's avatar
Thomas Witkowski committed
634
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
635
636
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
637

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
638
639
640
641
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
642
    }
643
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
644
645
646
647
648
649
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.recv(rank);
    }
650

651
    stdMpi.startCommunication();
652

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
653
654
655
656
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
657
658

      int counter = 0;
659
      vector<double> &recvData =  stdMpi.getRecvData(rank);
660
661
662
663

      for (int i = 0; i < vec.getSize(); i++) {
	DOFVector<double> &dofVec = *(vec.getDOFVector(i));

664
665
666
667
668
	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++];
	}
669
670
671
672
      }
    }
  }

673

674
  void MeshDistributor::fix3dMeshRefinement()
675
  {
676
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
677
678
679
680

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

681
682
683
684
    MPI::COMM_WORLD.Barrier();
    double first = MPI::Wtime();


685
686
687
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

688
689
690
691
692
693
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
694
      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
695
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
696
697
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
698
699
700
701
702
703
704
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

705

706
707
    // === Reset fixed edges, and start new search for edges that must  ===
    // === be fixed.                                                    ===
708

709
    FixRefinementPatch::connectedEdges.clear();
710

711
    // Traverse all edges
712
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
713
      // For each edge get all macro elements that contain this edge.
714
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
715
716
717
718

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

719
720

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

723
724
725
      // All elements on this edge.
      std::set<int> allEls;
      // Maps from each element index to the local edge index of the common edge.
726
      map<int, int> edgeNoInEl;
727
728
729

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
730
	  allEls.insert(edgeEls[i].elIndex);
731
732
733
	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
734
           
735
      TEST_EXIT_DBG(!allEls.empty())("Should not happen!\n");
736

737
738
739
740

      // If there is only one element in rank on this edge, there is nothing 
      // to do.
      if (allEls.size() == 1)
741
742
	continue;
      
743
744
745
746
747
      // 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);   
748
      
749
750
751
752
753
754
755
756
757
758
759
760
      // 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) {
761
	    pair<Element*, int> edge1 = 
762
	      make_pair(elObjDb.getElementPtr(*elIt), edgeNoInEl[*elIt]);
763

764
#if (DEBUG != 0)
765
766
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);
767
	    
768
	    WorldVector<double> c0, c1;
769
770
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
771
772
773
774
775
776
777
778
779
780
781
782
	    
	    // 	      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
783

784
785
786
787
788
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}	  
      }
    }
789

790

791
792
793
    MPI::COMM_WORLD.Barrier();
    MSG("Fix 3D mesh refinement needed %.5f seconds\n", MPI::Wtime() - first);
  }
794

795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814

  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;
815
816
	  }
	}
817
818
819
	
	if (found)
	  break;
820
      }
821
    } while (found);
822

823
824
825
    disconnectedEls.push_back(firstElem);
    if (otherElems.size())
      helpToFix(otherElems, edgeNoInEl, disconnectedEls);
826
827
828
  }


829
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
830
					   int level,
831
					   DofContainer& dofs)
832
833
834
835
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
836
837
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
838
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
839
840
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
841
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
842
843
844
845
846
847

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


848
849
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
850
851
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

852
    // Remove periodic boundaries in boundary manager on matrices and vectors.
853
854
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
855
856
857

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
858
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
859
860
861
862
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
863
864
865

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
866
867
868
  }


869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
  void MeshDistributor::removePeriodicBoundaryConditions(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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


891
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
892
  {
893
894
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
895
896
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
897
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
898
	boundaryMap.erase(it++);
899
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
900
901
902
903
904
	++it;      
    }    
  }


905
906
907
908
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

909
910
911
    if (mpiSize != 16)
      return;

912
    std::set<int> neighbours;
913
914
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
915
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
916
917
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
918
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
919
920
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
921
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
922
923
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
924
925
926
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
927
928
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
929
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
930
931
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
932
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
933
934
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
935
936
937
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
938
939
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
940
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);