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));
Thomas Witkowski's avatar
Thomas Witkowski committed
617
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
618
619
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
620

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
621
622
623
624
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
625
    }
626
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
627
628
629
630
631
632
    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);
    }
633

634
    stdMpi.startCommunication();
635

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

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

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

647
648
649
650
651
	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++];
	}
652
653
654
655
      }
    }
  }

656

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

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

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

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

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

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

      elInfo = stack.traverseNext(elInfo);
    }

684

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

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

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

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

699
700
701
702
703
704
705

      // 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". 

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

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

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

722
723
      // 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.
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
      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()) {
746
747
748
749
	valid3dMesh = false;

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

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

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

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

766
767
768
769
770
771
772
773
774
	    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
775
776
777
778
779
780
781
782
783

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

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


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

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

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


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

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

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

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


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


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

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


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

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

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

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

    levelData.init(neighbours);

    bool multiLevelTest = false;
    Parameters::get("parallel->multi level test", multiLevelTest);
    if (multiLevelTest) {
      switch (mpiRank) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
940
      case 0: case 1: case 2: case 3:
941
942
	{
	  std::set<int> m;