MeshDistributor.cc 61.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
      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
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
232

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

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

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

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

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

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

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

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

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

268
269
270
	}
      }
    }
271

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

    removeMacroElements();
275

276
277
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
278

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

286
    updateLocalGlobalNumbering();
287

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

294
295
    // === If in debug mode, make some tests. ===

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

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

304
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
305
306

    MSG("Debug mode tests finished!\n");
307
#endif
308

309
    // Create periodic DOF mapping, if there are periodic boundaries.
310
    createPeriodicMap();
311

312
313
314
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

315
316
317
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
318

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

327
      updateLocalGlobalNumbering();
328

329
330
      // === Update periodic mapping, if there are periodic boundaries. ===     

331
      createPeriodicMap();
332
333
334
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
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
    initialized = true;
341
342
  }

343

344
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
345
  {
346
    FUNCNAME("MeshDistributor::addProblemStat()");
347

348
349
350
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

351
352
353
354
355
356
357
358
359
360
    // === 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");
361
362
      }

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

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

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

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

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

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

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

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

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

451
    problemStat.push_back(probStat);
452

453
454
455
    // 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.
456

457
    if (initialized) {
458
      setRankDofs(probStat);
459
460
      removePeriodicBoundaryConditions(probStat);
    }
461
462
463
  }


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

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

    globalMeshDistributor->addProblemStat(probStat);
  }


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

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

503

504
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
505
  {
506
507
    FUNCNAME("MeshDistributor::synchVector()");

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

510
511
    for (DofComm::Iterator it(dofComm.getSendDofs()); 
	 !it.end(); it.nextRank()) {
512
      vector<double> dofs;
513

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

517
518
519
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
520

521
      stdMpi.send(it.getRank(), dofs);
522
    }
523
	   
524
    for (DofComm::Iterator it(dofComm.getRecvDofs()); !it.end(); it.nextRank())
525
      stdMpi.recv(it.getRank());
526

527
    stdMpi.startCommunication();
528

529
    for (DofComm::Iterator it(dofComm.getRecvDofs()); !it.end(); it.nextRank()) {
530
      int counter = 0;
531
532
533
534
535
536

      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++];
537
538
539
540
      }
    }
  }

541

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

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

      elInfo = stack.traverseNext(elInfo);
    }

565
566

    FixRefinementPatch::connectedEdges.clear();
567
    bool valid3dMesh = true;
568
569

    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
570
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
571
572
573
574
575

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

      std::set<int> el0, el1;
576
      map<int, int> edgeNoInEl;
577
578
579

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

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

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

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

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

	    WorldVector<double> c0, c1;
629
630
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
631

632
	    /*
633
634
635
636
	    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]);
637
	    */
638
639
640
641
642
643
644
645
646

	    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));
	  }
	}
      }
    }
647
648

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
649
650
651
  }


652
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
653
					   int level,
654
					   DofContainer& dofs)
655
656
657
658
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
659
660
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
661
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
662
663
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
664
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
665
666
667
668
669
670

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


671
  void MeshDistributor::setRankDofs(ProblemStatSeq *probStat)
672
  {
673
674
    FUNCNAME("MeshDistributor::setRankDofs()");

675
    int nComponents = probStat->getNumComponents();
676
677
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
678
679
680
681
682
683
684
	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
685
	  probStat->getSystemMatrix(i, j)->setDofMap(dofMap[rowFeSpace]);
686
687
	}
    
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
      TEST_EXIT_DBG(probStat->getRhsVector(i)->getFeSpace() ==
		    probStat->getSolution(i)->getFeSpace())
	("Should not happen!\n");
695
      
696
697
698
699
      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
700
701
      probStat->getRhsVector(i)->setDofMap(dofMap[feSpace]);
      probStat->getSolution(i)->setDofMap(dofMap[feSpace]);
702
703
704
705
    }    
  }


706
707
  void MeshDistributor::setRankDofs()
  {
708
709
    for (unsigned int i = 0; i < problemStat.size(); i++) 
      setRankDofs(problemStat[i]);
710
711
712
  }


713
714
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
715
716
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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

    // 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);
    }    
728
729
730

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
731
732
733
  }


734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
  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()));
    }
  }


756
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
757
  {
758
759
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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


770
771
772
773
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

774
775
776
    if (mpiSize != 16)
      return;

777
    std::set<int> neighbours;
778
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
832
833
834
835
    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;
    }    
836

837
    TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
838
839
840

    levelData.init(neighbours);

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

843
844
845
846
847
848
849
850
    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);
851
	  levelData.addLevel(m, 0);
852
853
854
855
856
857
	}
	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);
858
	  levelData.addLevel(m, 1);
859
860
861
862
863
864
	}
	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);
865
	  levelData.addLevel(m, 2);
866
867
868
869
870
871
	}
	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);
872
	  levelData.addLevel(m, 3);
873
874
875
876
877
878
879
	}
	break;
      }
    }
  }


880
  void MeshDistributor::checkMeshChange(bool tryRepartition)
881
  {
882
    FUNCNAME("MeshDistributor::checkMeshChange()");
883

884
885
    double first = MPI::Wtime();

886
887