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

56
57
namespace AMDiS {

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

62
63
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

64
65
66
  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
67

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

73
  MeshDistributor::MeshDistributor()
74
    : problemStat(0),
75
      initialized(false),
76
      name("parallel"),
77
78
79
      mesh(NULL),
      refineManager(NULL),
      partitioner(NULL),
80
      deserialized(false),
81
      writeSerializationFile(false),
82
      repartitioningAllowed(false),
83
      repartitionIthChange(20),
84
85
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
86
      repartitioningFailed(0),
87
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
88
      lastMeshChangeIndex(0),
89
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
90
      boundaryDofInfo(1),
91
      meshAdaptivity(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
92
93
      hasPeriodicBoundary(false),
      printTimings(false)
94
  {
95
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
96

97
    mpiComm = MPI::COMM_WORLD;
98
99
    mpiRank = mpiComm.Get_rank();
    mpiSize = mpiComm.Get_size();
100

101
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
102
103
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
104
    Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
Thomas Witkowski's avatar
Thomas Witkowski committed
105
    
106
    string partStr = "parmetis";
107
    Parameters::get(name + "->partitioner", partStr);
108
109
110
111

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

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

120
121
122
    if (partStr == "checker")
      partitioner = new CheckerPartitioner(&mpiComm);

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
130
131
    Parameters::get(name + "->print timings", printTimings);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
135

Thomas Witkowski's avatar
Thomas Witkowski committed
136
137
138
139
140
141
142
  MeshDistributor::~MeshDistributor()
  {
    if (partitioner)
      delete partitioner;
  }


143
  void MeshDistributor::initParallelization()
144
  {
145
    FUNCNAME("MeshDistributor::initParallelization()");
146

147
148
149
    if (initialized)
      return;

Thomas Witkowski's avatar
Thomas Witkowski committed
150
151
152
153
    double first = MPI::Wtime();
    MSG("Initialization phase 1 needed %.5f seconds\n", 
	first - ParallelProblemStatBase::initTimeStamp);

154
155
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");
156
    TEST_EXIT(feSpaces.size() > 0)
157
      ("No FE space has been defined for the mesh distributor!\n");
158
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
159

160
161
162
163
164
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
165
      doNext = false;
166
167
168
169
170
171
      for (unsigned int i = 0; i < feSpaces.size() - 1; i++) {
	if (feSpaces[i]->getBasisFcts()->getDegree() >
	    feSpaces[i + 1]->getBasisFcts()->getDegree()) {
	  const FiniteElemSpace *tmp = feSpaces[i + 1];
	  feSpaces[i + 1] = feSpaces[i];
	  feSpaces[i] = tmp;
172
173
	  doNext = true;
	}
174
      }
175
176
    } while (doNext);

177
    elObjDb.setFeSpace(feSpaces[0]);
178

Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
181
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

182
183
184
    // 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).
185
    if (deserialized) {
Thomas Witkowski's avatar
Thomas Witkowski committed
186
187
      createMeshLevelStructure();
      
188
189
      updateMacroElementInfo();

190
      removePeriodicBoundaryConditions();
191

192
      elObjDb.createMacroElementInfo(allMacroElements);
193

Thomas Witkowski's avatar
Thomas Witkowski committed
194
      updateLocalGlobalNumbering();
195

Thomas Witkowski's avatar
Thomas Witkowski committed
196
197
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
198
#if (DEBUG != 0)
199
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
200
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], 
201
				    *(dofMaps[0]), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
202
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
203
#endif    
204
205

      initialized = true;
206
      return;
207
    }
208

Thomas Witkowski's avatar
Thomas Witkowski committed
209
210
211
    
    // Create a very first pariitioning of the mesh.
    createInitialPartitioning();
212

213

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

219
    if (mpiRank == 0) {
220
#if (DEBUG != 0)    
221
      int writePartMesh = 1;
222
223
224
#else
      int writePartMesh = 0;
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
225
      Parameters::get("parallel->debug->write mesh partitioning", writePartMesh);
226

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

233
    // Create interior boundary information.
234
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
235

236
237
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

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

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

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

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

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

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

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

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

271
272
273
	}
      }
    }
274

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

    removeMacroElements();
278

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

Thomas Witkowski's avatar
Thomas Witkowski committed
281

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

289
    updateLocalGlobalNumbering();
290

291
292
293
    // === 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.                         ===
294

295
296
    check3dValidMesh();

297

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

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

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

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

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

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

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

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

328
      updateLocalGlobalNumbering();
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

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
335
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
337
338
      elObjDb.clear();

339
    initialized = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
340
    MSG("Init parallelization needed %.5f seconds\n", MPI::Wtime() - first);
341
342
  }

343

Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
  void MeshDistributor::createInitialPartitioning()
  {
    FUNCNAME("MeshDistributor::createInitialPartitioning()");

    // 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();

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

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

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

    string filename = "";
    Parameters::get("parallel->partitioner->read meta arh", filename);
    if (filename == "" || ArhReader::readMetaData(filename) != mpiSize) {
      // and now partition the mesh    
      bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
      partitioner->createPartitionMap(partitionMap);
    }
  }


  void MeshDistributor::setInitialElementWeights() 
  {
    FUNCNAME("MeshDistributor::setInitialElementWeights()");

    elemWeights.clear();
      
    string filename = "";
    Parameters::get(mesh->getName() + "->macro weights", filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

      ifstream infile;
      infile.open(filename.c_str(), ifstream::in);
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;

	elemWeights[elNum] = elWeight;
      }

      infile.close();

      return;
    }

    filename = "";
    Parameters::get("parallel->partitioner->read meta arh", filename);
    if (filename != "") {
      map<int, int> arhElInRank;
      map<int, int> arhElCodeSize;

      int nProc = ArhReader::readMetaData(filename, arhElInRank, arhElCodeSize);
      for (map<int, int>::iterator it = arhElCodeSize.begin(); 
	   it != arhElCodeSize.end(); ++it)
	elemWeights[it->first] = it->second;

      return;
    }

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      elemWeights[elInfo->getElement()->getIndex()] = 1.0;
      elInfo = stack.traverseNext(elInfo);    
    }
  }


425
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
426
  {
427
    FUNCNAME("MeshDistributor::addProblemStat()");
428

429
430
431
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

432

433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
    // === Add all FE spaces from stationary problem. ===

    vector<const FiniteElemSpace*> newFeSpaces = probStat->getFeSpaces();
    for (int i = 0; i < static_cast<int>(newFeSpaces.size()); i++)
      if (find(feSpaces.begin(), feSpaces.end(), newFeSpaces[i]) == 
	  feSpaces.end())
	feSpaces.push_back(newFeSpaces[i]);
    

    // === Add mesh of stationary problem and create a corresponding  ===
    // === refinement manager object.                                 ===
    
    if (mesh != NULL) {
      TEST_EXIT(mesh == probStat->getMesh())
	("Does not yet support for different meshes!\n");
    } else {
      mesh = probStat->getMesh();
    }
451
452
453
454
455
456
457
458
459
460
    
    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());
461
    }
462
463
464
    
    partitioner->setMesh(mesh);
    
465

466
467
468
    // === Check whether the stationary problem should be serialized. ===


469
470
    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
471
    Parameters::get(probStat->getName() + "->output->write serialization",  
472
		    writeSerialization);
473
    if (writeSerialization && !writeSerializationFile) {
474
      string filename = "";
475
      Parameters::get(name + "->output->serialization filename", filename);
476
477
478
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
479
480

      int tsModulo = -1;
481
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
482
		      tsModulo);
483
      
484
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
485
486
      writeSerializationFile = true;
    }    
487

488
489
490

    // === Check whether the stationary problem should be deserialized. ===

491
    int readSerialization = 0;
492
    Parameters::get(probStat->getName() + "->input->read serialization", 
493
		    readSerialization);
494
    if (readSerialization) {
495
      string filename = "";
496
      Parameters::get(probStat->getName() + "->input->serialization filename", 
497
		      filename);
498
      filename += ".p" + lexical_cast<string>(mpiRank);
499
      MSG("Start deserialization with %s\n", filename.c_str());
500
      ifstream in(filename.c_str());
501
502
503
504

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

505
506
507
508
509
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

510
      probStat->deserialize(in);
511
      in.close();
512
513
      MSG("Deserialization from file: %s\n", filename.c_str());

514
      filename = "";
515
      Parameters::get(name + "->input->serialization filename", filename);
516
517
518
519
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
520
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
521
522
523
524
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
525
526
527
528
529

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
530
531
532
533
534
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
535
536
    }

537
    problemStat.push_back(probStat);
538

539
540
541

    // === If the mesh distributor is already initialized, don't forget to  ===
    // === remove the periodic boundary conditions on these objects.        ===
542

543
    if (initialized)
544
      removePeriodicBoundaryConditions(probStat);
545
546
547
  }


548
549
550
551
552
553
554
555
556
557
558
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


559
  void MeshDistributor::exitParallelization()
560
  {}
561

562
563
564
565
566
567
568
569
570
571
572
573

  void MeshDistributor::registerDofMap(ParallelDofMapping &dofMap)
  {
    FUNCNAME("MeshDistributor::registerDofMap()");

    TEST_EXIT(find(dofMaps.begin(), dofMaps.end(), &dofMap) ==
	      dofMaps.end())
      ("Parallel DOF mapping already registerd in mesh distributor object!\n");

    dofMaps.push_back(&dofMap);
  }

574
  
575
  void MeshDistributor::testForMacroMesh()
576
  {
577
    FUNCNAME("MeshDistributor::testForMacroMesh()");
578
579
580
581
582
583
584

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
589
590
591
592
593
594
595
596
597
598
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

599

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
600
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
601
  {
602
603
    FUNCNAME("MeshDistributor::synchVector()");

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
604
    TEST_EXIT_DBG(level >= 0 && level <= 1)("Wrong level number!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
605

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
606
607
608
609
610
611
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

    StdMpi<vector<double> > stdMpi(levelComm);

    for (DofComm::Iterator it(dc.getSendDofs()); 
612
	 !it.end(); it.nextRank()) {
613
      vector<double> dofs;
614

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

Thomas Witkowski's avatar
Thomas Witkowski committed
618
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
619
620
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
621

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
622
623
624
625
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
626
    }
627
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
628
629
630
631
632
633
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.recv(rank);
    }
634

635
    stdMpi.startCommunication();
636

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
637
638
639
640
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
641
642

      int counter = 0;
643
      vector<double> &recvData =  stdMpi.getRecvData(rank);
644
645
646
647

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

648
649
650
651
652
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof()) {
	  TEST_EXIT_DBG(counter < recvData.size())
	    ("Recv data from rank %d has only %d entries!\n", rank, recvData.size());
	  dofVec[it.getDofIndex()] = recvData[counter++];
	}
653
654
655
656
      }
    }
  }

657

658
659
660
661
662
663
664
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

665
666
667
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

668
669
670
671
672
673
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
674
      for (int i = 0; i < 6; i++) {
675
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
676
677
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
678
679
680
681
682
683
684
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

685

686
687
688
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

689
    FixRefinementPatch::connectedEdges.clear();
690
    bool valid3dMesh = true;
691

692
    // Traverse all edges
693
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
694
      // For each edge get all macro elements that contain this edge.
695
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
696
697
698
699

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

700
701
702
703
704
705
706

      // We have now to check, if the current edge connects two macro elements,
      // which are otherwise not connected. The basic idea to check this is very
      // simple: We take the first macro element in rank's subdomain that contain
      // this edge and add it to the set variable "el0". All other macro elements
      // which share this edge are added to "el1". 

707
      std::set<int> el0, el1;
708
      map<int, int> edgeNoInEl;
709
710
711

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
712
713
714
715
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
716
717
718
719

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
720
           
721
722
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

723
724
      // If el1 is empty, there is only on macro element in the mesh which
      // contains this edge. Hence, we can continue with checking another edge.
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
      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()) {
747
748
749
750
	valid3dMesh = false;

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

752
753
754
	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 = 
755
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
756
	    pair<Element*, int> edge1 = 
757
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
758

759
#if (DEBUG != 0)
760
761
762
763
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
764
765
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
766

767
768
769
770
771
772
773
774
775
	    MSG("FOUND EDGE %d/%d <-> %d/%d\n", 
		edge0.first->getIndex(), edge0.second,
		edge1.first->getIndex(), edge1.second);

	    MSG("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]);
#endif
776
777
778
779
780
781
782
783
784

	    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));
	  }
	}
      }
    }
785
786

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
787
788
789
  }


790
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
791
					   int level,
792
					   DofContainer& dofs)
793
794
795
796
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
797
798
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
799
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
800
801
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
802
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
803
804
805
806
807
808

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


809
810
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
811
812
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

813
    // Remove periodic boundaries in boundary manager on matrices and vectors.
814
815
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
816
817
818

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
819
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
820
821
822
823
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
824
825
826

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
827
828
829
  }


830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
  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()));
    }
  }


852
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
853
  {
854
855
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
856
857
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
858
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
859
	boundaryMap.erase(it++);
860
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
861
862
863
864
865
	++it;      
    }    
  }


866
867
868
869
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

870
871
872
    if (mpiSize != 16)
      return;

873
    std::set<int> neighbours;
874
875
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
876
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
877
878
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
879
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
880
881
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
882
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
883
884
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
885
886
887
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
888
889
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
890
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
891
892
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
893
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
894
895
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
896
897
898
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
899
900
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
901
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);
902
903
      break;
    case 8:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
904
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(9);  neighbours.insert(10); neighbours.insert(11);
905
906
      break;
    case 9:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
907
908
909
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(6); 
      neighbours.insert(8); neighbours.insert(12);
      neighbours.insert(10); neighbours.insert(11); neighbours.insert(14); 
910
911
      break;
    case 10:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
912
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(11);
913
914
      break;
    case 11:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
915
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12);  neighbours.insert(10); neighbours.insert(14);
916
917
      break;
    case 12:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
918
919
920
      neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); 
      neighbours.insert(9); neighbours.insert(13);
      neighbours.insert(11); neighbours.insert(14); neighbours.insert(15); 
921
922
      break;
    case 13:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
923
      neighbours.insert(6); neighbours.insert(7); neighbours.insert(12);  neighbours.insert(14); neighbours.insert(15);
924
925
      break;
    case 14:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
926
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13);  neighbours.insert(11); neighbours.insert(15);
927
928
      break;
    case 15:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
929
      neighbours.insert(12); neighbours.insert(13); neighbours.insert(14);
930
931
      break;
    }    
932

933
    TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
934
935
936
937
938
939
940

    levelData.init(neighbours);

    bool multiLevelTest = false;