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

Thomas Witkowski's avatar
Thomas Witkowski committed
212
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
    // === In 3D we have to make some test, if the resulting mesh is valid.  ===
    // === If it is not valid, there is no possiblity yet to fix this        ===
    // === problem, just exit with an error message.                         ===
297

298
299
    check3dValidMesh();

300

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

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

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

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

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

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

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

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

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

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

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

346

Thomas Witkowski's avatar
Thomas Witkowski committed
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
372
  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");
    }
373
374

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


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


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

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

436

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

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


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

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

492
493
494

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

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

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

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

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

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

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

541
    problemStat.push_back(probStat);
542

543
544
545

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

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


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

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

    globalMeshDistributor->addProblemStat(probStat);
  }


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

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

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

578
  
579
  void MeshDistributor::testForMacroMesh()
580
  {
581
    FUNCNAME("MeshDistributor::testForMacroMesh()");
582
583
584
585
586
587
588

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
593
594
595
596
597
598
599
600
601
602
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

603

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
604
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
605
  {
606
607
    FUNCNAME("MeshDistributor::synchVector()");

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
610
611
612
613
614
615
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

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

    for (DofComm::Iterator it(dc.getSendDofs()); 
616
	 !it.end(); it.nextRank()) {
617
      vector<double> dofs;
618

619
620
      for (int i = 0; i < vec.getSize(); i++) {
	DOFVector<double> &dofVec = *(vec.getDOFVector(i));
Thomas Witkowski's avatar
Thomas Witkowski committed
621
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
622
623
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
624

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
625
626
627
628
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
629
    }
630
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
631
632
633
634
635
636
    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);
    }
637

638
    stdMpi.startCommunication();
639

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
640
641
642
643
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
644
645

      int counter = 0;
646
      vector<double> &recvData =  stdMpi.getRecvData(rank);
647
648
649
650

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

651
652
653
654
655
	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++];
	}
656
657
658
659
      }
    }
  }

660

661
662
663
664
665
666
667
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

668
669
670
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

671
672
673
674
675
676
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
677
      for (int i = 0; i < 6; i++) {
678
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
679
680
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
681
682
683
684
685
686
687
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

688

689
690
691
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

692
    FixRefinementPatch::connectedEdges.clear();
693
    bool valid3dMesh = true;
694

695
    // Traverse all edges
696
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
697
      // For each edge get all macro elements that contain this edge.
698
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
699
700
701
702

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

703
704
705
706
707
708
709

      // We have now to check, if the current edge connects two macro elements,
      // which are otherwise not connected. The basic idea to check this is very
      // simple: We take the first macro element in rank's subdomain that contain
      // this edge and add it to the set variable "el0". All other macro elements
      // which share this edge are added to "el1". 

710
      std::set<int> el0, el1;
711
      map<int, int> edgeNoInEl;
712
713
714

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
715
716
717
718
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
719
720
721
722

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
723
           
724
725
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

726
727
      // If el1 is empty, there is only on macro element in the mesh which
      // contains this edge. Hence, we can continue with checking another edge.
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
      if (el1.empty())
	continue;
      
      bool found = false;      
      do {
	found = false;
	for (std::set<int>::iterator elIt = el0.begin(); elIt != el0.end(); ++elIt) {
	  for (int i = 0; i < 4; i++) {
	    int neighEl = macroElementNeighbours[*elIt][i];
	    if (neighEl != -1 && el1.count(neighEl)) {
	      el0.insert(neighEl);
	      el1.erase(neighEl);
	      found = true;
	    }
	  }
	  
	  if (found)
	    break;
	}
      } while (found);
      
      if (!el1.empty()) {
750
751
752
753
	valid3dMesh = false;

	MSG_DBG("Unvalid 3D mesh with %d (%d - %d) elements on this edge!\n", 
		edgeEls.size(), el0.size(), el1.size());
754

755
756
757
	for (std::set<int>::iterator elIt0 = el0.begin(); elIt0 != el0.end(); ++elIt0) {
	  for (std::set<int>::iterator elIt1 = el1.begin(); elIt1 != el1.end(); ++elIt1) {
	    pair<Element*, int> edge0 = 
758
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
759
	    pair<Element*, int> edge1 = 
760
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
761

762
#if (DEBUG != 0)
763
764
765
766
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    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
	    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]);
#endif
779
780
781
782
783
784
785
786
787

	    TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n");

	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge0, edge1));
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}
      }
    }
788
789

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
790
791
792
  }


793
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
794
					   int level,
795
					   DofContainer& dofs)
796
797
798
799
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
800
801
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
802
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
803
804
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
805
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
806
807
808
809
810
811

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


812
813
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
814
815
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

816
    // Remove periodic boundaries in boundary manager on matrices and vectors.
817
818
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
819
820
821

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
822
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
823
824
825
826
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
827
828
829

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
830
831
832
  }


833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
  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()));
    }
  }


855
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
856
  {
857
858
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
859
860
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
861
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
862
	boundaryMap.erase(it++);
863
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
864
865
866
867
868
	++it;      
    }    
  }


869
870
871
872
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

873
874
875
    if (mpiSize != 16)
      return;

876
    std::set<int> neighbours;
877
878
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
879
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
880
881
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
882
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
883
884
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
885
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
886
887
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
888
889
890
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
891
892
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
893
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
894
895
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
896
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
897
898
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
899
900
901
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
902
903
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
904
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);
905
906
      break;
    case 8:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
907
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(9);  neighbours.insert(10); neighbours.insert(11);
908
909
      break;
    case 9:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
910
911
912
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(6); 
      neighbours.insert(8); neighbours.insert(12);
      neighbours.insert(10); neighbours.insert(11); neighbours.insert(14); 
913
914
      break;
    case 10:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
915
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(11);
916
917
      break;
    case 11:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
918
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12);  neighbours.insert(10); neighbours.insert(14);
919
920
      break;
    case 12:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
921
922
923
      neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); 
      neighbours.insert(9); neighbours.insert(13);
      neighbours.insert(11); neighbours.insert(14); neighbours.insert(15); 
924
925
      break;
    case 13:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
926
      neighbours.insert(6); neighbours.insert(7); neighbours.insert(12);  neighbours.insert(14); neighbours.insert(15);
927
928
      break;
    case 14:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
929
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13);  neighbours.insert(11); neighbours.insert(15);
930
931
      break;
    case 15:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
932
      neighbours.insert(12); neighbours.insert(13); neighbours.insert(14);
933
934
      break;
    }    
935