MeshDistributor.cc 68.9 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
      initialPartitioner(NULL),
81
      deserialized(false),
82
      writeSerializationFile(false),
83
      repartitioningAllowed(false),
84
      repartitionIthChange(20),
85
      repartitioningWaitAfterFail(20),
86
87
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
88
      repartitioningFailed(0),
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
      hasPeriodicBoundary(false),
95
96
      printTimings(false),
      printMemoryUsage(false)
97
  {
98
    FUNCNAME("MeshDistributor::MeshDistributor()");
Thomas Witkowski's avatar
Thomas Witkowski committed
99

100
    MPI::Intracomm &mpiComm = MPI::COMM_WORLD;
101
    mpiRank = mpiComm.Get_rank();
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);
108

Thomas Witkowski's avatar
Thomas Witkowski committed
109
    
110
111
    // === Create partitioner object. ===

112
    string partStr = "parmetis";
113
    Parameters::get(name + "->partitioner", partStr);
114
115

    if (partStr == "parmetis") 
116
      partitioner = new ParMetisPartitioner("parallel->partitioner", &mpiComm);
117

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

126
    if (partStr == "checker")
127
      partitioner = new CheckerPartitioner("parallel->partitioner", &mpiComm);
128

Thomas Witkowski's avatar
Thomas Witkowski committed
129
    if (partStr == "simple")
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
      partitioner = new SimplePartitioner("parallel->partitioner", &mpiComm);


    // === Create initial partitioner object. ===

    partStr = "";
    Parameters::get(name + "->initial partitioner", partStr);
    if (partStr == "") {
      initialPartitioner = partitioner;
    } else {    
      if (partStr == "checker") {
	initialPartitioner = 
	  new CheckerPartitioner("parallel->initial partitioner", &mpiComm);
      } else {
	ERROR_EXIT("Not yet supported, but very easy to implement!\n");
      }
    }


    // === And read some more parameters. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
150

151
    int tmp = 0;
152
    Parameters::get(name + "->box partitioning", tmp);
Thomas Witkowski's avatar
Thomas Witkowski committed
153
    partitioner->setBoxPartitioning(static_cast<bool>(tmp));
154
    initialPartitioner->setBoxPartitioning(static_cast<bool>(tmp));
Thomas Witkowski's avatar
Thomas Witkowski committed
155

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

159
160
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
161
162
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
163

Thomas Witkowski's avatar
Thomas Witkowski committed
164
165
  MeshDistributor::~MeshDistributor()
  {
166
    if (partitioner) {
Thomas Witkowski's avatar
Thomas Witkowski committed
167
      delete partitioner;
168
169
      partitioner = NULL;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
172
  }


173
  void MeshDistributor::initParallelization()
174
  {
175
    FUNCNAME("MeshDistributor::initParallelization()");
176

177
178
179
    if (initialized)
      return;

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

184
    TEST_EXIT(MPI::COMM_WORLD.Get_size() > 1)
185
      ("Parallelization does not work with only one process!\n");
186
    TEST_EXIT(feSpaces.size() > 0)
187
      ("No FE space has been defined for the mesh distributor!\n");
188
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
189

190
191
192
193
194
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
195
      doNext = false;
196
197
198
199
200
201
      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;
202
203
	  doNext = true;
	}
204
      }
205
206
    } while (doNext);

207
    elObjDb.setFeSpace(feSpaces[0]);
208

Thomas Witkowski's avatar
Thomas Witkowski committed
209

210
211
212
    // 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).
213
    if (deserialized) {
214
215
      updateMacroElementInfo();

216
      removePeriodicBoundaryConditions();
217

218
      elObjDb.createMacroElementInfo(allMacroElements);
219

Thomas Witkowski's avatar
Thomas Witkowski committed
220
      updateLocalGlobalNumbering();
221

Thomas Witkowski's avatar
Thomas Witkowski committed
222
223
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
224
#if (DEBUG != 0)
225
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
226
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], 
227
				    *(dofMaps[0]), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
228
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
229
#endif    
230
231

      initialized = true;
232
      return;
233
    }
234

235

Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
    // Create a very first pariitioning of the mesh.
    createInitialPartitioning();
238

239

240
#if (DEBUG != 0)    
241
242
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
243
244
#endif

245
    if (mpiRank == 0) {
246
#if (DEBUG != 0)    
247
      int writePartMesh = 1;
248
249
250
#else
      int writePartMesh = 0;
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
251
      Parameters::get("parallel->debug->write mesh partitioning", writePartMesh);
252

253
      if (writePartMesh > 0) {
254
	debug::writeElementIndexMesh(mesh , debugOutputDir + "elementIndex");
255
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part");
256
      }
257
    }
258

259
    // Create interior boundary information.
260
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
261

262
263
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

264
265
266
267
268
    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))) {
269
270
271

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

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

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
277
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
278
279
280
281
282
283
284
285
286
	}
      }
    }

    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))) {
287
288
289

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

290
291
292
293
294
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

297
298
299
	}
      }
    }
300

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
301
302
303
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
304

305
306
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
307

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

315

316
    updateLocalGlobalNumbering();
317

318

319
320
321
    // === In 3D we have to fix the mesh to allow local refinements. ===

    fix3dMeshRefinement();
322

323

324
325
    // === If in debug mode, make some tests. ===

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

329
    ParallelDebug::testAllElements(*this);
330
    debug::testSortedDofs(mesh, elMap);
331
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
332
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
333

334
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
335
336

    MSG("Debug mode tests finished!\n");
337
#endif
338

339
340
341
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

342
343
344
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
345

346
    // === Global refinements. ===
347
    
Thomas Witkowski's avatar
Thomas Witkowski committed
348
    int globalRefinement = 0;
349
    Parameters::get(mesh->getName() + "->global refinements", globalRefinement);
350
    
Thomas Witkowski's avatar
Thomas Witkowski committed
351
    if (globalRefinement > 0) {
352
      refineManager->globalRefine(mesh, globalRefinement);
353
     
354
      updateLocalGlobalNumbering();
355

356
357
358
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
359
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
360

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
361
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
362
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
363
364
      elObjDb.clear();

365
    initialized = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
366
    MSG("Init parallelization needed %.5f seconds\n", MPI::Wtime() - first);
367
368
  }

369

Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
372
373
374
375
376
377
378
379
380
381
382
383
  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
384
385
    bool useInitialPartitioning = 
      initialPartitioner->createInitialPartitioning();
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
389

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

390
    if (!useInitialPartitioning) {   
Thomas Witkowski's avatar
Thomas Witkowski committed
391
      // and now partition the mesh    
392
393
      bool partitioningSucceed = 
	initialPartitioner->partition(elemWeights, INITIAL);
Thomas Witkowski's avatar
Thomas Witkowski committed
394
395
      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
    }
396

397
398
399
400
401
    initialPartitioner->createPartitionMap(partitionMap);

    if (initialPartitioner != partitioner) {
      *partitioner = *initialPartitioner;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
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
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
  }


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


456
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
457
  {
458
    FUNCNAME("MeshDistributor::addProblemStat()");
459

460
461
462
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

463

464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
    // === 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();
    }
482
483
484
485
486
487
488
489
490
491
    
    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());
492
    }
493
494
    
    partitioner->setMesh(mesh);
495
    initialPartitioner->setMesh(mesh);
496
    
497

498
499
500
    // === Check whether the stationary problem should be serialized. ===


501
502
    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
503
    Parameters::get(probStat->getName() + "->output->write serialization",  
504
		    writeSerialization);
505
    if (writeSerialization && !writeSerializationFile) {
506
      string filename = "";
507
      Parameters::get(name + "->output->serialization filename", filename);
508
509
510
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
511
512

      int tsModulo = -1;
513
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
514
		      tsModulo);
515
      
516
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
517
518
      writeSerializationFile = true;
    }    
519

520
521
522

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

523
    int readSerialization = 0;
524
    Parameters::get(probStat->getName() + "->input->read serialization", 
525
		    readSerialization);
526
    if (readSerialization) {
527
      string filename = "";
528
      Parameters::get(probStat->getName() + "->input->serialization filename", 
529
		      filename);
530
      filename += ".p" + lexical_cast<string>(mpiRank);
531
      MSG("Start deserialization with %s\n", filename.c_str());
532
      ifstream in(filename.c_str());
533
534
535
536

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

537
538
539
540
541
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

542
      probStat->deserialize(in);
543
      in.close();
544
545
      MSG("Deserialization from file: %s\n", filename.c_str());

546
      filename = "";
547
      Parameters::get(name + "->input->serialization filename", filename);
548
549
550
551
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
552
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
553
554
555
556
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
557
558
559
560
561

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
562
563
564
565
566
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
567
568
    }

569
    problemStat.push_back(probStat);
570

571
572
573

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

575
    if (initialized)
576
      removePeriodicBoundaryConditions(probStat);
577
578
579
  }


580
581
582
583
584
585
586
587
588
589
590
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


591
  void MeshDistributor::exitParallelization()
592
  {}
593

594
595
596
597
598
599
600
601
602

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

603
    dofMaps.push_back(&dofMap);   
604
605
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
606
607
608
609
610
611
612
613

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

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

614
615
    if (it != dofMaps.end())
      dofMaps.erase(it);
Thomas Witkowski's avatar
Thomas Witkowski committed
616
617
618
  }


619
  void MeshDistributor::testForMacroMesh()
620
  {
621
    FUNCNAME("MeshDistributor::testForMacroMesh()");
622
623
624
625
626
627
628

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
633
634
635
636
637
638
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

639
    TEST_EXIT(nMacroElements >= MPI::COMM_WORLD.Get_size())
640
641
642
      ("The mesh has less macro elements than number of mpi processes!\n");
  }

643

644
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
645
  {
646
    FUNCNAME("MeshDistributor::synchVector()");
647
648
649
650
651
652
653
   
    int nLevels = levelData.getNumberOfLevels();
    for (int level = nLevels - 1; level >= 0; level--) {
      MPI::Intracomm &mpiComm = levelData.getMpiComm(level);
      
      if (mpiComm == MPI::COMM_SELF)
	continue;
654

655
      StdMpi<vector<double> > stdMpi(mpiComm);
656

657
658
659
660
661
662
663
664
665
666
667
668
669
670
      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);
671
      }
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
      
      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++];
	  }
700
	}
701
702
703
704
      }
    }
  }

705

706
  void MeshDistributor::fix3dMeshRefinement()
707
  {
708
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
709
710
711
712

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

713
714
715
716
    MPI::COMM_WORLD.Barrier();
    double first = MPI::Wtime();


717
718
719
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

720
721
722
723
724
725
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
726
      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
727
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
728
729
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
730
731
732
733
734
735
736
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

737

738
739
    // === Reset fixed edges, and start new search for edges that must  ===
    // === be fixed.                                                    ===
740

741
    FixRefinementPatch::connectedEdges.clear();
742

743
    // Traverse all edges
744
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
745
      // For each edge get all macro elements that contain this edge.
746
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
747
748
749
750

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

751
752

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

755
756
757
      // All elements on this edge.
      std::set<int> allEls;
      // Maps from each element index to the local edge index of the common edge.
758
      map<int, int> edgeNoInEl;
759
760
761

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
762
	  allEls.insert(edgeEls[i].elIndex);
763
764
765
	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
766
           
767
      TEST_EXIT_DBG(!allEls.empty())("Should not happen!\n");
768

769
770
771
772

      // If there is only one element in rank on this edge, there is nothing 
      // to do.
      if (allEls.size() == 1)
773
774
	continue;
      
775
776
777
778
779
      // 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);   
780
      
781
782
783
784
785
786
787
788
789
790
791
792
      // 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) {
793
	    pair<Element*, int> edge1 = 
794
	      make_pair(elObjDb.getElementPtr(*elIt), edgeNoInEl[*elIt]);
795

796
#if (DEBUG != 0)
797
798
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);
799
	    
800
	    WorldVector<double> c0, c1;
801
802
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
803
804
805
806
807
808
809
810
811
812
813
814
	    
	    // 	      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
815

816
817
818
819
820
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}	  
      }
    }
821

822

823
824
825
    MPI::COMM_WORLD.Barrier();
    MSG("Fix 3D mesh refinement needed %.5f seconds\n", MPI::Wtime() - first);
  }
826

827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846

  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;
847
848
	  }
	}
849
850
851
	
	if (found)
	  break;
852
      }
853
    } while (found);
854

855
856
857
    disconnectedEls.push_back(firstElem);
    if (otherElems.size())
      helpToFix(otherElems, edgeNoInEl, disconnectedEls);
858
859
860
  }


861
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
862
					   int level,
863
					   DofContainer& dofs)
864
865
866
867
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
868
    for (DofComm::Iterator it(dofComm[level].getSendDofs(), feSpace); 
869
	 !it.end(); it.nextRank())
870
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
871
    for (DofComm::Iterator it(dofComm[level].getRecvDofs(), feSpace); 
872
	 !it.end(); it.nextRank())
873
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
874
875
876
877
878
879

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


880
881
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
882
883
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

884
    // Remove periodic boundaries in boundary manager on matrices and vectors.
885
886
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
887
888
889

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
890
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
891
892
893
894
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
895
896
897

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
898
899
900
  }


901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
  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()));
    }
  }


923
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
924
  {
925
926
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
927
928
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
929
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
930
	boundaryMap.erase(it++);
931
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
932
933
934
935
936
	++it;      
    }    
  }


937
938
  void MeshDistributor::createMeshLevelStructure()
  {
939
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");