MeshDistributor.cc 64.1 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
  MeshDistributor::MeshDistributor()
72
    : problemStat(0),
73
      initialized(false),
74
      name("parallel"),
75
76
77
      mesh(NULL),
      refineManager(NULL),
      partitioner(NULL),
78
      deserialized(false),
79
      writeSerializationFile(false),
80
      repartitioningAllowed(false),
81
      repartitionIthChange(20),
82
83
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
84
      repartitioningFailed(0),
85
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
86
      lastMeshChangeIndex(0),
87
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
88
      boundaryDofInfo(1),
89
90
      meshAdaptivity(true),
      hasPeriodicBoundary(false)
91
  {
92
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
93

94
    mpiComm = MPI::COMM_WORLD;
95
96
    mpiRank = mpiComm.Get_rank();
    mpiSize = mpiComm.Get_size();
97

98
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
99
100
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
101
    Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
102

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

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

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

117
118
119
    if (partStr == "checker")
      partitioner = new CheckerPartitioner(&mpiComm);

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
130

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


138
  void MeshDistributor::initParallelization()
139
  {
140
    FUNCNAME("MeshDistributor::initParallelization()");
141

142
143
144
    if (initialized)
      return;

145
146
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");
147
    TEST_EXIT(feSpaces.size() > 0)
148
      ("No FE space has been defined for the mesh distributor!\n");
149
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
150

151
152
153
154
155
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

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

168
    elObjDb.setFeSpace(feSpaces[0]);
169

Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
172
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

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

181
      removePeriodicBoundaryConditions();
182

183
      elObjDb.createMacroElementInfo(allMacroElements);
184

Thomas Witkowski's avatar
Thomas Witkowski committed
185
      updateLocalGlobalNumbering();
186

Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
189
#if (DEBUG != 0)
190
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
191
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], 
192
				    *(dofMaps[0]), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
193
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
194
#endif    
195
196

      initialized = true;
197
      return;
198
    }
199

200

201
202
203
    // 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.
204
205
    testForMacroMesh();

206
207
    // For later mesh repartitioning, we need to store some information about
    // the macro mesh.
208
209
    createMacroElementInfo();

210
    // create an initial partitioning of the mesh
211
    partitioner->createInitialPartitioning();
212

213
    // set the element weights, which are 1 at the very first begin
214
    setInitialElementWeights();
215
216

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

221

222
#if (DEBUG != 0)    
223
224
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
225
226
#endif

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

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

241
    // Create interior boundary information.
242
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
243

244
245
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

246
247
248
249
250
    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))) {
251
252
253

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

254
255
256
257
258
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

    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))) {
269
270
271

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

272
273
274
275
276
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

279
280
281
	}
      }
    }
282

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
283
284
285
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
286

287
288
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
289

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

297
    updateLocalGlobalNumbering();
298

299
300
301
    // === 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.                         ===
302

303
304
    check3dValidMesh();

305

306
307
    // === If in debug mode, make some tests. ===

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

311
    ParallelDebug::testAllElements(*this);
312
    debug::testSortedDofs(mesh, elMap);
313
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
314
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
315

316
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
317
318

    MSG("Debug mode tests finished!\n");
319
#endif
320

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

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

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

336
      updateLocalGlobalNumbering();
337

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

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
343
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
344
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
345
346
      elObjDb.clear();

347
    initialized = true;
348
349
  }

350

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

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

358

359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
    // === Add all FE spaces from stationary problem. ===

    vector<const FiniteElemSpace*> newFeSpaces = probStat->getFeSpaces();
    for (int i = 0; i < static_cast<int>(newFeSpaces.size()); i++)
      if (find(feSpaces.begin(), feSpaces.end(), newFeSpaces[i]) == 
	  feSpaces.end())
	feSpaces.push_back(newFeSpaces[i]);
    

    // === Add mesh of stationary problem and create a corresponding  ===
    // === refinement manager object.                                 ===
    
    if (mesh != NULL) {
      TEST_EXIT(mesh == probStat->getMesh())
	("Does not yet support for different meshes!\n");
    } else {
      mesh = probStat->getMesh();
    }
377
378
379
380
381
382
383
384
385
386
    
    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());
387
    }
388
389
390
    
    partitioner->setMesh(mesh);
    
391

392
393
394
    // === Check whether the stationary problem should be serialized. ===


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

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

414
415
416

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

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

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

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

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

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

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

463
    problemStat.push_back(probStat);
464

465
466
467

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

469
    if (initialized)
470
      removePeriodicBoundaryConditions(probStat);
471
472
473
  }


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

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

    globalMeshDistributor->addProblemStat(probStat);
  }


485
  void MeshDistributor::exitParallelization()
486
  {}
487

488
489
490
491
492
493
494
495
496
497
498
499

  void MeshDistributor::registerDofMap(ParallelDofMapping &dofMap)
  {
    FUNCNAME("MeshDistributor::registerDofMap()");

    TEST_EXIT(find(dofMaps.begin(), dofMaps.end(), &dofMap) ==
	      dofMaps.end())
      ("Parallel DOF mapping already registerd in mesh distributor object!\n");

    dofMaps.push_back(&dofMap);
  }

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

525

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
548
549
550
551
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
552
    }
553
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
554
555
556
557
558
559
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.recv(rank);
    }
560

561
    stdMpi.startCommunication();
562

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

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

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

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

583

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

611

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

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

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

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

626
627
628
629
630
631
632

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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


735
736
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
737
738
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

739
    // Remove periodic boundaries in boundary manager on matrices and vectors.
740
741
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
742
743
744

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
745
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
746
747
748
749
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
750
751
752

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
753
754
755
  }


756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
  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()));
    }
  }


778
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
779
  {
780
781
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
782
783
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
784
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
785
	boundaryMap.erase(it++);
786
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
787
788
789
790
791
	++it;      
    }    
  }


792
793
794
795
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

796
797
798
    if (mpiSize != 16)
      return;

799
    std::set<int> neighbours;
800
801
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
802
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
803
804
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
805
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
806
807
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
808
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
809
810
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
811
812
813
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
814
815
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
816
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
817
818
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
819
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
820
821
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
822
823
824
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
825
826
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
827
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);
828
829
      break;
    case 8:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
830
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(9);  neighbours.insert(10); neighbours.insert(11);
831
832
      break;
    case 9:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
833
834
835
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(6); 
      neighbours.insert(8); neighbours.insert(12);
      neighbours.insert(10); neighbours.insert(11); neighbours.insert(14); 
836
837
      break;
    case 10:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
838
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(11);
839
840
      break;
    case 11:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
841
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12);  neighbours.insert(10); neighbours.insert(14);
842
843
      break;
    case 12:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
844
845
846
      neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); 
      neighbours.insert(9); neighbours.insert(13);
      neighbours.insert(11); neighbours.insert(14); neighbours.insert(15); 
847
848
      break;
    case 13:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
849
      neighbours.insert(6); neighbours.insert(7); neighbours.insert(12);  neighbours.insert(14); neighbours.insert(15);
850
851
      break;
    case 14:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
852
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13);  neighbours.insert(11); neighbours.insert(15);
853
854
      break;
    case 15:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
855
      neighbours.insert(12); neighbours.insert(13); neighbours.insert(14);
856
857
      break;
    }    
858

859
    TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
860
861
862
863
864
865
866

    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
867
      case 0: case 1: case 2: case 3:
868
869
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
870
	  m.insert(0); m.insert(1); m.insert(2); m.insert(3);
871
	  levelData.addLevel(m, 0);
872
873
	}
	break;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
874
      case 4: case 5: case 6: case 7:
875
876
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
877
	  m.insert(4); m.insert(5); m.insert(6); m.insert(7);
878
	  levelData.addLevel(m, 1);