MeshDistributor.cc 68.1 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
  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
Thomas Witkowski's avatar
Thomas Witkowski committed
360
    bool useInitialPartitioning = partitioner->createInitialPartitioning();
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
363
364

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

Thomas Witkowski's avatar
Thomas Witkowski committed
365
    if (!useInitialPartitioning) {
Thomas Witkowski's avatar
Thomas Witkowski committed
366
367
368
369
      // and now partition the mesh    
      bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
    }
370
371

    partitioner->createPartitionMap(partitionMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
372
373
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
  }


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


426
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
427
  {
428
    FUNCNAME("MeshDistributor::addProblemStat()");
429

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

433

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

467
468
469
    // === Check whether the stationary problem should be serialized. ===


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

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

489
490
491

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

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

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

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

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

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

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

538
    problemStat.push_back(probStat);
539

540
541
542

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

544
    if (initialized)
545
      removePeriodicBoundaryConditions(probStat);
546
547
548
  }


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

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

    globalMeshDistributor->addProblemStat(probStat);
  }


560
  void MeshDistributor::exitParallelization()
561
  {}
562

563
564
565
566
567
568
569
570
571
572
573
574

  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
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589

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


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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

614

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

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

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

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

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

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
636
637
638
639
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
640
    }
641
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
642
643
644
645
646
647
    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);
    }
648

649
    stdMpi.startCommunication();
650

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

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

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

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

671

672
  void MeshDistributor::fix3dMeshRefinement()
673
  {
674
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
675
676
677
678

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

679
680
681
682
    MPI::COMM_WORLD.Barrier();
    double first = MPI::Wtime();


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

686
687
688
689
690
691
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

703

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

707
    FixRefinementPatch::connectedEdges.clear();
708

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

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

717
718

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

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

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

735
736
737
738

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

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

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

788

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

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

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

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


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

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

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


846
847
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
848
849
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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

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

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


867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
  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()));
    }
  }


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

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


903
904
905
906
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

907
908
909
    if (mpiSize != 16)
      return;

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