MeshDistributor.cc 64.9 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),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
92
93
      boundaryDofInfo(1),
      meshAdaptivity(true)
94
  {
95
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
96

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

101
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
102
103
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
104
    Parameters::get(name + "->log main rank", Msg::outputMainRank);
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
105
    Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
106

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
134

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


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

146
147
148
    if (initialized)
      return;

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

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

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

    bool doNext = false;
    do {
161
162
      doNext = false;
      for (unsigned int i = 0; i < feSpaces.size() - 1; i++) {
163
164
165
166
167
168
169
	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;
	}
170
      }
171
172
    } while (doNext);

173
    elObjDb.setFeSpace(feSpaces[0]);
174

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

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

186
      removePeriodicBoundaryConditions();
187

188
      elObjDb.createMacroElementInfo(allMacroElements);
189

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

192
193
      setRankDofs();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
196
#if (DEBUG != 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
197
198
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], dofMap, 
				    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(feSpaces[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

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

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

361
    initialized = true;
362
363
  }

364

365
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
366
  {
367
    FUNCNAME("MeshDistributor::addProblemStat()");
368

369
370
371
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

372
373
374
375
376
377
378
379
380
381
    // === 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");
382
383
      }

384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
      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());
402
    }
403
404
405
    
    partitioner->setMesh(mesh);
    
406
407
408

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
409
    Parameters::get(probStat->getName() + "->output->write serialization",  
410
		    writeSerialization);
411
    if (writeSerialization && !writeSerializationFile) {
412
      string filename = "";
413
      Parameters::get(name + "->output->serialization filename", filename);
414
415
416
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
417
418

      int tsModulo = -1;
419
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
420
		      tsModulo);
421
      
422
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
423
424
      writeSerializationFile = true;
    }    
425
426

    int readSerialization = 0;
427
    Parameters::get(probStat->getName() + "->input->read serialization", 
428
		    readSerialization);
429
    if (readSerialization) {
430
      string filename = "";
431
      Parameters::get(probStat->getName() + "->input->serialization filename", 
432
		      filename);
433
      filename += ".p" + lexical_cast<string>(mpiRank);
434
      MSG("Start deserialization with %s\n", filename.c_str());
435
      ifstream in(filename.c_str());
436
437
438
439

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

440
441
442
443
444
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

445
      probStat->deserialize(in);
446
      in.close();
447
448
      MSG("Deserialization from file: %s\n", filename.c_str());

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

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
465
466
467
468
469
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
470
471
    }

472
    problemStat.push_back(probStat);
473

474
475
476
    // 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.
477

478
    if (initialized) {
479
      setRankDofs(probStat);
480
481
      removePeriodicBoundaryConditions(probStat);
    }
482
483
484
  }


485
486
487
488
489
490
491
492
493
494
495
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


496
  void MeshDistributor::exitParallelization()
497
  {}
498

499
  
500
  void MeshDistributor::testForMacroMesh()
501
  {
502
    FUNCNAME("MeshDistributor::testForMacroMesh()");
503
504
505
506
507
508
509

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
514
515
516
517
518
519
520
521
522
523
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

524

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
525
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
526
  {
527
528
    FUNCNAME("MeshDistributor::synchVector()");

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
531
532
533
534
535
536
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

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

    for (DofComm::Iterator it(dc.getSendDofs()); 
537
	 !it.end(); it.nextRank()) {
538
      vector<double> dofs;
539

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

Thomas Witkowski's avatar
Thomas Witkowski committed
543
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
544
545
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
546

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
547
548
549
550
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
551
    }
552
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
553
554
555
556
557
558
    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);
    }
559

560
    stdMpi.startCommunication();
561

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
562
563
564
565
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
566
567

      int counter = 0;
568
      vector<double> &recvData =  stdMpi.getRecvData(rank);
569
570
571
572

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

573
574
575
576
577
	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++];
	}
578
579
580
581
      }
    }
  }

582

583
584
585
586
587
588
589
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

590
591
592
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

593
594
595
596
597
598
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
599
      for (int i = 0; i < 6; i++) {
600
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
601
602
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
603
604
605
606
607
608
609
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

610

611
612
613
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

614
    FixRefinementPatch::connectedEdges.clear();
615
    bool valid3dMesh = true;
616

617
    // Traverse all edges
618
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
619
      // For each edge get all macro elements that contain this edge.
620
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
621
622
623
624

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

625
626
627
628
629
630
631

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

632
      std::set<int> el0, el1;
633
      map<int, int> edgeNoInEl;
634
635
636

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
637
638
639
640
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
641
642
643
644

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
645
           
646
647
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

648
649
      // 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.
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
      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()) {
672
673
674
675
	valid3dMesh = false;

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

677
678
679
	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 = 
680
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
681
	    pair<Element*, int> edge1 = 
682
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
683

684
#if (DEBUG != 0)
685
686
687
688
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
689
690
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
691

692
693
694
695
696
697
698
699
700
	    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
701
702
703
704
705
706
707
708
709

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

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
712
713
714
  }


715
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
716
					   int level,
717
					   DofContainer& dofs)
718
719
720
721
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
722
723
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
724
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
725
726
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
727
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
728
729
730
731
732
733

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


734
  void MeshDistributor::setRankDofs(ProblemStatSeq *probStat)
735
  {
736
737
    FUNCNAME("MeshDistributor::setRankDofs()");

738
    int nComponents = probStat->getNumComponents();
739
740
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
741
742
743
744
745
746
747
	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
748
	  probStat->getSystemMatrix(i, j)->setDofMap(dofMap[rowFeSpace]);
749
750
	}
    
751

752
      TEST_EXIT_DBG(probStat->getRhs()->getDOFVector(i))("No RHS vector!\n");
753
754
      TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))
	("No solution vector!\n");
755
756
757
      TEST_EXIT_DBG(probStat->getRhsVector(i)->getFeSpace() ==
		    probStat->getSolution(i)->getFeSpace())
	("Should not happen!\n");
758
      
759
760
761
762
      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
763
764
      probStat->getRhsVector(i)->setDofMap(dofMap[feSpace]);
      probStat->getSolution(i)->setDofMap(dofMap[feSpace]);
765
766
767
768
    }    
  }


769
770
  void MeshDistributor::setRankDofs()
  {
771
772
    for (unsigned int i = 0; i < problemStat.size(); i++) 
      setRankDofs(problemStat[i]);
773
774
775
  }


776
777
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
778
779
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

780
    // Remove periodic boundaries in boundary manager on matrices and vectors.
781
782
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
783
784
785

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
786
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
787
788
789
790
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
791
792
793

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
794
795
796
  }


797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
  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()));
    }
  }


819
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
820
  {
821
822
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
823
824
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
825
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
826
	boundaryMap.erase(it++);
827
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
828
829
830
831
832
	++it;      
    }    
  }


833
834
835
836
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

837
838
839
    if (mpiSize != 16)
      return;

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