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

72
  MeshDistributor::MeshDistributor()
73
    : problemStat(0),
74
      initialized(false),
75
      name("parallel"),
76
      feSpaces(0),
77
78
79
80
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
81
      deserialized(false),
82
      writeSerializationFile(false),
83
      repartitioningAllowed(false),
84
      repartitionIthChange(20),
85
86
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
87
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
88
      lastMeshChangeIndex(0),
89
      createBoundaryDofFlag(0),
90
      boundaryDofInfo(1)
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);
101
    Parameters::get(name + "->log main rank", Msg::outputMainRank);
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

131
  void MeshDistributor::initParallelization()
132
  {
133
    FUNCNAME("MeshDistributor::initParallelization()");
134

135
136
137
    if (initialized)
      return;

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

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

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

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

162
    elObjDb.setFeSpace(feSpaces[0]);
163

164
165
166
    // 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).
167
    if (deserialized) {
168
169
      updateMacroElementInfo();

170
      setRankDofs();
171

172
      removePeriodicBoundaryConditions();
173

174
      elObjDb.createMacroElementInfo(allMacroElements);
175

176
177
      createBoundaryDofs();

Thomas Witkowski's avatar
Thomas Witkowski committed
178
#if (DEBUG != 0)
179
      ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
180
#endif    
181
182

      initialized = true;
183
      return;
184
    }
185

186

187
188
189
190
191
    // 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.
    testForMacroMesh();

192
193
194
195
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

196
    // create an initial partitioning of the mesh
197
    partitioner->createInitialPartitioning();
198

199
    // set the element weights, which are 1 at the very first begin
200
    setInitialElementWeights();
201
202

    // and now partition the mesh    
203
204
    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
205
    partitioner->createPartitionMap(partitionMap);
206

207

208
#if (DEBUG != 0)    
209
210
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
211
212
#endif

213
    if (mpiRank == 0) {
214
#if (DEBUG != 0)    
215
      int writePartMesh = 1;
216
217
218
#else
      int writePartMesh = 0;
#endif
219
      Parameters::get("dbg->write part mesh", writePartMesh);
220

221
      if (writePartMesh > 0) {
222
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex.vtu");
223
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
224
      }
225
    }
226

227
228
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
229

230
    // Create interior boundary information.
231
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
232

233
#if (DEBUG != 0)    
234
    ParallelDebug::printBoundaryInfo(*this);
235
#endif
236
237
238
    
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

239
240
241
242
243
    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))) {
244
245
246

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

247
248
249
250
251
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

    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))) {
262
263
264

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

265
266
267
268
269
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

272
273
274
	}
      }
    }
275

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
276
277
278
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
279

280
281
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
282

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

290
    updateLocalGlobalNumbering();
291

292
293
294
    // === 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.                         ===
295
296
297
    
    check3dValidMesh();

298
299
    // === If in debug mode, make some tests. ===

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

303
    ParallelDebug::testAllElements(*this);
304
    debug::testSortedDofs(mesh, elMap);
305
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
306
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
307

308
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
309
310

    MSG("Debug mode tests finished!\n");
311
#endif
312

313
    // Create periodic DOF mapping, if there are periodic boundaries.
314
    createPeriodicMap();
315

316
317
318
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

319
320
321
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
322

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

331
      updateLocalGlobalNumbering();
332

333
334
      // === Update periodic mapping, if there are periodic boundaries. ===     

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

341
    // Set DOF rank information to all matrices and vectors.
342
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
343

344
    initialized = true;
345
346
  }

347

348
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
349
  {
350
    FUNCNAME("MeshDistributor::addProblemStat()");
351

352
353
354
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

355
356
357
358
359
360
361
362
363
364
    // === 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");
365
366
      }

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

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

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

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

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

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

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

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

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

455
    problemStat.push_back(probStat);
456

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

461
    if (initialized) {
462
      setRankDofs(probStat);
463
464
      removePeriodicBoundaryConditions(probStat);
    }
465
466
467
  }


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

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

    globalMeshDistributor->addProblemStat(probStat);
  }


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

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

507

508
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
509
  {
510
511
    FUNCNAME("MeshDistributor::synchVector()");

512
    StdMpi<vector<double> > stdMpi(mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
513

514
515
    for (DofComm::Iterator it(dofComm.getSendDofs()); 
	 !it.end(); it.nextRank()) {
516
      vector<double> dofs;
517

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

521
522
523
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
524

525
      stdMpi.send(it.getRank(), dofs);
526
    }
527
	   
528
    for (DofComm::Iterator it(dofComm.getRecvDofs()); !it.end(); it.nextRank())
529
      stdMpi.recv(it.getRank());
530

531
    stdMpi.startCommunication();
532

533
    for (DofComm::Iterator it(dofComm.getRecvDofs()); !it.end(); it.nextRank()) {
534
      int counter = 0;
535
536
537
538
539
540

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

	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
	  dofVec[it.getDofIndex()] = stdMpi.getRecvData(it.getRank())[counter++];
541
542
543
544
      }
    }
  }

545

546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
      for (int i = 0; i < 4; i++) {
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
561
	allEdges.insert(elObjDb.getEdgeLocalMap(elData));
562
563
564
565
566
567
568
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

569
570

    FixRefinementPatch::connectedEdges.clear();
571
    bool valid3dMesh = true;
572
573

    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
574
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
575
576
577
578
579

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

      std::set<int> el0, el1;
580
      map<int, int> edgeNoInEl;
581
582
583

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
584
585
586
587
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
588
589
590
591

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
            
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

      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()) {
617
618
619
620
	valid3dMesh = false;

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

622
623
624
	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 = 
625
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
626
	    pair<Element*, int> edge1 = 
627
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
628
629
630
631
632

	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
633
634
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
635

636
	    /*
637
638
639
640
	    MSG_DBG("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]);
641
	    */
642
643
644
645
646
647
648
649
650

	    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));
	  }
	}
      }
    }
651
652

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
653
654
655
  }


656
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
657
					   int level,
658
					   DofContainer& dofs)
659
660
661
662
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
663
664
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
665
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
666
667
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
668
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
669
670
671
672
673
674

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


675
  void MeshDistributor::setRankDofs(ProblemStatSeq *probStat)
676
  {
677
678
    FUNCNAME("MeshDistributor::setRankDofs()");

679
    int nComponents = probStat->getNumComponents();
680
681
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
682
683
684
685
686
687
688
	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
689
	  probStat->getSystemMatrix(i, j)->setDofMap(dofMap[rowFeSpace]);
690
691
	}
    
692

693
      TEST_EXIT_DBG(probStat->getRhs()->getDOFVector(i))("No RHS vector!\n");
694
695
      TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))
	("No solution vector!\n");
696
697
698
      TEST_EXIT_DBG(probStat->getRhsVector(i)->getFeSpace() ==
		    probStat->getSolution(i)->getFeSpace())
	("Should not happen!\n");
699
      
700
701
702
703
      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
704
705
      probStat->getRhsVector(i)->setDofMap(dofMap[feSpace]);
      probStat->getSolution(i)->setDofMap(dofMap[feSpace]);
706
707
708
709
    }    
  }


710
711
  void MeshDistributor::setRankDofs()
  {
712
713
    for (unsigned int i = 0; i < problemStat.size(); i++) 
      setRankDofs(problemStat[i]);
714
715
716
  }


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

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

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

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


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


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

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


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

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

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

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

    levelData.init(neighbours);

845
846
    MSG("INIT MESH LEVEL %d\n", levelData.getLevelNumber());

847
848
849
850
851
852
853
854
    bool multiLevelTest = false;
    Parameters::get("parallel->multi level test", multiLevelTest);
    if (multiLevelTest) {
      switch (mpiRank) {
      case 0: case 1: case 4: case 5:
	{
	  std::set<int> m;
	  m.insert(0); m.insert(1); m.insert(4); m.insert(5);
855
	  levelData.addLevel(m, 0);
856
857
858
859
860
861
	}
	break;
      case 2: case 3: case 6: case 7:
	{
	  std::set<int> m;
	  m.insert(2); m.insert(3); m.insert(6); m.insert(7);
862
	  levelData.addLevel(m, 1);
863
864
865
866
867
868
	}
	break;
      case 8: case 9: case 12: case 13:
	{
	  std::set<int> m;
	  m.insert(8); m.insert(9); m.insert(12); m.insert(13);
869
	  levelData.addLevel(m, 2);
870
871
872
873
874
875
	}
	break;
      case 10: case 11: case 14: case 15:
	{
	  std::set<int> m;
	  m.insert(10); m.insert(11); m.insert(14); m.insert(15);
876
	  levelData.addLevel(m, 3);
877
878
879
880
881
882
883
	}
	break;
      }
    }
  }


884
  void MeshDistributor::checkMeshChange(bool tryRepartition)
885
  {
886
    FUNCNAME("MeshDistributor::checkMeshChange()");