MeshDistributor.cc 75.2 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
91
      sebastianMode(false),
      boundaryDofInfo(1)
92
  {
93
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
94

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

99
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
100
101
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
102
    Parameters::get(name + "->log main rank", Msg::outputMainRank);
103

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
131

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

136
137
138
    if (initialized)
      return;

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

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

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

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

163
    elObjDb.setMesh(feSpaces[0]->getMesh());
164

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

171
      setRankDofs();
172

173
      removePeriodicBoundaryConditions();
174

175
176
      macroElIndexMap.clear();
      macroElIndexTypeMap.clear();
177

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

184
185
      createBoundaryDofs();

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

      initialized = true;
191
      return;
192
    }
193

194

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

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

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

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

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

215

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

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

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

235
236
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
237

238
    // Create interior boundary information.
239
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
240

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

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

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

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

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

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

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

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

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

280
281
282
	}
      }
    }
283

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

    removeMacroElements();
287

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

Thomas Witkowski's avatar
Thomas Witkowski committed
290

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

298
    updateLocalGlobalNumbering();
299

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

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

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

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

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

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

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

324
325
326
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

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

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

339
      updateLocalGlobalNumbering();
340

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

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

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

352
    initialized = true;
353
354
  }

355

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

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

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

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

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

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

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

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

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

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

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

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

463
    problemStat.push_back(probStat);
464

465
466
467
    // If the mesh distributor is already initialized, don't forget to set rank
    // DOFs object to the matrices and vectors of the added stationary problem,
    // and to remove the periodic boundary conditions on these objects.
468

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


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

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

    globalMeshDistributor->addProblemStat(probStat);
  }


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

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

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

515

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

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

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

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

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

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

538
    stdMpi.startCommunication();
539

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

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

552

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

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

      elInfo = stack.traverseNext(elInfo);
    }

576
577

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

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

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

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

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

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

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

629
630
631
632
633
634
635
636
637
638
639
	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;
640
641
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
642

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

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

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


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

    DofContainerSet dofSet;
670
    for (DofComm::Iterator it(sendDofs, level, feSpace); !it.end(); it.nextRank())
671
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
672
    for (DofComm::Iterator it(recvDofs, level, feSpace); !it.end(); it.nextRank())
673
      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
	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
694
	  probStat->getSystemMatrix(i, j)->setDofMap(dofMap[rowFeSpace]);
695
696
	}
    
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
      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
709
710
      probStat->getRhsVector(i)->setDofMap(dofMap[feSpace]);
      probStat->getSolution(i)->setDofMap(dofMap[feSpace]);
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
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

783
784
785
    if (mpiSize != 16)
      return;

786
    std::set<int> neighbours;
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
    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;
    }    
845

846
    TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
847
848
849

    levelData.init(neighbours);

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

852
853
854
855
856
857
858
859
    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);
860
	  levelData.addLevel(m, 0);
861
862
863
864
865
866
	}
	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);
867
	  levelData.addLevel(m, 1);
868
869
870
871
872
873
	}
	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);
874
	  levelData.addLevel(m, 2);
875
876
877
878
879
880
	}
	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);
881
	  levelData.addLevel(m, 3);
882
883
884
885
886
887
888
	}
	break;
      }
    }
  }


889
  void MeshDistributor::checkMeshChange(bool tryRepartition)
890
  {
891
    FUNCNAME("MeshDistributor::checkMeshChange()");