MeshDistributor.cc 67.3 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
      repartitioningWaitAfterFail(20),
85
86
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
87
      repartitioningFailed(0),
88

89
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
90
      lastMeshChangeIndex(0),
91
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
92
      boundaryDofInfo(1),
93
      meshAdaptivity(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
94
95
      hasPeriodicBoundary(false),
      printTimings(false)
96
  {
97
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
98

99
    mpiComm = MPI::COMM_WORLD;
100
101
    mpiRank = mpiComm.Get_rank();
    mpiSize = mpiComm.Get_size();
102

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
138

Thomas Witkowski's avatar
Thomas Witkowski committed
139
140
141
142
143
144
145
  MeshDistributor::~MeshDistributor()
  {
    if (partitioner)
      delete partitioner;
  }


146
  void MeshDistributor::initParallelization()
147
  {
148
    FUNCNAME("MeshDistributor::initParallelization()");
149

150
151
152
    if (initialized)
      return;

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

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

163
164
165
166
167
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
168
      doNext = false;
169
170
171
172
173
174
      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;
175
176
	  doNext = true;
	}
177
      }
178
179
    } while (doNext);

180
    elObjDb.setFeSpace(feSpaces[0]);
181

Thomas Witkowski's avatar
Thomas Witkowski committed
182
183
184
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

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

193
      removePeriodicBoundaryConditions();
194

195
      elObjDb.createMacroElementInfo(allMacroElements);
196

Thomas Witkowski's avatar
Thomas Witkowski committed
197
      updateLocalGlobalNumbering();
198

Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
      elObjDb.setData(partitionMap, levelData);

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

      initialized = true;
209
      return;
210
    }
211

212

Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
    // Create a very first pariitioning of the mesh.
    createInitialPartitioning();
215

216

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

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

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

236
    // Create interior boundary information.
237
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
238

239
240
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

241
242
243
244
245
    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))) {
246
247
248

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

249
250
251
252
253
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

    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))) {
264
265
266

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

267
268
269
270
271
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

274
275
276
	}
      }
    }
277

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
278
279
280
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
281

282
283
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
284

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

292
    updateLocalGlobalNumbering();
293

294

295
296
297
    // === In 3D we have to fix the mesh to allow local refinements. ===

    fix3dMeshRefinement();
298

299

300
301
    // === If in debug mode, make some tests. ===

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

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

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

    MSG("Debug mode tests finished!\n");
313
#endif
314

315
316
317
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

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

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

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

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

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

345

Thomas Witkowski's avatar
Thomas Witkowski committed
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
  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");
    }
372
373

    partitioner->createPartitionMap(partitionMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
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
425
426
427
  }


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


428
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
429
  {
430
    FUNCNAME("MeshDistributor::addProblemStat()");
431

432
433
434
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

435

436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
    // === 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();
    }
454
455
456
457
458
459
460
461
462
463
    
    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());
464
    }
465
466
467
    
    partitioner->setMesh(mesh);
    
468

469
470
471
    // === Check whether the stationary problem should be serialized. ===


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

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

491
492
493

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

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

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

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

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

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

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

540
    problemStat.push_back(probStat);
541

542
543
544

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

546
    if (initialized)
547
      removePeriodicBoundaryConditions(probStat);
548
549
550
  }


551
552
553
554
555
556
557
558
559
560
561
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


562
  void MeshDistributor::exitParallelization()
563
  {}
564

565
566
567
568
569
570
571
572
573
574
575
576

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

577
  
578
  void MeshDistributor::testForMacroMesh()
579
  {
580
    FUNCNAME("MeshDistributor::testForMacroMesh()");
581
582
583
584
585
586
587

    int nMacroElements = 0;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

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

602

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

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

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

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

    for (DofComm::Iterator it(dc.getSendDofs()); 
615
	 !it.end(); it.nextRank()) {
616
      vector<double> dofs;
617

618
619
      for (int i = 0; i < vec.getSize(); i++) {
	DOFVector<double> &dofVec = *(vec.getDOFVector(i));
Thomas Witkowski's avatar
Thomas Witkowski committed
620
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
621
622
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
623

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

637
    stdMpi.startCommunication();
638

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

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

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

650
651
652
653
654
	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++];
	}
655
656
657
658
      }
    }
  }

659

660
  void MeshDistributor::fix3dMeshRefinement()
661
  {
662
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
663
664
665
666

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

667
668
669
670
    MPI::COMM_WORLD.Barrier();
    double first = MPI::Wtime();


671
672
673
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

674
675
676
677
678
679
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
680
      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
681
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
682
683
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
684
685
686
687
688
689
690
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

691

692
693
    // === Reset fixed edges, and start new search for edges that must  ===
    // === be fixed.                                                    ===
694

695
    FixRefinementPatch::connectedEdges.clear();
696

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

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

705
706

      // We have now to check, if the current edge connects two macro elements,
707
      // which are otherwise not connected. 
708

709
710
711
      // All elements on this edge.
      std::set<int> allEls;
      // Maps from each element index to the local edge index of the common edge.
712
      map<int, int> edgeNoInEl;
713
714
715

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
716
	  allEls.insert(edgeEls[i].elIndex);
717
718
719
	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
720
           
721
      TEST_EXIT_DBG(!allEls.empty())("Should not happen!\n");
722

723
724
725
726

      // If there is only one element in rank on this edge, there is nothing 
      // to do.
      if (allEls.size() == 1)
727
728
	continue;
      
729
730
731
732
733
      // Create a set of element index sets. All element within one set are 
      // neighbours of each other, but can not reach any other element of the
      // other index sets by neighbourhood relations.
      vector<std::set<int> > disconnectedEls;
      helpToFix(allEls, edgeNoInEl, disconnectedEls);   
734
      
735
736
737
738
739
740
741
742
743
744
745
746
      // Fix the edges
      for (int dc0 = 0; dc0 < static_cast<int>(disconnectedEls.size()); dc0++) {
	for (int dc1 = 0; dc1 < static_cast<int>(disconnectedEls.size()); dc1++) {
	  if (dc0 == dc1)
	    continue;   
	  
	  int elIdx = *(disconnectedEls[dc1].begin());
	  pair<Element*, int> edge0 = 
	    make_pair(elObjDb.getElementPtr(elIdx), edgeNoInEl[elIdx]);
	  
	  for (std::set<int>::iterator elIt = disconnectedEls[dc0].begin();
	       elIt != disconnectedEls[dc0].end(); ++elIt) {
747
	    pair<Element*, int> edge1 = 
748
	      make_pair(elObjDb.getElementPtr(*elIt), edgeNoInEl[*elIt]);
749

750
#if (DEBUG != 0)
751
752
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);
753
	    
754
	    WorldVector<double> c0, c1;
755
756
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
757
758
759
760
761
762
763
764
765
766
767
768
	    
	    // 	      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]);
	    
	    TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n");
#endif
769

770
771
772
773
774
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}	  
      }
    }
775

776

777
778
779
    MPI::COMM_WORLD.Barrier();
    MSG("Fix 3D mesh refinement needed %.5f seconds\n", MPI::Wtime() - first);
  }
780

781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800

  void MeshDistributor::helpToFix(std::set<int> &elems, 
				  map<int, int> &edgeNoInEl,
				  vector<std::set<int> > &disconnectedEls)
  {
    std::set<int> firstElem;
    firstElem.insert(*(elems.begin()));

    std::set<int> otherElems(++(elems.begin()), elems.end());
    
    bool found = false;      
    do {
      found = false;
      for (std::set<int>::iterator elIt = firstElem.begin(); elIt != firstElem.end(); ++elIt) {
	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
	  int neighEl = macroElementNeighbours[*elIt][i];
	  if (neighEl != -1 && otherElems.count(neighEl)) {
	    firstElem.insert(neighEl);
	    otherElems.erase(neighEl);
	    found = true;
801
802
	  }
	}
803
804
805
	
	if (found)
	  break;
806
      }
807
    } while (found);
808

809
810
811
    disconnectedEls.push_back(firstElem);
    if (otherElems.size())
      helpToFix(otherElems, edgeNoInEl, disconnectedEls);
812
813
814
  }


815
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
816
					   int level,
817
					   DofContainer& dofs)
818
819
820
821
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
822
823
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
824
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
825
826
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
827
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
828
829
830
831
832
833

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


834
835
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
836
837
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

838
    // Remove periodic boundaries in boundary manager on matrices and vectors.
839
840
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
841
842
843

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
844
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
845
846
847
848
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
849
850
851

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
852
853
854
  }


855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
  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()));
    }
  }


877
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
878
  {
879
880
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
881
882
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
883
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
884
	boundaryMap.erase(it++);
885
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
886
887
888
889
890
	++it;      
    }    
  }


891
892
893
894
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

895
896
897
    if (mpiSize != 16)
      return;

898
    std::set<int> neighbours;
899
900
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
901
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
902
903
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
904
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
905
906
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
907
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
908
909
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
910
911
912
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
913
914
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
915
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
916
917
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
918
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
919
920
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
921
922
923
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
924
925
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
926
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);
927
928
      break;
    case 8:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
929
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(9);  neighbours.insert(10); neighbours.insert(11);