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::MeshDistributor()");
Thomas Witkowski's avatar
Thomas Witkowski committed
98

99
    MPI::Intracomm &mpiComm = MPI::COMM_WORLD;
100
    mpiRank = mpiComm.Get_rank();
101

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

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

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

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

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

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

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

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

    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
138
139
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
140

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


148
  void MeshDistributor::initParallelization()
149
  {
150
    FUNCNAME("MeshDistributor::initParallelization()");
151

152
153
154
    if (initialized)
      return;

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

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

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

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

182
    elObjDb.setFeSpace(feSpaces[0]);
183

Thomas Witkowski's avatar
Thomas Witkowski committed
184

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) {
189
190
      updateMacroElementInfo();

191
      removePeriodicBoundaryConditions();
192

193
      elObjDb.createMacroElementInfo(allMacroElements);
194

Thomas Witkowski's avatar
Thomas Witkowski committed
195
      updateLocalGlobalNumbering();
196

Thomas Witkowski's avatar
Thomas Witkowski committed
197
198
      elObjDb.setData(partitionMap, levelData);

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

      initialized = true;
207
      return;
208
    }
209

210

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

214

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

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

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

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

237
238
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

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

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

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

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

    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))) {
262
263
264

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

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

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

272
273
274
	}
      }
    }
275

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

    removeMacroElements();
279

280
281
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
282

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

290

291
    updateLocalGlobalNumbering();
292

293

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

    fix3dMeshRefinement();
297

298

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

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

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

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

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

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

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

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

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

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

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

344

Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
347
348
349
350
351
352
353
354
355
356
357
358
  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
359
    bool useInitialPartitioning = partitioner->createInitialPartitioning();
Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
362
363

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

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

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


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


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

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

432

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

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


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

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

488
489
490

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

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

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

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

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

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

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

537
    problemStat.push_back(probStat);
538

539
540
541

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

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


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

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

    globalMeshDistributor->addProblemStat(probStat);
  }


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

562
563
564
565
566
567
568
569
570

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

571
    dofMaps.push_back(&dofMap);   
572
573
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588

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


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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

609
    TEST_EXIT(nMacroElements >= MPI::COMM_WORLD.Get_size())
610
611
612
      ("The mesh has less macro elements than number of mpi processes!\n");
  }

613

614
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
615
  {
616
    FUNCNAME("MeshDistributor::synchVector()");
617
618
619
620
621
622
623
   
    int nLevels = levelData.getNumberOfLevels();
    for (int level = nLevels - 1; level >= 0; level--) {
      MPI::Intracomm &mpiComm = levelData.getMpiComm(level);
      
      if (mpiComm == MPI::COMM_SELF)
	continue;
624

625
      StdMpi<vector<double> > stdMpi(mpiComm);
626

627
628
629
630
631
632
633
634
635
636
637
638
639
640
      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);
641
      }
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
      
      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++];
	  }
670
	}
671
672
673
674
      }
    }
  }

675

676
  void MeshDistributor::fix3dMeshRefinement()
677
  {
678
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
679
680
681
682

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

683
684
685
686
    MPI::COMM_WORLD.Barrier();
    double first = MPI::Wtime();


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

690
691
692
693
694
695
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

707

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

711
    FixRefinementPatch::connectedEdges.clear();
712

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

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

721
722

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

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

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

739
740
741
742

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

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

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

792

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

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

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

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


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

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

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


850
851
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
852
853
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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

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

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


871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
  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()));
    }
  }


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

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


907
908
  void MeshDistributor::createMeshLevelStructure()
  {
909
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");
910

911
912
913
914
915
916
917
918
    int levelMode = -1;
    Parameters::get("parallel->level mode", levelMode);

    TEST_EXIT(levelMode >= 0)("Wrong level mode %d!\n", levelMode);

    if (levelMode == 0) {
      std::set<int> neighbours;
      levelData.init(neighbours);
919
      return;
920
    }
921

922
923
924
925
926
927
    if (levelMode == 1) {
      std::set<int> neighbours;
      levelData.init(neighbours);
      levelData.addLevelMode1();
      return;
    }
928

929
930
    TEST_EXIT(MPI::COMM_WORLD.Get_size() == 16)
      ("Only mpiSize == 16 supported!\n");
931

932
    if (levelMode == 2) {
933

934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
      std::set<int> neighbours;
      switch (mpiRank) {
      case 0:
	neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
	break;
      case 1:
	neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
	break;
      case 2:
	neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
	break;
      case 3:
	neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
	neighbours.insert(2); neighbours.insert(6);
	neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
	break;
      case 4:
	neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
	break;
      case 5:
	neighbours.insert(4); neighbours.insert(6); neighbours.insert