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


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

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

54
55
namespace AMDiS {

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

60
61
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

62
63
64
  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
65

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

71
  bool MeshDistributor::sebastianMode = false;
72

73
  MeshDistributor::MeshDistributor()
74
    : problemStat(0),
75
      initialized(false),
76
      name("parallel"),
77
      feSpaces(0),
78
79
80
81
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
82
      deserialized(false),
83
      writeSerializationFile(false),
84
      repartitioningAllowed(false),
85
      repartitionIthChange(20),
86
87
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
88
      repartitioningFailed(0),
89
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
90
      lastMeshChangeIndex(0),
91
      createBoundaryDofFlag(0),
92
      boundaryDofInfo(1)
93
  {
94
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
95

96
    mpiComm = MPI::COMM_WORLD;
97
98
    mpiRank = mpiComm.Get_rank();
    mpiSize = mpiComm.Get_size();
99

100
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
101
102
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
103
    Parameters::get(name + "->log main rank", Msg::outputMainRank);
104

105
    string partStr = "parmetis";
106
    Parameters::get(name + "->partitioner", partStr);
107
108
109
110

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

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

119
120
121
    if (partStr == "checker")
      partitioner = new CheckerPartitioner(&mpiComm);

Thomas Witkowski's avatar
Thomas Witkowski committed
122
123
124
    if (partStr == "simple")
      partitioner = new SimplePartitioner(&mpiComm);

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

129
    TEST_EXIT(partitioner)("Could not create partitioner \"%s\"!\n", partStr.c_str());
130
131
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
132

Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
135
136
137
138
139
  MeshDistributor::~MeshDistributor()
  {
    if (partitioner)
      delete partitioner;
  }


140
  void MeshDistributor::initParallelization()
141
  {
142
    FUNCNAME("MeshDistributor::initParallelization()");
143

144
145
146
    if (initialized)
      return;

147
148
149
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");

150
151
    TEST_EXIT(feSpaces.size() > 0)
      ("No FE space has been defined for the mesh distributor!\n");
152
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
153

154
155
156
157
158
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
159
160
      doNext = false;
      for (unsigned int i = 0; i < feSpaces.size() - 1; i++) {
161
162
163
164
165
166
167
	if (feSpaces[i]->getBasisFcts()->getDegree() >
	    feSpaces[i + 1]->getBasisFcts()->getDegree()) {
	  const FiniteElemSpace *tmp = feSpaces[i + 1];
	  feSpaces[i + 1] = feSpaces[i];
	  feSpaces[i] = tmp;
	  doNext = true;
	}
168
      }
169
170
    } while (doNext);

171
    elObjDb.setFeSpace(feSpaces[0]);
172

Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
175
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

176
177
178
    // 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).
179
    if (deserialized) {
Thomas Witkowski's avatar
Thomas Witkowski committed
180
181
      createMeshLevelStructure();
      
182
183
      updateMacroElementInfo();

184
      removePeriodicBoundaryConditions();
185

186
      elObjDb.createMacroElementInfo(allMacroElements);
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188
      updateLocalGlobalNumbering();
189

190
191
      setRankDofs();

Thomas Witkowski's avatar
Thomas Witkowski committed
192
193
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
194
#if (DEBUG != 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
195
196
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], dofMap, 
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
197
#endif    
198
199

      initialized = true;
200
      return;
201
    }
202

203

204
205
206
    // 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.
207
208
    testForMacroMesh();

209
210
    // For later mesh repartitioning, we need to store some information about
    // the macro mesh.
211
212
    createMacroElementInfo();

213
    // create an initial partitioning of the mesh
214
    partitioner->createInitialPartitioning();
215

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

    // and now partition the mesh    
220
221
    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
222
    partitioner->createPartitionMap(partitionMap);
223

224

225
#if (DEBUG != 0)    
226
227
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
228
229
#endif

230
    if (mpiRank == 0) {
231
#if (DEBUG != 0)    
232
      int writePartMesh = 1;
233
234
235
#else
      int writePartMesh = 0;
#endif
236
      Parameters::get("dbg->write part mesh", writePartMesh);
237

238
      if (writePartMesh > 0) {
239
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex");
240
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
241
      }
242
    }
243

244
    // Create interior boundary information.
245
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
246

247
248
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

249
250
251
252
253
    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))) {
254
255
256

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

257
258
259
260
261
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
262
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
263
264
265
266
267
268
269
270
271
	}
      }
    }

    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))) {
272
273
274

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

275
276
277
278
279
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

282
283
284
	}
      }
    }
285

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
286
287
288
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
289

290
291
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
292

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

300
    updateLocalGlobalNumbering();
301

302
303
304
    // === 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.                         ===
305

306
307
    check3dValidMesh();

308

309
310
    // === If in debug mode, make some tests. ===

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

314
    ParallelDebug::testAllElements(*this);
315
    debug::testSortedDofs(mesh, elMap);
316
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
317
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
318

319
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
320
321

    MSG("Debug mode tests finished!\n");
322
#endif
323

324
    // Create periodic DOF mapping, if there are periodic boundaries.
325
    createPeriodicMap();
326

327
328
329
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

330
331
332
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
333

334
    // === Global refinements. ===
335
    
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    int globalRefinement = 0;
337
    Parameters::get(mesh->getName() + "->global refinements", globalRefinement);
338
    
Thomas Witkowski's avatar
Thomas Witkowski committed
339
    if (globalRefinement > 0) {
340
      refineManager->globalRefine(mesh, globalRefinement);
341

342
      updateLocalGlobalNumbering();
343

344
345
      // === Update periodic mapping, if there are periodic boundaries. ===     

346
      createPeriodicMap();
347
348
349
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
350
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
351

352
    // Set DOF rank information to all matrices and vectors.
353
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
354

355
    initialized = true;
356
357
  }

358

359
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
360
  {
361
    FUNCNAME("MeshDistributor::addProblemStat()");
362

363
364
365
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

366
367
368
369
370
371
372
373
374
375
    // === Add FE spaces from stationary problem to mesh distributor. ===

    for (unsigned int i = 0; i < probStat->getFeSpaces().size(); i++) {
      bool foundFeSpace = false;
      for (unsigned int j = 0; j < feSpaces.size(); j++) {
	if (feSpaces[j] == probStat->getFeSpaces()[i])
	  foundFeSpace = true;

	TEST_EXIT(feSpaces[j]->getMesh() == probStat->getFeSpaces()[i]->getMesh())
	  ("MeshDistributor does not yet support usage of different meshes!\n");
376
377
      }

378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
      if (!foundFeSpace)
	feSpaces.push_back(probStat->getFeSpaces()[i]);
    }

    TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n");

    mesh = feSpaces[0]->getMesh();
    info = probStat->getInfo();
    
    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());
396
    }
397
398
399
    
    partitioner->setMesh(mesh);
    
400
401
402

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
403
    Parameters::get(probStat->getName() + "->output->write serialization",  
404
		    writeSerialization);
405
    if (writeSerialization && !writeSerializationFile) {
406
      string filename = "";
407
      Parameters::get(name + "->output->serialization filename", filename);
408
409
410
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
411
412

      int tsModulo = -1;
413
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
414
		      tsModulo);
415
      
416
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
417
418
      writeSerializationFile = true;
    }    
419
420

    int readSerialization = 0;
421
    Parameters::get(probStat->getName() + "->input->read serialization", 
422
		    readSerialization);
423
    if (readSerialization) {
424
      string filename = "";
425
      Parameters::get(probStat->getName() + "->input->serialization filename", 
426
		      filename);
427
      filename += ".p" + lexical_cast<string>(mpiRank);
428
      MSG("Start deserialization with %s\n", filename.c_str());
429
      ifstream in(filename.c_str());
430
431
432
433

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

434
435
436
437
438
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

439
      probStat->deserialize(in);
440
      in.close();
441
442
      MSG("Deserialization from file: %s\n", filename.c_str());

443
      filename = "";
444
      Parameters::get(name + "->input->serialization filename", filename);
445
446
447
448
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
449
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
450
451
452
453
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
454
455
456
457
458

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
459
460
461
462
463
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
464
465
    }

466
    problemStat.push_back(probStat);
467

468
469
470
    // If the mesh distributor is already initialized, don't forget to set rank
    // DOFs object to the matrices and vectors of the added stationary problem,
    // and to remove the periodic boundary conditions on these objects.
471

472
    if (initialized) {
473
      setRankDofs(probStat);
474
475
      removePeriodicBoundaryConditions(probStat);
    }
476
477
478
  }


479
480
481
482
483
484
485
486
487
488
489
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


490
  void MeshDistributor::exitParallelization()
491
  {}
492

493
  
494
  void MeshDistributor::testForMacroMesh()
495
  {
496
    FUNCNAME("MeshDistributor::testForMacroMesh()");
497
498
499
500
501
502
503

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
508
509
510
511
512
513
514
515
516
517
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

518

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
519
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
520
  {
521
522
    FUNCNAME("MeshDistributor::synchVector()");

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
525
526
527
528
529
530
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

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

    for (DofComm::Iterator it(dc.getSendDofs()); 
531
	 !it.end(); it.nextRank()) {
532
      vector<double> dofs;
533

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

Thomas Witkowski's avatar
Thomas Witkowski committed
537
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
538
539
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
540

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
541
542
543
544
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
545
    }
546
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
547
548
549
550
551
552
    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);
    }
553

554
    stdMpi.startCommunication();
555

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
556
557
558
559
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
560
561

      int counter = 0;
562
      vector<double> &recvData =  stdMpi.getRecvData(rank);
563
564
565
566

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

567
568
569
570
571
	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++];
	}
572
573
574
575
      }
    }
  }

576

577
578
579
580
581
582
583
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

584
585
586
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

587
588
589
590
591
592
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
593
      for (int i = 0; i < 6; i++) {
594
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
595
596
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
597
598
599
600
601
602
603
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

604

605
606
607
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

608
    FixRefinementPatch::connectedEdges.clear();
609
    bool valid3dMesh = true;
610

611
    // Traverse all edges
612
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
613
      // For each edge get all macro elements that contain this edge.
614
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
615
616
617
618

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

619
620
621
622
623
624
625

      // 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". 

626
      std::set<int> el0, el1;
627
      map<int, int> edgeNoInEl;
628
629
630

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
631
632
633
634
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
635
636
637
638

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
639
           
640
641
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

642
643
      // 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.
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
      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()) {
666
667
668
669
	valid3dMesh = false;

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

671
672
673
	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 = 
674
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
675
	    pair<Element*, int> edge1 = 
676
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
677

678
#if (DEBUG != 0)
679
680
681
682
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
683
684
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
685

686
687
688
689
690
691
692
693
694
	    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
695
696
697
698
699
700
701
702
703

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

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
706
707
708
  }


709
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
710
					   int level,
711
					   DofContainer& dofs)
712
713
714
715
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
716
717
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
718
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
719
720
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
721
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
722
723
724
725
726
727

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


728
  void MeshDistributor::setRankDofs(ProblemStatSeq *probStat)
729
  {
730
731
    FUNCNAME("MeshDistributor::setRankDofs()");

732
    int nComponents = probStat->getNumComponents();
733
734
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
735
736
737
738
739
740
741
	if (probStat->getSystemMatrix(i, j)) {
	  const FiniteElemSpace *rowFeSpace = 
	    probStat->getSystemMatrix(i, j)->getRowFeSpace();

	  TEST_EXIT(find(feSpaces.begin(), feSpaces.end(), rowFeSpace) != feSpaces.end())
	    ("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
742
	  probStat->getSystemMatrix(i, j)->setDofMap(dofMap[rowFeSpace]);
743
744
	}
    
745

746
      TEST_EXIT_DBG(probStat->getRhs()->getDOFVector(i))("No RHS vector!\n");
747
748
      TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))
	("No solution vector!\n");
749
750
751
      TEST_EXIT_DBG(probStat->getRhsVector(i)->getFeSpace() ==
		    probStat->getSolution(i)->getFeSpace())
	("Should not happen!\n");
752
      
753
754
755
756
      const FiniteElemSpace *feSpace = probStat->getRhsVector(i)->getFeSpace();
      TEST_EXIT(find(feSpaces.begin(), feSpaces.end(), feSpace) != feSpaces.end())
	("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
757
758
      probStat->getRhsVector(i)->setDofMap(dofMap[feSpace]);
      probStat->getSolution(i)->setDofMap(dofMap[feSpace]);
759
760
761
762
    }    
  }


763
764
  void MeshDistributor::setRankDofs()
  {
765
766
    for (unsigned int i = 0; i < problemStat.size(); i++) 
      setRankDofs(problemStat[i]);
767
768
769
  }


770
771
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
772
773
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

774
    // Remove periodic boundaries in boundary manager on matrices and vectors.
775
776
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
777
778
779

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
780
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
781
782
783
784
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
785
786
787

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
788
789
790
  }


791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
  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()));
    }
  }


813
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
814
  {
815
816
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
817
818
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
819
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
820
	boundaryMap.erase(it++);
821
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
822
823
824
825
826
	++it;      
    }    
  }


827
828
829
830
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

831
832
833
    if (mpiSize != 16)
      return;

834
    std::set<int> neighbours;
835
836
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
837
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
838
839
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
840
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
841
842
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
843
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
844
845
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
846
847
848
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
849
850
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
851
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
852
853
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
854
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
855
856
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
857
858
859
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
860
861
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
862
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);
863
864
      break;
    case 8:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
865
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(9);  neighbours.insert(10); neighbours.insert(11);
866
867
      break;
    case 9:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
868
869
870
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(6); 
      neighbours.insert(8); neighbours.insert(12);
      neighbours.insert(10); neighbours.insert(11); neighbours.insert(14); 
871
872
      break;
    case 10:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
873
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(11);
874
875
      break;
    case 11:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
876
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12);  neighbours.insert(10); neighbours.insert(14);
877
878
      break;
    case 12:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
879
880
881
      neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); 
      neighbours.insert(9); neighbours.insert(13);
      neighbours.insert(11); neighbours.insert(14); neighbours.insert(15); 
Thomas Witkowski's avatar