MeshDistributor.cc 68.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
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
89
      lastMeshChangeIndex(0),
90
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
91
      boundaryDofInfo(1),
92
      meshAdaptivity(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
93
      hasPeriodicBoundary(false),
94
95
      printTimings(false),
      printMemoryUsage(false)
96
  {
97
    FUNCNAME("MeshDistributor::MeshDistributor()");
Thomas Witkowski's avatar
Thomas Witkowski committed
98

99
    MPI::Intracomm &mpiComm = MPI::COMM_WORLD;
100
    mpiRank = mpiComm.Get_rank();
101

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

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

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

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

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

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

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

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

    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
139
140
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
141

Thomas Witkowski's avatar
Thomas Witkowski committed
142
143
144
145
146
147
148
  MeshDistributor::~MeshDistributor()
  {
    if (partitioner)
      delete partitioner;
  }


149
  void MeshDistributor::initParallelization()
150
  {
151
    FUNCNAME("MeshDistributor::initParallelization()");
152

153
154
155
    if (initialized)
      return;

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

160
    TEST_EXIT(MPI::COMM_WORLD.Get_size() > 1)
161
      ("Parallelization does not work with only one process!\n");
162
    TEST_EXIT(feSpaces.size() > 0)
163
      ("No FE space has been defined for the mesh distributor!\n");
164
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
165

166
167
168
169
170
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

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

183
    elObjDb.setFeSpace(feSpaces[0]);
184

Thomas Witkowski's avatar
Thomas Witkowski committed
185

186
187
188
    // 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).
189
    if (deserialized) {
190
191
      updateMacroElementInfo();

192
      removePeriodicBoundaryConditions();
193

194
      elObjDb.createMacroElementInfo(allMacroElements);
195

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

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

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

      initialized = true;
208
      return;
209
    }
210

211

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

215

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

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

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

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

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

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

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

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

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

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

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

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

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

273
274
275
	}
      }
    }
276

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

    removeMacroElements();
280

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

Thomas Witkowski's avatar
Thomas Witkowski committed
283

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

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
  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
Thomas Witkowski's avatar
Thomas Witkowski committed
360
    bool useInitialPartitioning = partitioner->createInitialPartitioning();
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
363
364

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

Thomas Witkowski's avatar
Thomas Witkowski committed
365
    if (!useInitialPartitioning) {
Thomas Witkowski's avatar
Thomas Witkowski committed
366
367
368
369
      // and now partition the mesh    
      bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
    }
370
371

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


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


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

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

433

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

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


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

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

489
490
491

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

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

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

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

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

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

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

538
    problemStat.push_back(probStat);
539

540
541
542

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

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


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

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

    globalMeshDistributor->addProblemStat(probStat);
  }


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

563
564
565
566
567
568
569
570
571

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

572
    dofMaps.push_back(&dofMap);   
573
574
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
575
576
577
578
579
580
581
582

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

    vector<ParallelDofMapping*>::iterator it = 
      find(dofMaps.begin(), dofMaps.end(), &dofMap);

583
584
    if (it != dofMaps.end())
      dofMaps.erase(it);
Thomas Witkowski's avatar
Thomas Witkowski committed
585
586
587
  }


588
  void MeshDistributor::testForMacroMesh()
589
  {
590
    FUNCNAME("MeshDistributor::testForMacroMesh()");
591
592
593
594
595
596
597

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
602
603
604
605
606
607
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

608
    TEST_EXIT(nMacroElements >= MPI::COMM_WORLD.Get_size())
609
610
611
      ("The mesh has less macro elements than number of mpi processes!\n");
  }

612

613
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
614
  {
615
    FUNCNAME("MeshDistributor::synchVector()");
616
617
618
619
620
621
622
   
    int nLevels = levelData.getNumberOfLevels();
    for (int level = nLevels - 1; level >= 0; level--) {
      MPI::Intracomm &mpiComm = levelData.getMpiComm(level);
      
      if (mpiComm == MPI::COMM_SELF)
	continue;
623

624
      StdMpi<vector<double> > stdMpi(mpiComm);
625

626
627
628
629
630
631
632
633
634
635
636
637
638
639
      for (DofComm::Iterator it(dofComm[level].getSendDofs()); 
	   !it.end(); it.nextRank()) {
	vector<double> dofs;
	
	for (int i = 0; i < vec.getSize(); i++) {
	  DOFVector<double> &dofVec = *(vec.getDOFVector(i));
	  for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
	    dofs.push_back(dofVec[it.getDofIndex()]);
	}
	
	int rank = it.getRank();
	// if (level > 0)
	//   rank = levelData.mapRank(rank, 0, level);
	stdMpi.send(rank, dofs);
640
      }
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
      
      for (DofComm::Iterator it(dofComm[level].getRecvDofs()); 
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	// if (level > 0)
	//   rank = levelData.mapRank(rank, 0, level);
	stdMpi.recv(rank);
      }
      
      stdMpi.startCommunication();
      
      for (DofComm::Iterator it(dofComm[level].getRecvDofs()); 
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	// if (level > 0)
	//   rank = levelData.mapRank(rank, 0, level);
	
	int counter = 0;
	vector<double> &recvData =  stdMpi.getRecvData(rank);
	
	for (int i = 0; i < vec.getSize(); i++) {
	  DOFVector<double> &dofVec = *(vec.getDOFVector(i));
	  
	  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++];
	  }
669
	}
670
671
672
673
      }
    }
  }

674

675
  void MeshDistributor::fix3dMeshRefinement()
676
  {
677
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
678
679
680
681

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

682
683
684
685
    MPI::COMM_WORLD.Barrier();
    double first = MPI::Wtime();


686
687
688
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

689
690
691
692
693
694
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
695
      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
696
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
697
698
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
699
700
701
702
703
704
705
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

706

707
708
    // === Reset fixed edges, and start new search for edges that must  ===
    // === be fixed.                                                    ===
709

710
    FixRefinementPatch::connectedEdges.clear();
711

712
    // Traverse all edges
713
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
714
      // For each edge get all macro elements that contain this edge.
715
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
716
717
718
719

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

720
721

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

724
725
726
      // All elements on this edge.
      std::set<int> allEls;
      // Maps from each element index to the local edge index of the common edge.
727
      map<int, int> edgeNoInEl;
728
729
730

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
731
	  allEls.insert(edgeEls[i].elIndex);
732
733
734
	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
735
           
736
      TEST_EXIT_DBG(!allEls.empty())("Should not happen!\n");
737

738
739
740
741

      // If there is only one element in rank on this edge, there is nothing 
      // to do.
      if (allEls.size() == 1)
742
743
	continue;
      
744
745
746
747
748
      // 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);   
749
      
750
751
752
753
754
755
756
757
758
759
760
761
      // 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) {
762
	    pair<Element*, int> edge1 = 
763
	      make_pair(elObjDb.getElementPtr(*elIt), edgeNoInEl[*elIt]);
764

765
#if (DEBUG != 0)
766
767
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);
768
	    
769
	    WorldVector<double> c0, c1;
770
771
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
772
773
774
775
776
777
778
779
780
781
782
783
	    
	    // 	      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
784

785
786
787
788
789
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}	  
      }
    }
790

791

792
793
794
    MPI::COMM_WORLD.Barrier();
    MSG("Fix 3D mesh refinement needed %.5f seconds\n", MPI::Wtime() - first);
  }
795

796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815

  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;
816
817
	  }
	}
818
819
820
	
	if (found)
	  break;
821
      }
822
    } while (found);
823

824
825
826
    disconnectedEls.push_back(firstElem);
    if (otherElems.size())
      helpToFix(otherElems, edgeNoInEl, disconnectedEls);
827
828
829
  }


830
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
831
					   int level,
832
					   DofContainer& dofs)
833
834
835
836
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
837
    for (DofComm::Iterator it(dofComm[level].getSendDofs(), feSpace); 
838
	 !it.end(); it.nextRank())
839
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
840
    for (DofComm::Iterator it(dofComm[level].getRecvDofs(), feSpace); 
841
	 !it.end(); it.nextRank())
842
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
843
844
845
846
847
848

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


849
850
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
851
852
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

853
    // Remove periodic boundaries in boundary manager on matrices and vectors.
854
855
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
856
857
858

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
859
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
860
861
862
863
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
864
865
866

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
867
868
869
  }


870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
  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()));
    }
  }


892
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
893
  {
894
895
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
896
897
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
898
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
899
	boundaryMap.erase(it++);
900
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
901
902
903
904
905
	++it;      
    }    
  }


906
907
  void MeshDistributor::createMeshLevelStructure()
  {
908
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");
909

910
    int levelMode = 0;
911
912
913
914
915
916
917
    Parameters::get("parallel->level mode", levelMode);

    TEST_EXIT(levelMode >= 0)("Wrong level mode %d!\n", levelMode);

    if (levelMode == 0) {
      std::set<int> neighbours;
      levelData.init(neighbours);
918
      return;
919
    }
920

921
922
923
924
925
926
    if (levelMode == 1) {
      std::set<int> neighbours;
      levelData.init(neighbours);
      levelData.addLevelMode1();
      return;
    }
927

928
929
    TEST_EXIT(MPI::COMM_WORLD.Get_size() == 16)
      ("Only mpiSize == 16 supported!\n");
930

931
    if (levelMode == 2) {
932

933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991