MeshDistributor.cc 71.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
32
33
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
#include "io/VtkWriter.h"
34
35
36
37
38
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
39
40
#include "DOFMatrix.h"
#include "DOFVector.h"
41
#include "SystemVector.h"
42
#include "ElementDofIterator.h"
43
44
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
45
#include "VertexVector.h"
46
#include "MeshStructure.h"
47
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
48
#include "ProblemInstat.h"
49
#include "RefinementManager3d.h"
50
#include "Debug.h"
51

52
53
namespace AMDiS {

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

58
59
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

60
61
62
  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
63

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

69

70
  MeshDistributor::MeshDistributor()
71
    : problemStat(0),
72
      initialized(false),
73
      name("parallel"),
74
75
76
77
78
      feSpace(NULL),
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
79
      nRankDofs(0),
80
      rStartDofs(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
81
      nOverallDofs(0),
82
      deserialized(false),
83
      writeSerializationFile(false),
84
      repartitioningAllowed(false),
85
      repartitionIthChange(20),
86
87
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
88
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
89
90
      lastMeshChangeIndex(0),
      createBoundaryDofFlag(0)
91
  {
92
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
93

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

    int tmp = 0;
99
    Parameters::get(name + "->repartitioning", tmp);
100
    repartitioningAllowed = (tmp > 0);
101

102
103
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
104
105

    tmp = 0;
106
    Parameters::get(name + "->log main rank", tmp);
107
    Msg::outputMainRank = (tmp > 0);
108

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

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

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

123
124
125
    if (partStr == "checker")
      partitioner = new CheckerPartitioner(&mpiComm);

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
136

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

141
142
143
    if (initialized)
      return;

144
145
146
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");

147
148
    TEST_EXIT(feSpace)("No FE space has been defined for the mesh distributor!\n");
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
149

150
151
    elObjects.setFeSpace(feSpace);

152
153
154
    // 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).
155
    if (deserialized) {
156
157
      updateMacroElementInfo();

158
      setRankDofs();
159

160
      removePeriodicBoundaryConditions();
161

162
163
      macroElIndexMap.clear();
      macroElIndexTypeMap.clear();
164

165
      for (vector<MacroElement*>::iterator it = allMacroElements.begin();
166
	   it != allMacroElements.end(); ++it) {
167
168
169
170
	macroElIndexMap[(*it)->getIndex()] = (*it)->getElement();
	macroElIndexTypeMap[(*it)->getIndex()] = (*it)->getElType();
      }

171
172
      createBoundaryDofs();

Thomas Witkowski's avatar
Thomas Witkowski committed
173
#if (DEBUG != 0)
174
      ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
175
#endif    
176
177

      initialized = true;
178
      return;
179
    }
180

181

182
183
184
185
186
    // 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();

187
188
189
190
    // For later mesh repartitioning, we need to store some information about the
    // macro mesh.
    createMacroElementInfo();

191
    // create an initial partitioning of the mesh
192
    partitioner->createInitialPartitioning();
193

194
    // set the element weights, which are 1 at the very first begin
195
    setInitialElementWeights();
196
197

    // and now partition the mesh    
198
199
    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
200
    partitioner->createPartitionMap(partitionMap);
201

202

203
#if (DEBUG != 0)    
204
205
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
206
207
#endif

208
    if (mpiRank == 0) {
209
#if (DEBUG != 0)    
210
      int writePartMesh = 1;
211
212
213
#else
      int writePartMesh = 0;
#endif
214
      Parameters::get("dbg->write part mesh", writePartMesh);
215

216
      if (writePartMesh > 0) {
217
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex.vtu");
218
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
219
      }
220
    }
221

222

223
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
224

225
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
226

227
#if (DEBUG != 0)
228
    ParallelDebug::printBoundaryInfo(*this);
229
230
#endif

231
232
233
    
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

234
235
236
237
238
    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))) {
239
240
241

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

242
243
244
245
246
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
247
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
248
249
250
251
252
253
254
255
256
	}
      }
    }

    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))) {
257
258
259

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

260
261
262
263
264
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
265
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
266
267
268
	}
      }
    }
269

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
270
271
272
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
273

274
275
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
276

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

284
    updateLocalGlobalNumbering();
285

286
287
288
289
290
291
    // === 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.                                           ===
    
    check3dValidMesh();

292
293
    // === If in debug mode, make some tests. ===

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

297
    ParallelDebug::testAllElements(*this);
298
    debug::testSortedDofs(mesh, elMap);
299
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
300
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
301

302
    debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh");   
303
304

    MSG("Debug mode tests finished!\n");
305
#endif
306

307
    // === Create periodic DOF mapping, if there are periodic boundaries. ===
308

309
    createPeriodicMap();
310

311
312
313
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
314

315
    // === Global refinements. ===
316
    
Thomas Witkowski's avatar
Thomas Witkowski committed
317
    int globalRefinement = 0;
318
    Parameters::get(mesh->getName() + "->global refinements", globalRefinement);
319
    
Thomas Witkowski's avatar
Thomas Witkowski committed
320
    if (globalRefinement > 0) {
321
      refineManager->globalRefine(mesh, globalRefinement);
322

323
      updateLocalGlobalNumbering();
324
325

     
326
327
      // === Update periodic mapping, if there are periodic boundaries. ===     

328
      createPeriodicMap();
329
330
331
332

#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
333
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
334

335

Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
    /// === Set DOF rank information to all matrices and vectors. ===

338
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
339

340

Thomas Witkowski's avatar
Thomas Witkowski committed
341
342
    // === Remove periodic boundary conditions in sequential problem definition. ===

343
    removePeriodicBoundaryConditions();
344
345

    initialized = true;
346
347
  }

348

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

    if (feSpace != NULL) {
354
      vector<FiniteElemSpace*> feSpaces = probStat->getFeSpaces();
355
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
356
357
358
359
360
// 	MSG("MESH %p <-> %p   BF %p <-> %p\n",
// 	    feSpace->getMesh(),
// 	    feSpaces[i]->getMesh(),
// 	    feSpace->getBasisFcts(),
// 	    feSpaces[i]->getBasisFcts());
361
362
363
364
	TEST_EXIT(feSpace == feSpaces[i])
	  ("Parallelizaton is not supported for multiple FE spaces!\n");
      }
    } else {
365
      feSpace = probStat->getFeSpace(0);
366
      mesh = feSpace->getMesh();
367
      info = probStat->getInfo();
368
369
370
      
      TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
	("Only meshes with one DOFAdmin are supported!\n");
371
      TEST_EXIT(mesh->getDofAdmin(0).getNumberOfPreDofs(VERTEX) == 0)
372
373
374
375
376
377
378
379
380
381
382
383
384
	("Wrong pre dof number for DOFAdmin!\n");
      
      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
      partitioner->setMesh(mesh);
386
387
388
389
    }

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

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

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

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

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

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

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

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

453
    problemStat.push_back(probStat);
454
455
456
457
458

    // If the mesh distributor is already initialized, don't forgett to set rank
    // DOFs object to the matrices and vectors of the added stationary problem.

    if (initialized)
459
      setRankDofs(probStat);
460
461
462
  }


463
464
465
466
467
468
469
470
471
472
473
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


474
  void MeshDistributor::exitParallelization()
475
  {}
476

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

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
492
493
494
495
496
497
498
499
500
501
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

502

503
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
504
  {
505
    int nComponents = vec.getSize();
506
    StdMpi<vector<double> > stdMpi(mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
507
508
509

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt) {
510
      vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
511
512
513
514
515
516
517
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nComponents * nSendDofs);
      
      for (int i = 0; i < nComponents; i++) {
	DOFVector<double> *dofvec = vec.getDOFVector(i);
	for (int j = 0; j < nSendDofs; j++)
	  dofs.push_back((*dofvec)[*((sendIt->second)[j])]);
518
519
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
520
      stdMpi.send(sendIt->first, dofs);
521
522
523
    }

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
524
525
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size() * nComponents);
526

527
    stdMpi.startCommunication();
528
529

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
530
531
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
532
533

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
534
535
536
537
538
      for (int i = 0; i < nComponents; i++) {
	DOFVector<double> *dofvec = vec.getDOFVector(i);
 	for (int j = 0; j < nRecvDofs; j++)
	  (*dofvec)[*(recvIt->second)[j]] = 
	    stdMpi.getRecvData(recvIt->first)[counter++];
539
540
541
542
      }
    }
  }

543

544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
  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);
	allEdges.insert(elObjects.getEdgeLocalMap(elData));
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

567
568

    FixRefinementPatch::connectedEdges.clear();
569
    bool valid3dMesh = true;
570
571
572
573
574
575
576
577

    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
      vector<ElementObjectData>& edgeEls = elObjects.getElements(*it);

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

      std::set<int> el0, el1;
578
      map<int, int> edgeNoInEl;
579
580
581

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

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
            
      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()) {
615
616
617
618
	valid3dMesh = false;

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

620
621
622
623
624
625
626
	for (std::set<int>::iterator elIt0 = el0.begin(); elIt0 != el0.end(); ++elIt0) {
	  for (std::set<int>::iterator elIt1 = el1.begin(); elIt1 != el1.end(); ++elIt1) {

	    TEST_EXIT_DBG(macroElIndexMap.count(*elIt0))("Should not happen!\n");
	    TEST_EXIT_DBG(macroElIndexMap.count(*elIt1))("Should not happen!\n");
	    TEST_EXIT_DBG(edgeNoInEl.count(*elIt0))("Should not happen!\n");
	    TEST_EXIT_DBG(edgeNoInEl.count(*elIt1))("Should not happen!\n");
627

628
629
630
631
632
633
634
635
636
637
638
639
	    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;
	    mesh->getDofIndexCoords(dofEdge0.first, feSpace, c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpace, c1);

640
641
642
643
	    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]);
644
645
646
647
648
649
650
651
652

	    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));
	  }
	}
      }
    }
653
654

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


658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
  void MeshDistributor::getAllBoundaryDofs(DofContainer& dofs)
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt)
      for (DofContainer::iterator dofIt = sendIt->second.begin();
	   dofIt != sendIt->second.end(); ++dofIt)
	dofSet.insert(*dofIt);

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      for (DofContainer::iterator dofIt = recvIt->second.begin();
	   dofIt != recvIt->second.end(); ++dofIt)
	dofSet.insert(*dofIt);
    
    dofs.clear();
    dofs.insert(dofs.begin(), dofSet.begin(), dofSet.end());
  }


681
  void MeshDistributor::setRankDofs(ProblemStatSeq *probStat)
682
  {
683
    int nComponents = probStat->getNumComponents();
684
685
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
686
687
	if (probStat->getSystemMatrix(i, j))
	  probStat->getSystemMatrix(i, j)->setRankDofs(isRankDof);
688

689
      TEST_EXIT_DBG(probStat->getRhs()->getDOFVector(i))("No RHS vector!\n");
690
691
      TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))
	("No solution vector!\n");
692
      
693
694
      probStat->getRhs()->getDOFVector(i)->setRankDofs(isRankDof);
      probStat->getSolution()->getDOFVector(i)->setRankDofs(isRankDof);
695
696
697
698
    }    
  }


699
700
  void MeshDistributor::setRankDofs()
  {
701
702
    for (unsigned int i = 0; i < problemStat.size(); i++) 
      setRankDofs(problemStat[i]);
703
704
705
  }


706
707
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
708
709
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

710
    // Remove periodic boundaries in boundary manager on matrices and vectors.
711
712
    for (unsigned int i = 0; i < problemStat.size(); i++) {
      int nComponents = problemStat[i]->getNumComponents();
713
714
715

      for (int j = 0; j < nComponents; j++) {
	for (int k = 0; k < nComponents; k++) {
716
	  DOFMatrix* mat = problemStat[i]->getSystemMatrix(j, k);
717
	  if (mat && mat->getBoundaryManager())
718
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
719
720
	}
	
721
722
	if (problemStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(problemStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
723
	
724
725
	if (problemStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(problemStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
726
727
728
729
730
731
732
733
734
735
      }
    }

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

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


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
745
746
747
748
749
750
751
752
753
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


754
  void MeshDistributor::checkMeshChange(bool tryRepartition)
755
  {
756
    FUNCNAME("MeshDistributor::checkMeshChange()");
757

758
759
    double first = MPI::Wtime();

760
761
762
763
764
    // === 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);
765

766
    if (recvAllValues == 0)
767
768
      return;

769
770
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
771
   
772
    do {
773
774
      bool meshChanged = false;

775
776
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
777
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
778
779

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
780
781
	if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	    (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
782
 	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
783
784

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

789
790
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it) {
	if (it.getRank() == mpiRank) {
791
792
793
794
795
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE)) {
	    MeshStructure elCode;
	    elCode.init(it->rankObj);
	    
796
	    MeshManipulation mm(feSpace);
797
	    meshChanged |= mm.fitElementToMeshCode(elCode, it->neighObj);
798
	  }
799
800
801
802
803
804
	} else {
	  if ((mesh->getDim() == 2 && it->rankObj.subObj == EDGE) || 
	      (mesh->getDim() == 3 && it->rankObj.subObj == FACE))
	    allBound[it.getRank()].push_back(*it);	
	}
      }
805

806

807
      // === Check the boundaries and adapt mesh if necessary. ===
808
      MSG_DBG("Run checkAndAdaptBoundary ...\n");
809

810
      meshChanged |= checkAndAdaptBoundary(allBound);
811
812
813

      // === Check on all ranks if at least one rank's mesh has changed. ===

814
      int sendValue = static_cast<int>(meshChanged);
815
816
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
817
818

      MSG("Mesh changed on %d ranks!\n", recvAllValues);
819
    } while (recvAllValues != 0);
820

821
#if (DEBUG != 0)
822
    debug::writeMesh(feSpace, -1, debugOutputDir + "mesh");
823
824
#endif

825

826
    // === Because the mesh has been changed, update the DOF numbering and mappings. ===
827

828
829
    updateLocalGlobalNumbering();

830

831
    // === Update periodic mapping, if there are periodic boundaries. ===
832

Thomas Witkowski's avatar
Thomas Witkowski committed
833
    createPeriodicMap();
834

835
836
837
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
838
839


840
841
    // === The mesh has changed, so check if it is required to repartition the mesh. ===

842
    nMeshChangesAfterLastRepartitioning++;
843

844
845
846
847

    INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n", 
		  MPI::Wtime() - first);

848
849
    if (tryRepartition &&
	repartitioningAllowed && 
850
851
852
	nMeshChangesAfterLastRepartitioning >= repartitionIthChange) {
      repartitionMesh();
      nMeshChangesAfterLastRepartitioning = 0;
853
854
855
856
    } else {
      MSG_DBG("Repartitioning not tried because tryRepartitioning = %d   reparttitioningAllowed = %d   nMeshChange =%d    repartitionIthChange = %d\n",
	      tryRepartition, repartitioningAllowed, 
	      nMeshChangesAfterLastRepartitioning, repartitionIthChange);
857
    }
858
859
860