MeshDistributor.cc 76.4 KB
Newer Older
1
//
2
3
4
5
6
7
8
9
10
11
12
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include <algorithm>
14
15
#include <iostream>
#include <fstream>
16
17
#include <limits>
#include <stdint.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
18
#include <boost/lexical_cast.hpp>
19
20
#include <boost/filesystem.hpp>

21
#include "parallel/MeshDistributor.h"
22
#include "parallel/MeshManipulation.h"
23
#include "parallel/ParallelDebug.h"
24
#include "parallel/StdMpi.h"
25
#include "parallel/MeshPartitioner.h"
26
#include "parallel/ParMetisPartitioner.h"
27
#include "parallel/ZoltanPartitioner.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28
#include "parallel/SimplePartitioner.h"
29
#include "parallel/CheckerPartitioner.h"
30
#include "parallel/MpiHelper.h"
31
#include "parallel/DofComm.h"
32
33
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
34
#include "io/MacroWriter.h"
35
#include "io/VtkWriter.h"
36
37
38
39
40
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
41
42
#include "DOFMatrix.h"
#include "DOFVector.h"
43
#include "SystemVector.h"
44
#include "ElementDofIterator.h"
45
46
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
47
#include "VertexVector.h"
48
#include "MeshStructure.h"
49
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
50
#include "ProblemInstat.h"
51
#include "RefinementManager3d.h"
52
#include "Debug.h"
53

54
55
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
56
  using boost::lexical_cast;
57
  using namespace boost::filesystem;
58
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
59

60
61
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

62
63
64
  const Flag MeshDistributor::BOUNDARY_SUBOBJ_SORTED              = 0X01L;
  const Flag MeshDistributor::BOUNDARY_FILL_INFO_SEND_DOFS        = 0X02L;
  const Flag MeshDistributor::BOUNDARY_FILL_INFO_RECV_DOFS        = 0X04L;
Thomas Witkowski's avatar
Thomas Witkowski committed
65

66
67
68
69
70
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

71

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
90
      createBoundaryDofFlag(0),
      sebastianMode(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);
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.setMesh(feSpaces[0]->getMesh());
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
175
      macroElIndexMap.clear();
      macroElIndexTypeMap.clear();
176

177
      for (vector<MacroElement*>::iterator it = allMacroElements.begin();
178
	   it != allMacroElements.end(); ++it) {
179
180
	macroElIndexMap.insert(make_pair((*it)->getIndex(), (*it)->getElement()));
	macroElIndexTypeMap.insert(make_pair((*it)->getIndex(), (*it)->getElType()));
181
182
      }

183
184
      createBoundaryDofs();

Thomas Witkowski's avatar
Thomas Witkowski committed
185
#if (DEBUG != 0)
186
      ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
187
#endif    
188
189

      initialized = true;
190
      return;
191
    }
192

193

194
195
196
197
198
    // 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();

199
200
201
202
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

203
    // create an initial partitioning of the mesh
204
    partitioner->createInitialPartitioning();
205

206
    // set the element weights, which are 1 at the very first begin
207
    setInitialElementWeights();
208
209

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

214

215
#if (DEBUG != 0)    
216
217
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
218
219
#endif

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

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

234

235
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
236

237
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
238

239
#if (DEBUG != 0)    
240
    ParallelDebug::printBoundaryInfo(*this);
241
#endif
242
243
244
    
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

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

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

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

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

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

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

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

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

278
279
280
	}
      }
    }
281

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

    removeMacroElements();
285

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

Thomas Witkowski's avatar
Thomas Witkowski committed
288

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

296
    updateLocalGlobalNumbering();
297

298
299
300
    // === 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.                         ===
301
302
303
    
    check3dValidMesh();

304
305
    // === If in debug mode, make some tests. ===

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

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

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

    MSG("Debug mode tests finished!\n");
317
#endif
318

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

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

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

329
330
331
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

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

340
      updateLocalGlobalNumbering();
341

342
343
      // === Update periodic mapping, if there are periodic boundaries. ===     

344
      createPeriodicMap();
345
346
347
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
348
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
349

350
    // Set DOF rank information to all matrices and vectors.
351
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
352

353
    initialized = true;
354
355
  }

356

357
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
358
  {
359
    FUNCNAME("MeshDistributor::addProblemStat()");
360

361
362
363
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

364
365
366
367
368
369
370
371
372
373
    // === 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");
374
375
      }

376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
      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());
394
    }
395
396
397
    
    partitioner->setMesh(mesh);
    
398
399
400

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

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

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

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

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

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

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

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

464
    problemStat.push_back(probStat);
465

466
467
468
    // 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.
469

470
    if (initialized) {
471
      setRankDofs(probStat);
472
473
      removePeriodicBoundaryConditions(probStat);
    }
474
475
476
  }


477
478
479
480
481
482
483
484
485
486
487
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


488
  void MeshDistributor::exitParallelization()
489
  {}
490

491
  
492
  void MeshDistributor::testForMacroMesh()
493
  {
494
    FUNCNAME("MeshDistributor::testForMacroMesh()");
495
496
497
498
499
500
501

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
506
507
508
509
510
511
512
513
514
515
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

516

517
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
518
  {
519
520
    FUNCNAME("MeshDistributor::synchVector()");

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

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

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

529
530
531
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
532

533
      stdMpi.send(it.getRank(), dofs);
534
    }
535
536
537
	   
    for (DofComm::Iterator it(recvDofs); !it.end(); it.nextRank())
      stdMpi.recv(it.getRank());
538

539
    stdMpi.startCommunication();
540

541
    for (DofComm::Iterator it(recvDofs); !it.end(); it.nextRank()) {
542
      int counter = 0;
543
544
545
546
547
548

      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++];
549
550
551
552
      }
    }
  }

553

554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
  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);
569
	allEdges.insert(elObjDb.getEdgeLocalMap(elData));
570
571
572
573
574
575
576
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

577
578

    FixRefinementPatch::connectedEdges.clear();
579
    bool valid3dMesh = true;
580
581

    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
582
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
583
584
585
586
587

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

      std::set<int> el0, el1;
588
      map<int, int> edgeNoInEl;
589
590
591

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
592
593
594
595
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
596
597
598
599

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
            
      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()) {
625
626
627
628
	valid3dMesh = false;

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

630
631
632
633
634
635
636
637
638
639
640
	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 = 
	      make_pair(macroElIndexMap[*elIt0], edgeNoInEl[*elIt0]);
	    pair<Element*, int> edge1 = 
	      make_pair(macroElIndexMap[*elIt1], edgeNoInEl[*elIt1]);

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

	    WorldVector<double> c0, c1;
641
642
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
643

644
	    /*
645
646
647
648
	    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]);
649
	    */
650
651
652
653
654
655
656
657
658

	    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));
	  }
	}
      }
    }
659
660

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
661
662
663
  }


664
665
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
					   DofContainer& dofs)
666
667
668
669
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
670
671
672
673
    for (DofComm::Iterator it(sendDofs, feSpace); !it.end(); it.nextRank())
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
    for (DofComm::Iterator it(recvDofs, feSpace); !it.end(); it.nextRank())
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
674
675
676
677
678
679

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


680
  void MeshDistributor::setRankDofs(ProblemStatSeq *probStat)
681
  {
682
683
    FUNCNAME("MeshDistributor::setRankDofs()");

684
    int nComponents = probStat->getNumComponents();
685
686
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
687
688
689
690
691
692
693
694
695
696
	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");

	  probStat->getSystemMatrix(i, j)->setRankDofs(dofFeData[rowFeSpace].isRankDof);
	}
    
697

698
      TEST_EXIT_DBG(probStat->getRhs()->getDOFVector(i))("No RHS vector!\n");
699
700
      TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))
	("No solution vector!\n");
701
702
703
      TEST_EXIT_DBG(probStat->getRhsVector(i)->getFeSpace() ==
		    probStat->getSolution(i)->getFeSpace())
	("Should not happen!\n");
704
      
705
706
707
708
709
710
      const FiniteElemSpace *feSpace = probStat->getRhsVector(i)->getFeSpace();
      TEST_EXIT(find(feSpaces.begin(), feSpaces.end(), feSpace) != feSpaces.end())
	("Should not happen!\n");

      probStat->getRhsVector(i)->setRankDofs(dofFeData[feSpace].isRankDof);
      probStat->getSolution(i)->setRankDofs(dofFeData[feSpace].isRankDof);
711
712
713
714
    }    
  }


715
716
  void MeshDistributor::setRankDofs()
  {
717
718
    for (unsigned int i = 0; i < problemStat.size(); i++) 
      setRankDofs(problemStat[i]);
719
720
721
  }


722
723
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
724
725
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

726
    // Remove periodic boundaries in boundary manager on matrices and vectors.
727
728
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
729
730
731
732
733
734
735
736

    // 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);
    }    
737
738
739

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
740
741
742
  }


743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
  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()));
    }
  }


765
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
766
  {
767
768
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
769
770
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
771
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
772
	boundaryMap.erase(it++);
773
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
774
775
776
777
778
	++it;      
    }    
  }


779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
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
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

    if (mpiSize != 16)
      return;

    std::set<int> neighbours;
    for (InteriorBoundary::iterator it(rankIntBoundary); !it.end(); ++it)
      neighbours.insert(it.getRank());

    for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
      neighbours.insert(it.getRank());

    levelData.init(neighbours);

    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);
	  levelData.addLevel(m);
	}
	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);
	  levelData.addLevel(m);
	}
	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);
	  levelData.addLevel(m);
	}
	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);
	  levelData.addLevel(m);
	}
	break;
      }
    }
  }


832
  void MeshDistributor::checkMeshChange(bool tryRepartition)
833
  {
834
    FUNCNAME("MeshDistributor::checkMeshChange()");
835

836
837
    double first = MPI::Wtime();

838
839
840
    int skip = 0;
    Parameters::get("parallel->debug->skip check mesh change", skip);

841
842
843
844
845
    // === If mesh has not been changed on all ranks, return. ===

    int recvAllValues = 0;
    int sendValue = static_cast<int>(mesh->getChangeIndex() != lastMeshChangeIndex);
    mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
846

847
    if (recvAllValues == 0)
848
849
      return;

850
851
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
852
853

    if (skip == 0)
854
    do {
855
856
      bool meshChanged = false;

857
858
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
859
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
860

861
      for (InteriorBoundary::iterator it(rankIntBoundary); !it.end(); ++it)
862
863
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
864
 	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
865
866

      for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
867
868
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
869
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
870

871
872
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it) {
	if (it.getRank() == mpiRank) {
873
874
875
876
877
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE)) {
	    MeshStructure elCode;
	    elCode.init(it->rankObj);
	    
878
	    MeshManipulation mm(mesh);
879
	    meshChanged |= mm.fitElementToMeshCode(elCode, it->neighObj);
880
	  }
881
882
883
884
885
886
	} else {
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
	    allBound[it.getRank()].push_back(*it);	
	}
      }
887

Thomas Witkowski's avatar