MeshDistributor.cc 63.4 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
      componentSpaces(0),
78
79
80
81
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
82
83
      dofMap(FESPACE_WISE),
      dofMapSd(FESPACE_WISE),
84
      deserialized(false),
85
      writeSerializationFile(false),
86
      repartitioningAllowed(false),
87
      repartitionIthChange(20),
88
89
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
90
      repartitioningFailed(0),
91
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
92
      lastMeshChangeIndex(0),
93
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
94
95
      boundaryDofInfo(1),
      meshAdaptivity(true)
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);
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
106
    Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
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));

132
    TEST_EXIT(partitioner)("Could not create partitioner \"%s\"!\n", partStr.c_str());
133
134
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
135

Thomas Witkowski's avatar
Thomas Witkowski committed
136
137
138
139
140
141
142
  MeshDistributor::~MeshDistributor()
  {
    if (partitioner)
      delete partitioner;
  }


143
  void MeshDistributor::initParallelization()
144
  {
145
    FUNCNAME("MeshDistributor::initParallelization()");
146

147
148
149
    if (initialized)
      return;

150
151
152
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");

153
    TEST_EXIT(componentSpaces.size() > 0)
154
      ("No FE space has been defined for the mesh distributor!\n");
155
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
156

157
158
159
160
161
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

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

174
    elObjDb.setFeSpace(uniqueFeSpaces[0]);
175

Thomas Witkowski's avatar
Thomas Witkowski committed
176
177
178
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

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

187
      removePeriodicBoundaryConditions();
188

189
      elObjDb.createMacroElementInfo(allMacroElements);
190

Thomas Witkowski's avatar
Thomas Witkowski committed
191
      updateLocalGlobalNumbering();
192

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

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

      initialized = true;
202
      return;
203
    }
204

205

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

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

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

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

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

226

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

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

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

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

249
250
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

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

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

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

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

    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))) {
274
275
276

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

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

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

284
285
286
	}
      }
    }
287

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

    removeMacroElements();
291

292
293
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
294

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

302
    updateLocalGlobalNumbering();
303

304
305
306
    // === 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.                         ===
307

308
309
    check3dValidMesh();

310

311
312
    // === If in debug mode, make some tests. ===

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

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

321
    debug::writeMesh(uniqueFeSpaces[0], -1, debugOutputDir + "macro_mesh");   
322
323

    MSG("Debug mode tests finished!\n");
324
#endif
325

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

329
330
331
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

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

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

344
      updateLocalGlobalNumbering();
345

346
347
      // === Update periodic mapping, if there are periodic boundaries. ===     

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

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
354
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
355
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
356
357
      elObjDb.clear();

358
    initialized = true;
359
360
  }

361

362
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
363
  {
364
    FUNCNAME("MeshDistributor::addProblemStat()");
365

366
367
368
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

369
370
    TEST_EXIT(componentSpaces.size() == 0)
      ("Parallelization of coupled problems is deactived at the moment!\n");
371

372
373
374
    componentSpaces = probStat->getComponentFeSpaces();
    uniqueFeSpaces = probStat->getFeSpaces();
    mesh = uniqueFeSpaces[0]->getMesh();
375
376
377
378
379
380
381
382
383
384
385
    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());
386
    }
387
388
389
    
    partitioner->setMesh(mesh);
    
390
391
392

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
393
    Parameters::get(probStat->getName() + "->output->write serialization",  
394
		    writeSerialization);
395
    if (writeSerialization && !writeSerializationFile) {
396
      string filename = "";
397
      Parameters::get(name + "->output->serialization filename", filename);
398
399
400
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
401
402

      int tsModulo = -1;
403
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
404
		      tsModulo);
405
      
406
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
407
408
      writeSerializationFile = true;
    }    
409
410

    int readSerialization = 0;
411
    Parameters::get(probStat->getName() + "->input->read serialization", 
412
		    readSerialization);
413
    if (readSerialization) {
414
      string filename = "";
415
      Parameters::get(probStat->getName() + "->input->serialization filename", 
416
		      filename);
417
      filename += ".p" + lexical_cast<string>(mpiRank);
418
      MSG("Start deserialization with %s\n", filename.c_str());
419
      ifstream in(filename.c_str());
420
421
422
423

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

424
425
426
427
428
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

429
      probStat->deserialize(in);
430
      in.close();
431
432
      MSG("Deserialization from file: %s\n", filename.c_str());

433
      filename = "";
434
      Parameters::get(name + "->input->serialization filename", filename);
435
436
437
438
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
439
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
440
441
442
443
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
444
445
446
447
448

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
449
450
451
452
453
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
454
455
    }

456
    problemStat.push_back(probStat);
457

458
459
460
    // 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.
461

462
    if (initialized)
463
      removePeriodicBoundaryConditions(probStat);
464
465
466
  }


467
468
469
470
471
472
473
474
475
476
477
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


478
  void MeshDistributor::exitParallelization()
479
  {}
480

481
  
482
  void MeshDistributor::testForMacroMesh()
483
  {
484
    FUNCNAME("MeshDistributor::testForMacroMesh()");
485
486
487
488
489
490
491

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
496
497
498
499
500
501
502
503
504
505
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

506

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
507
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
508
  {
509
510
    FUNCNAME("MeshDistributor::synchVector()");

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
513
514
515
516
517
518
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

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

    for (DofComm::Iterator it(dc.getSendDofs()); 
519
	 !it.end(); it.nextRank()) {
520
      vector<double> dofs;
521

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

Thomas Witkowski's avatar
Thomas Witkowski committed
525
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
526
527
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
528

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
529
530
531
532
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
533
    }
534
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
535
536
537
538
539
540
    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);
    }
541

542
    stdMpi.startCommunication();
543

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
544
545
546
547
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
548
549

      int counter = 0;
550
      vector<double> &recvData =  stdMpi.getRecvData(rank);
551
552
553
554

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

555
556
557
558
559
	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++];
	}
560
561
562
563
      }
    }
  }

564

565
566
567
568
569
570
571
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

572
573
574
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

575
576
577
578
579
580
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
581
      for (int i = 0; i < 6; i++) {
582
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
583
584
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
585
586
587
588
589
590
591
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

592

593
594
595
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

596
    FixRefinementPatch::connectedEdges.clear();
597
    bool valid3dMesh = true;
598

599
    // Traverse all edges
600
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
601
      // For each edge get all macro elements that contain this edge.
602
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
603
604
605
606

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

607
608
609
610
611
612
613

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

614
      std::set<int> el0, el1;
615
      map<int, int> edgeNoInEl;
616
617
618

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
619
620
621
622
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
623
624
625
626

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
627
           
628
629
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

630
631
      // 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.
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
      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()) {
654
655
656
657
	valid3dMesh = false;

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

659
660
661
	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 = 
662
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
663
	    pair<Element*, int> edge1 = 
664
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
665

666
#if (DEBUG != 0)
667
668
669
670
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
671
672
	    mesh->getDofIndexCoords(dofEdge0.first, uniqueFeSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, uniqueFeSpaces[0], c1);
673

674
675
676
677
678
679
680
681
682
	    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
683
684
685
686
687
688
689
690
691

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

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
694
695
696
  }


697
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
698
					   int level,
699
					   DofContainer& dofs)
700
701
702
703
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
704
705
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
706
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
707
708
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
709
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
710
711
712
713
714
715

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


716
717
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
718
719
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

720
    // Remove periodic boundaries in boundary manager on matrices and vectors.
721
722
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
723
724
725

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
726
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
727
728
729
730
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
731
732
733

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
734
735
736
  }


737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
  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()));
    }
  }


759
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
760
  {
761
762
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
763
764
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
765
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
766
	boundaryMap.erase(it++);
767
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
768
769
770
771
772
	++it;      
    }    
  }


773
774
775
776
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

777
778
779
    if (mpiSize != 16)
      return;

780
    std::set<int> neighbours;
781
782
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
783
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
784
785
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
786
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
787
788
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
789
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
790
791
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
792
793
794
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
795
796
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
797
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
798
799
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
800
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
801
802
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
803
804
805
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
806
807
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
808
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);
809
810
      break;
    case 8:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
811
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(9);  neighbours.insert(10); neighbours.insert(11);
812
813
      break;
    case 9:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
814
815
816
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(6); 
      neighbours.insert(8); neighbours.insert(12);
      neighbours.insert(10); neighbours.insert(11); neighbours.insert(14); 
817
818
      break;
    case 10:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
819
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(11);
820
821
      break;
    case 11:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
822
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12);  neighbours.insert(10); neighbours.insert(14);
823
824
      break;
    case 12:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
825
826
827
      neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); 
      neighbours.insert(9); neighbours.insert(13);
      neighbours.insert(11); neighbours.insert(14); neighbours.insert(15); 
828
829
      break;
    case 13:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
830
      neighbours.insert(6); neighbours.insert(7); neighbours.insert(12);  neighbours.insert(14); neighbours.insert(15);
831
832
      break;
    case 14:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
833
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13);  neighbours.insert(11); neighbours.insert(15);
834
835
      break;
    case 15:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
836
      neighbours.insert(12); neighbours.insert(13); neighbours.insert(14);
837
838
      break;
    }    
839

840
    TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
841
842
843
844
845
846
847

    levelData.init(neighbours);

    bool multiLevelTest = false;
    Parameters::get("parallel->multi level test", multiLevelTest);
    if (multiLevelTest) {
      switch (mpiRank) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
848
      case 0: case 1: case 2: case 3:
849
850
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
851
	  m.insert(0); m.insert(1); m.insert(2); m.insert(3);
852
	  levelData.addLevel(m, 0);
853
854
	}
	break;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
855
      case 4: case 5: case 6: case 7:
856
857
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
858
	  m.insert(4); m.insert(5); m.insert(6); m.insert(7);
859
	  levelData.addLevel(m, 1);
860
861
	}
	break;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed