MeshDistributor.cc 68.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


Thomas Witkowski's avatar
Thomas Witkowski committed
22
#include <algorithm>
23
24
#include <iostream>
#include <fstream>
25
26
#include <limits>
#include <stdint.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
27
#include <boost/lexical_cast.hpp>
28
29
#include <boost/filesystem.hpp>

30
#include "parallel/MeshDistributor.h"
31
#include "parallel/MeshManipulation.h"
32
#include "parallel/ParallelDebug.h"
33
#include "parallel/StdMpi.h"
34
#include "parallel/MeshPartitioner.h"
35
#include "parallel/ParMetisPartitioner.h"
36
#include "parallel/ZoltanPartitioner.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
37
#include "parallel/SimplePartitioner.h"
38
#include "parallel/CheckerPartitioner.h"
39
#include "parallel/MpiHelper.h"
40
#include "parallel/DofComm.h"
41
#include "parallel/ParallelProblemStat.h"
42
43
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
44
#include "io/MacroWriter.h"
45
#include "io/VtkWriter.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
46
#include "io/ArhReader.h"
47
48
49
50
51
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
52
53
#include "DOFMatrix.h"
#include "DOFVector.h"
54
#include "SystemVector.h"
55
#include "ElementDofIterator.h"
56
57
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
58
#include "VertexVector.h"
59
#include "MeshStructure.h"
60
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
61
#include "ProblemInstat.h"
62
#include "RefinementManager3d.h"
63
#include "Debug.h"
64

65
namespace AMDiS { namespace Parallel {
66

Thomas Witkowski's avatar
Thomas Witkowski committed
67
  using boost::lexical_cast;
68
  using namespace boost::filesystem;
69
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
70

71
72
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

73
74
75
  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
76

77
78
79
80
81
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

82
  MeshDistributor::MeshDistributor()
83
    : problemStat(0),
84
      initialized(false),
85
      name("parallel"),
86
87
88
      mesh(NULL),
      refineManager(NULL),
      partitioner(NULL),
89
      initialPartitioner(NULL),
90
      deserialized(false),
91
      writeSerializationFile(false),
92
      repartitioningAllowed(false),
93
      repartitionIthChange(20),
94
      repartitioningWaitAfterFail(20),
95
96
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
97
      repartitioningFailed(0),
98
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
99
      lastMeshChangeIndex(0),
100
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
101
      boundaryDofInfo(1),
102
      meshAdaptivity(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
103
      hasPeriodicBoundary(false),
104
105
      printTimings(false),
      printMemoryUsage(false)
106
  {
107
    FUNCNAME("MeshDistributor::MeshDistributor()");
Thomas Witkowski's avatar
Thomas Witkowski committed
108

109
    MPI::Intracomm &mpiComm = MPI::COMM_WORLD;
110
    mpiRank = mpiComm.Get_rank();
111

112
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
113
114
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
115
    Parameters::get(name + "->repartition wait after fail", repartitioningWaitAfterFail);
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
116
    Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
117

Thomas Witkowski's avatar
Thomas Witkowski committed
118
    
119
120
    // === Create partitioner object. ===

121
    string partStr = "parmetis";
122
    Parameters::get(name + "->partitioner", partStr);
123
124

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

Thomas Witkowski's avatar
Thomas Witkowski committed
127
128
    if (partStr == "zoltan") {
#ifdef HAVE_ZOLTAN
129
      partitioner = new ZoltanPartitioner("parallel->partitioner", &mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
130
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
131
      ERROR_EXIT("AMDiS was compiled without Zoltan support. Therefore you cannot make use of it!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
132
133
#endif
    }
134

135
    if (partStr == "checker")
136
      partitioner = new CheckerPartitioner("parallel->partitioner", &mpiComm);
137

Thomas Witkowski's avatar
Thomas Witkowski committed
138
    if (partStr == "simple")
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
      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
159

160
    int tmp = 0;
161
    Parameters::get(name + "->box partitioning", tmp);
Thomas Witkowski's avatar
Thomas Witkowski committed
162
    partitioner->setBoxPartitioning(static_cast<bool>(tmp));
163
    initialPartitioner->setBoxPartitioning(static_cast<bool>(tmp));
Thomas Witkowski's avatar
Thomas Witkowski committed
164

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

168
169
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();
170
171
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
172

Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
  MeshDistributor::~MeshDistributor()
  {
175
    if (partitioner) {
Thomas Witkowski's avatar
Thomas Witkowski committed
176
      delete partitioner;
177
178
      partitioner = NULL;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
181
  }


182
  void MeshDistributor::initParallelization()
183
  {
184
    FUNCNAME("MeshDistributor::initParallelization()");
185

186
187
188
    if (initialized)
      return;

Thomas Witkowski's avatar
Thomas Witkowski committed
189
190
    double first = MPI::Wtime();
    MSG("Initialization phase 1 needed %.5f seconds\n", 
191
	first - ParallelProblemStat::initTimeStamp);
Thomas Witkowski's avatar
Thomas Witkowski committed
192

193
    TEST_EXIT(MPI::COMM_WORLD.Get_size() > 1)
194
      ("Parallelization does not work with only one process!\n");
195
    TEST_EXIT(feSpaces.size() > 0)
196
      ("No FE space has been defined for the mesh distributor!\n");
197
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
198

199
200
201
202
203
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
204
      doNext = false;
205
206
207
208
209
210
      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;
211
212
	  doNext = true;
	}
213
      }
214
215
    } while (doNext);

216
    elObjDb.setFeSpace(feSpaces[0]);
217

Thomas Witkowski's avatar
Thomas Witkowski committed
218

219
220
221
    // 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).
222
    if (deserialized) {
223
224
      updateMacroElementInfo();

225
      removePeriodicBoundaryConditions();
226

227
      elObjDb.createMacroElementInfo(allMacroElements);
228

Thomas Witkowski's avatar
Thomas Witkowski committed
229
      updateLocalGlobalNumbering();
230

Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
233
#if (DEBUG != 0)
234
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
235
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], 
236
				    *(dofMaps[0]), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
237
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
238
#endif    
239
240

      initialized = true;
241
      return;
242
    }
243

244

Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
    // Create a very first pariitioning of the mesh.
    createInitialPartitioning();
247

248

249
#if (DEBUG != 0)    
250
251
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
252
253
#endif

254
    if (mpiRank == 0) {
255
#if (DEBUG != 0)    
256
      int writePartMesh = 1;
257
258
259
#else
      int writePartMesh = 0;
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
260
      Parameters::get("parallel->debug->write mesh partitioning", writePartMesh);
261

262
      if (writePartMesh > 0) {
263
	debug::writeElementIndexMesh(mesh , debugOutputDir + "elementIndex");
264
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part");
265
      }
266
    }
267

268
    // Create interior boundary information.
269
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
270

271
272
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

273
274
275
276
277
    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))) {
278
279
280

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

281
282
283
284
285
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
286
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
287
288
289
290
291
292
293
294
295
	}
      }
    }

    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))) {
296
297
298

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

299
300
301
302
303
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

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

306
307
308
	}
      }
    }
309

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
310
311
312
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
313

314
315
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
316

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

324

325
    updateLocalGlobalNumbering();
326

327

328
329
330
    // === In 3D we have to fix the mesh to allow local refinements. ===

    fix3dMeshRefinement();
331

332

333
334
    // === If in debug mode, make some tests. ===

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

338
    ParallelDebug::testAllElements(*this);
339
    debug::testSortedDofs(mesh, elMap);
340
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
341
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
342

343
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
344
345

    MSG("Debug mode tests finished!\n");
346
#endif
347

348
349
350
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

351
352
353
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
354

355
    // === Global refinements. ===
356
    
Thomas Witkowski's avatar
Thomas Witkowski committed
357
    int globalRefinement = 0;
358
    Parameters::get(mesh->getName() + "->global refinements", globalRefinement);
359
    
Thomas Witkowski's avatar
Thomas Witkowski committed
360
    if (globalRefinement > 0) {
361
      refineManager->globalRefine(mesh, globalRefinement);
362
     
363
      updateLocalGlobalNumbering();
364

365
366
367
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
368
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
369

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
370
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
371
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
372
373
      elObjDb.clear();

374
    initialized = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
375
    MSG("Init parallelization needed %.5f seconds\n", MPI::Wtime() - first);
376
377
  }

378

Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
381
382
383
384
385
386
387
388
389
390
391
392
  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
393
394
    bool useInitialPartitioning = 
      initialPartitioner->createInitialPartitioning();
Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
397
398

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

399
    if (!useInitialPartitioning) {   
Thomas Witkowski's avatar
Thomas Witkowski committed
400
      // and now partition the mesh    
401
402
      bool partitioningSucceed = 
	initialPartitioner->partition(elemWeights, INITIAL);
Thomas Witkowski's avatar
Thomas Witkowski committed
403
404
      TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
    }
405

406
407
408
409
410
    initialPartitioner->createPartitionMap(partitionMap);

    if (initialPartitioner != partitioner) {
      *partitioner = *initialPartitioner;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
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
456
457
458
459
460
461
462
463
464
  }


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


465
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
466
  {
467
    FUNCNAME("MeshDistributor::addProblemStat()");
468

469
470
471
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

472

473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
    // === 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();
    }
491
492
493
494
495
496
497
498
499
500
    
    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());
501
    }
502
503
    
    partitioner->setMesh(mesh);
504
    initialPartitioner->setMesh(mesh);
505
    
506

507
508
509
    // === Check whether the stationary problem should be serialized. ===


510
511
    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
512
    Parameters::get(probStat->getName() + "->output->write serialization",  
513
		    writeSerialization);
514
    if (writeSerialization && !writeSerializationFile) {
515
      string filename = "";
516
      Parameters::get(name + "->output->serialization filename", filename);
517
518
519
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
520
521

      int tsModulo = -1;
522
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
523
		      tsModulo);
524
      
525
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
526
527
      writeSerializationFile = true;
    }    
528

529
530
531

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

532
    int readSerialization = 0;
533
    Parameters::get(probStat->getName() + "->input->read serialization", 
534
		    readSerialization);
535
    if (readSerialization) {
536
      string filename = "";
537
      Parameters::get(probStat->getName() + "->input->serialization filename", 
538
		      filename);
539
      filename += ".p" + lexical_cast<string>(mpiRank);
540
      MSG("Start deserialization with %s\n", filename.c_str());
541
      ifstream in(filename.c_str());
542
543
544
545

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

546
547
548
549
550
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

551
      probStat->deserialize(in);
552
      in.close();
553
554
      MSG("Deserialization from file: %s\n", filename.c_str());

555
      filename = "";
556
      Parameters::get(name + "->input->serialization filename", filename);
557
558
559
560
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
561
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
562
563
564
565
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
566
567
568
569
570

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
571
572
573
574
575
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
576
577
    }

578
    problemStat.push_back(probStat);
579

580
581
582

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

584
    if (initialized)
585
      removePeriodicBoundaryConditions(probStat);
586
587
588
  }


589
590
591
592
593
594
595
596
597
598
599
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


600
  void MeshDistributor::exitParallelization()
601
  {}
602

603
604
605
606
607
608
609
610
611

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

612
    dofMaps.push_back(&dofMap);   
613
614
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
615
616
617
618
619
620

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

621
622
    if (it != dofMaps.end())
      dofMaps.erase(it);
Thomas Witkowski's avatar
Thomas Witkowski committed
623
624
625
  }


626
  void MeshDistributor::testForMacroMesh()
627
  {
628
    FUNCNAME("MeshDistributor::testForMacroMesh()");
629
630
631
632
633
634
635

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
640
641
642
643
644
645
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

646
    TEST_EXIT(nMacroElements >= MPI::COMM_WORLD.Get_size())
647
648
649
      ("The mesh has less macro elements than number of mpi processes!\n");
  }

650

651
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
652
  {
653
    FUNCNAME("MeshDistributor::synchVector()");
654
655
656
657
658
659
660
   
    int nLevels = levelData.getNumberOfLevels();
    for (int level = nLevels - 1; level >= 0; level--) {
      MPI::Intracomm &mpiComm = levelData.getMpiComm(level);
      
      if (mpiComm == MPI::COMM_SELF)
	continue;
661

662
      StdMpi<vector<double> > stdMpi(mpiComm);
663

664
665
666
667
668
669
670
671
672
673
674
675
676
677
      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);
678
      }
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
      
      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++];
	  }
707
	}
708
709
710
711
      }
    }
  }

712

713
  void MeshDistributor::fix3dMeshRefinement()
714
  {
715
    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
716
717
718
719

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

720
721
722
723
    MPI::COMM_WORLD.Barrier();
    double first = MPI::Wtime();


724
725
726
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

727
728
729
730
731
732
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
733
      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
734
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
735
736
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
737
738
739
740
741
742
743
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

744

745
746
    // === Reset fixed edges, and start new search for edges that must  ===
    // === be fixed.                                                    ===
747

748
    FixRefinementPatch::connectedEdges.clear();
749

750
    // Traverse all edges
751
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
752
      // For each edge get all macro elements that contain this edge.
753
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
754
755
756
757

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

758
759

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

762
763
764
      // All elements on this edge.
      std::set<int> allEls;
      // Maps from each element index to the local edge index of the common edge.
765
      map<int, int> edgeNoInEl;
766
767
768

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
769
	  allEls.insert(edgeEls[i].elIndex);
770
771
772
	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
773
           
774
      TEST_EXIT_DBG(!allEls.empty())("Should not happen!\n");
775

776
777
778
779

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

803
#if (DEBUG != 0)
804
805
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);
806
	    
807
	    WorldVector<double> c0, c1;
808
809
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
810
811
812
813
814
815
816
817
818
819
820
821
	    
	    // 	      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
822

823
824
825
826
827
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}	  
      }
    }
828

829

830
831
832
    MPI::COMM_WORLD.Barrier();
    MSG("Fix 3D mesh refinement needed %.5f seconds\n", MPI::Wtime() - first);
  }
833

834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853

  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;
854
855
	  }
	}
856
857
858
	
	if (found)
	  break;
859
      }
860
    } while (found);
861

862
863
864
    disconnectedEls.push_back(firstElem);
    if (otherElems.size())
      helpToFix(otherElems, edgeNoInEl, disconnectedEls);
865
866
867
  }


868
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
869
					   int level,
870
					   DofContainer& dofs)
871
872
  {
    DofContainerSet dofSet;
873
    for (DofComm::Iterator it(dofComm[level].getSendDofs(), feSpace); 
874
	 !it.end(); it.nextRank())
875
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
876
    for (DofComm::Iterator it(dofComm[level].getRecvDofs(), feSpace); 
877
	 !it.end(); it.nextRank())
878
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
879
880
881
882
883
884

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


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

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

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
901
902
903
  }


904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
  void MeshDistributor::removePeriodicBoundaryConditions(ProblemStatSeq *probStat)
  {
    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()));
    }
  }


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


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