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

133
  void MeshDistributor::initParallelization()
134
  {
135
    FUNCNAME("MeshDistributor::initParallelization()");
136

137
138
139
    if (initialized)
      return;

140
141
142
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");

143
144
    TEST_EXIT(feSpaces.size() > 0)
      ("No FE space has been defined for the mesh distributor!\n");
145
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
146

147
148
149
150
151
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
152
153
      doNext = false;
      for (unsigned int i = 0; i < feSpaces.size() - 1; i++) {
154
155
156
157
158
159
160
	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;
	}
161
      }
162
163
    } while (doNext);

164
    elObjDb.setFeSpace(feSpaces[0]);
165

Thomas Witkowski's avatar
Thomas Witkowski committed
166
167
168
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

169
170
171
    // 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).
172
    if (deserialized) {
Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
      createMeshLevelStructure();
      
175
176
      updateMacroElementInfo();

177
      removePeriodicBoundaryConditions();
178

179
      elObjDb.createMacroElementInfo(allMacroElements);
180

Thomas Witkowski's avatar
Thomas Witkowski committed
181
      updateLocalGlobalNumbering();
182

183
184
      setRankDofs();

Thomas Witkowski's avatar
Thomas Witkowski committed
185
186
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
187
#if (DEBUG != 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
188
189
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], dofMap, 
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
190
#endif    
191
192

      initialized = true;
193
      return;
194
    }
195

196

197
198
199
    // 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.
200
201
    testForMacroMesh();

202
203
    // For later mesh repartitioning, we need to store some information about
    // the macro mesh.
204
205
    createMacroElementInfo();

206
    // create an initial partitioning of the mesh
207
    partitioner->createInitialPartitioning();
208

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

    // and now partition the mesh    
213
214
    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
215
    partitioner->createPartitionMap(partitionMap);
216

217

218
#if (DEBUG != 0)    
219
220
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
221
222
#endif

223
    if (mpiRank == 0) {
224
#if (DEBUG != 0)    
225
      int writePartMesh = 1;
226
227
228
#else
      int writePartMesh = 0;
#endif
229
      Parameters::get("dbg->write part mesh", writePartMesh);
230

231
      if (writePartMesh > 0) {
232
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex");
233
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
234
      }
235
    }
236

237
    // Create interior boundary information.
238
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
239

240
241
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

242
243
244
245
246
    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))) {
247
248
249

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

250
251
252
253
254
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

    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))) {
265
266
267

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

268
269
270
271
272
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

275
276
277
	}
      }
    }
278

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
279
280
281
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
282

283
284
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
285

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

293
    updateLocalGlobalNumbering();
294

295
296
297
    // === 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.                         ===
298

299
300
    check3dValidMesh();

301

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

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

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

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

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

317
    // Create periodic DOF mapping, if there are periodic boundaries.
318
    createPeriodicMap();
319

320
321
322
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

323
324
325
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
326

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

335
      updateLocalGlobalNumbering();
336

337
338
      // === Update periodic mapping, if there are periodic boundaries. ===     

339
      createPeriodicMap();
340
341
342
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
343
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
344

345
    // Set DOF rank information to all matrices and vectors.
346
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
347

348
    initialized = true;
349
350
  }

351

352
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
353
  {
354
    FUNCNAME("MeshDistributor::addProblemStat()");
355

356
357
358
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

359
360
361
362
363
364
365
366
367
368
    // === 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");
369
370
      }

371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
      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());
389
    }
390
391
392
    
    partitioner->setMesh(mesh);
    
393
394
395

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

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

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

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

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

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

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

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

459
    problemStat.push_back(probStat);
460

461
462
463
    // 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.
464

465
    if (initialized) {
466
      setRankDofs(probStat);
467
468
      removePeriodicBoundaryConditions(probStat);
    }
469
470
471
  }


472
473
474
475
476
477
478
479
480
481
482
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


483
  void MeshDistributor::exitParallelization()
484
  {}
485

486
  
487
  void MeshDistributor::testForMacroMesh()
488
  {
489
    FUNCNAME("MeshDistributor::testForMacroMesh()");
490
491
492
493
494
495
496

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
501
502
503
504
505
506
507
508
509
510
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

511

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
512
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
513
  {
514
515
    FUNCNAME("MeshDistributor::synchVector()");

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
518
519
520
521
522
523
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

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

    for (DofComm::Iterator it(dc.getSendDofs()); 
524
	 !it.end(); it.nextRank()) {
525
      vector<double> dofs;
526

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

Thomas Witkowski's avatar
Thomas Witkowski committed
530
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
531
532
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
533

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

547
    stdMpi.startCommunication();
548

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
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);
553
554

      int counter = 0;
555
      vector<double> &recvData =  stdMpi.getRecvData(rank);
556
557
558
559

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

560
561
562
563
564
	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++];
	}
565
566
567
568
      }
    }
  }

569

570
571
572
573
574
575
576
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

577
578
579
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

580
581
582
583
584
585
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

597

598
599
600
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

601
    FixRefinementPatch::connectedEdges.clear();
602
    bool valid3dMesh = true;
603

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

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

612
613
614
615
616
617
618

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

619
      std::set<int> el0, el1;
620
      map<int, int> edgeNoInEl;
621
622
623

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
624
625
626
627
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
628
629
630
631

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
632
           
633
634
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

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

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

664
665
666
	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 = 
667
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
668
	    pair<Element*, int> edge1 = 
669
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
670

671
#if (DEBUG != 0)
672
673
674
675
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
676
677
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
678

679
680
681
682
683
684
685
686
687
	    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
688
689
690
691
692
693
694
695
696

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

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
699
700
701
  }


702
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
703
					   int level,
704
					   DofContainer& dofs)
705
706
707
708
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
709
710
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
711
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
712
713
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
714
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
715
716
717
718
719
720

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


721
  void MeshDistributor::setRankDofs(ProblemStatSeq *probStat)
722
  {
723
724
    FUNCNAME("MeshDistributor::setRankDofs()");

725
    int nComponents = probStat->getNumComponents();
726
727
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
728
729
730
731
732
733
734
	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
735
	  probStat->getSystemMatrix(i, j)->setDofMap(dofMap[rowFeSpace]);
736
737
	}
    
738

739
      TEST_EXIT_DBG(probStat->getRhs()->getDOFVector(i))("No RHS vector!\n");
740
741
      TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))
	("No solution vector!\n");
742
743
744
      TEST_EXIT_DBG(probStat->getRhsVector(i)->getFeSpace() ==
		    probStat->getSolution(i)->getFeSpace())
	("Should not happen!\n");
745
      
746
747
748
749
      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
750
751
      probStat->getRhsVector(i)->setDofMap(dofMap[feSpace]);
      probStat->getSolution(i)->setDofMap(dofMap[feSpace]);
752
753
754
755
    }    
  }


756
757
  void MeshDistributor::setRankDofs()
  {
758
759
    for (unsigned int i = 0; i < problemStat.size(); i++) 
      setRankDofs(problemStat[i]);
760
761
762
  }


763
764
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
765
766
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

767
    // Remove periodic boundaries in boundary manager on matrices and vectors.
768
769
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
770
771
772

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
773
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
774
775
776
777
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
778
779
780

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
781
782
783
  }


784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
  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()));
    }
  }


806
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
807
  {
808
809
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
810
811
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
812
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
813
	boundaryMap.erase(it++);
814
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
815
816
817
818
819
	++it;      
    }    
  }


820
821
822
823
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

824
825
826
    if (mpiSize != 16)
      return;

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