MeshDistributor.cc 63.8 KB
Newer Older
1
//
2
3
4
5
6
7
8
9
10
11
12
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include <algorithm>
14
15
#include <iostream>
#include <fstream>
16
17
#include <limits>
#include <stdint.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
18
#include <boost/lexical_cast.hpp>
19
20
#include <boost/filesystem.hpp>

21
#include "parallel/MeshDistributor.h"
22
#include "parallel/MeshManipulation.h"
23
#include "parallel/ParallelDebug.h"
24
#include "parallel/StdMpi.h"
25
#include "parallel/MeshPartitioner.h"
26
#include "parallel/ParMetisPartitioner.h"
27
#include "parallel/ZoltanPartitioner.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28
#include "parallel/SimplePartitioner.h"
29
#include "parallel/CheckerPartitioner.h"
30
#include "parallel/MpiHelper.h"
31
#include "parallel/DofComm.h"
32
33
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
34
#include "io/MacroWriter.h"
35
#include "io/VtkWriter.h"
36
37
38
39
40
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
41
42
#include "DOFMatrix.h"
#include "DOFVector.h"
43
#include "SystemVector.h"
44
#include "ElementDofIterator.h"
45
46
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
47
#include "VertexVector.h"
48
#include "MeshStructure.h"
49
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
50
#include "ProblemInstat.h"
51
#include "RefinementManager3d.h"
52
#include "Debug.h"
53

54
55
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
56
  using boost::lexical_cast;
57
  using namespace boost::filesystem;
58
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
59

60
61
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

62
63
64
  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
65

66
67
68
69
70
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

71
  bool MeshDistributor::sebastianMode = false;
72

73
  MeshDistributor::MeshDistributor()
74
    : problemStat(0),
75
      initialized(false),
76
      name("parallel"),
77
78
79
      mesh(NULL),
      refineManager(NULL),
      partitioner(NULL),
80
      deserialized(false),
81
      writeSerializationFile(false),
82
      repartitioningAllowed(false),
83
      repartitionIthChange(20),
84
85
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
86
      repartitioningFailed(0),
87
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
88
      lastMeshChangeIndex(0),
89
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
90
91
      boundaryDofInfo(1),
      meshAdaptivity(true)
92
  {
93
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
94

95
    mpiComm = MPI::COMM_WORLD;
96
97
    mpiRank = mpiComm.Get_rank();
    mpiSize = mpiComm.Get_size();
98

99
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
100
101
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
102
    Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
103

104
    string partStr = "parmetis";
105
    Parameters::get(name + "->partitioner", partStr);
106
107
108
109

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

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

118
119
120
    if (partStr == "checker")
      partitioner = new CheckerPartitioner(&mpiComm);

Thomas Witkowski's avatar
Thomas Witkowski committed
121
122
123
    if (partStr == "simple")
      partitioner = new SimplePartitioner(&mpiComm);

124
    int tmp = 0;
125
    Parameters::get(name + "->box partitioning", tmp);
Thomas Witkowski's avatar
Thomas Witkowski committed
126
127
    partitioner->setBoxPartitioning(static_cast<bool>(tmp));

128
    TEST_EXIT(partitioner)("Could not create partitioner \"%s\"!\n", partStr.c_str());
129
130
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
131

Thomas Witkowski's avatar
Thomas Witkowski committed
132
133
134
135
136
137
138
  MeshDistributor::~MeshDistributor()
  {
    if (partitioner)
      delete partitioner;
  }


139
  void MeshDistributor::initParallelization()
140
  {
141
    FUNCNAME("MeshDistributor::initParallelization()");
142

143
144
145
    if (initialized)
      return;

146
147
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");
148
    TEST_EXIT(feSpaces.size() > 0)
149
      ("No FE space has been defined for the mesh distributor!\n");
150
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
151

152
153
154
155
156
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
157
      doNext = false;
158
159
160
161
162
163
      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;
164
165
	  doNext = true;
	}
166
      }
167
168
    } while (doNext);

169
    elObjDb.setFeSpace(feSpaces[0]);
170

Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
173
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

174
175
176
    // 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).
177
    if (deserialized) {
Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
      createMeshLevelStructure();
      
180
181
      updateMacroElementInfo();

182
      removePeriodicBoundaryConditions();
183

184
      elObjDb.createMacroElementInfo(allMacroElements);
185

Thomas Witkowski's avatar
Thomas Witkowski committed
186
      updateLocalGlobalNumbering();
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188
189
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
190
#if (DEBUG != 0)
191
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
192
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], 
193
				    *(dofMaps[0]), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
194
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
195
#endif    
196
197

      initialized = true;
198
      return;
199
    }
200

201

202
203
204
    // 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.
205
206
    testForMacroMesh();

207
208
    // For later mesh repartitioning, we need to store some information about
    // the macro mesh.
209
210
    createMacroElementInfo();

211
    // create an initial partitioning of the mesh
212
    partitioner->createInitialPartitioning();
213

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

    // and now partition the mesh    
218
219
    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
220
    partitioner->createPartitionMap(partitionMap);
221

222

223
#if (DEBUG != 0)    
224
225
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
226
227
#endif

228
    if (mpiRank == 0) {
229
#if (DEBUG != 0)    
230
      int writePartMesh = 1;
231
232
233
#else
      int writePartMesh = 0;
#endif
234
      Parameters::get("dbg->write part mesh", writePartMesh);
235

236
      if (writePartMesh > 0) {
237
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex");
238
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
239
      }
240
    }
241

242
    // Create interior boundary information.
243
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
244

245
246
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

247
248
249
250
251
    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))) {
252
253
254

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

255
256
257
258
259
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
260
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
261
262
263
264
265
266
267
268
269
	}
      }
    }

    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))) {
270
271
272

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

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

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

280
281
282
	}
      }
    }
283

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
284
285
286
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
287

288
289
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
290

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

298
    updateLocalGlobalNumbering();
299

300
301
302
    // === In 3D we have to make some test, if the resulting mesh is valid.  ===
    // === If it is not valid, there is no possiblity yet to fix this        ===
    // === problem, just exit with an error message.                         ===
303

304
305
    check3dValidMesh();

306

307
308
    // === If in debug mode, make some tests. ===

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

312
    ParallelDebug::testAllElements(*this);
313
    debug::testSortedDofs(mesh, elMap);
314
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
315
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
316

317
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
318
319

    MSG("Debug mode tests finished!\n");
320
#endif
321

322
    // Create periodic DOF mapping, if there are periodic boundaries.
323
    createPeriodicMap();
324

325
326
327
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

328
329
330
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
331

332
    // === Global refinements. ===
333
    
Thomas Witkowski's avatar
Thomas Witkowski committed
334
    int globalRefinement = 0;
335
    Parameters::get(mesh->getName() + "->global refinements", globalRefinement);
336
    
Thomas Witkowski's avatar
Thomas Witkowski committed
337
    if (globalRefinement > 0) {
338
      refineManager->globalRefine(mesh, globalRefinement);
339

340
      updateLocalGlobalNumbering();
341

342
343
      // === Update periodic mapping, if there are periodic boundaries. ===     

344
      createPeriodicMap();
345
346
347
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
348
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
349

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
350
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
351
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
352
353
      elObjDb.clear();

354
    initialized = true;
355
356
  }

357

358
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
359
  {
360
    FUNCNAME("MeshDistributor::addProblemStat()");
361

362
363
364
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

365

366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
    // === 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();
    }
384
385
386
387
388
389
390
391
392
393
    
    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());
394
    }
395
396
397
    
    partitioner->setMesh(mesh);
    
398

399
400
401
    // === Check whether the stationary problem should be serialized. ===


402
403
    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
404
    Parameters::get(probStat->getName() + "->output->write serialization",  
405
		    writeSerialization);
406
    if (writeSerialization && !writeSerializationFile) {
407
      string filename = "";
408
      Parameters::get(name + "->output->serialization filename", filename);
409
410
411
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
412
413

      int tsModulo = -1;
414
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
415
		      tsModulo);
416
      
417
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
418
419
      writeSerializationFile = true;
    }    
420

421
422
423

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

424
    int readSerialization = 0;
425
    Parameters::get(probStat->getName() + "->input->read serialization", 
426
		    readSerialization);
427
    if (readSerialization) {
428
      string filename = "";
429
      Parameters::get(probStat->getName() + "->input->serialization filename", 
430
		      filename);
431
      filename += ".p" + lexical_cast<string>(mpiRank);
432
      MSG("Start deserialization with %s\n", filename.c_str());
433
      ifstream in(filename.c_str());
434
435
436
437

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

438
439
440
441
442
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

443
      probStat->deserialize(in);
444
      in.close();
445
446
      MSG("Deserialization from file: %s\n", filename.c_str());

447
      filename = "";
448
      Parameters::get(name + "->input->serialization filename", filename);
449
450
451
452
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
453
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
454
455
456
457
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
458
459
460
461
462

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
463
464
465
466
467
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
468
469
    }

470
    problemStat.push_back(probStat);
471

472
473
474

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

476
    if (initialized)
477
      removePeriodicBoundaryConditions(probStat);
478
479
480
  }


481
482
483
484
485
486
487
488
489
490
491
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


492
  void MeshDistributor::exitParallelization()
493
  {}
494

495
496
497
498
499
500
501
502
503
504
505
506

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

    TEST_EXIT(find(dofMaps.begin(), dofMaps.end(), &dofMap) ==
	      dofMaps.end())
      ("Parallel DOF mapping already registerd in mesh distributor object!\n");

    dofMaps.push_back(&dofMap);
  }

507
  
508
  void MeshDistributor::testForMacroMesh()
509
  {
510
    FUNCNAME("MeshDistributor::testForMacroMesh()");
511
512
513
514
515
516
517

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
522
523
524
525
526
527
528
529
530
531
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

532

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
533
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
534
  {
535
536
    FUNCNAME("MeshDistributor::synchVector()");

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
539
540
541
542
543
544
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

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

    for (DofComm::Iterator it(dc.getSendDofs()); 
545
	 !it.end(); it.nextRank()) {
546
      vector<double> dofs;
547

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

Thomas Witkowski's avatar
Thomas Witkowski committed
551
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
552
553
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
554

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
555
556
557
558
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
559
    }
560
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
561
562
563
564
565
566
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.recv(rank);
    }
567

568
    stdMpi.startCommunication();
569

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
570
571
572
573
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
574
575

      int counter = 0;
576
      vector<double> &recvData =  stdMpi.getRecvData(rank);
577
578
579
580

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

581
582
583
584
585
	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++];
	}
586
587
588
589
      }
    }
  }

590

591
592
593
594
595
596
597
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

598
599
600
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

601
602
603
604
605
606
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
607
      for (int i = 0; i < 6; i++) {
608
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
609
610
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
611
612
613
614
615
616
617
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

618

619
620
621
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

622
    FixRefinementPatch::connectedEdges.clear();
623
    bool valid3dMesh = true;
624

625
    // Traverse all edges
626
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
627
      // For each edge get all macro elements that contain this edge.
628
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
629
630
631
632

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

633
634
635
636
637
638
639

      // We have now to check, if the current edge connects two macro elements,
      // which are otherwise not connected. The basic idea to check this is very
      // simple: We take the first macro element in rank's subdomain that contain
      // this edge and add it to the set variable "el0". All other macro elements
      // which share this edge are added to "el1". 

640
      std::set<int> el0, el1;
641
      map<int, int> edgeNoInEl;
642
643
644

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
645
646
647
648
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
649
650
651
652

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
653
           
654
655
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

656
657
      // If el1 is empty, there is only on macro element in the mesh which
      // contains this edge. Hence, we can continue with checking another edge.
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
      if (el1.empty())
	continue;
      
      bool found = false;      
      do {
	found = false;
	for (std::set<int>::iterator elIt = el0.begin(); elIt != el0.end(); ++elIt) {
	  for (int i = 0; i < 4; i++) {
	    int neighEl = macroElementNeighbours[*elIt][i];
	    if (neighEl != -1 && el1.count(neighEl)) {
	      el0.insert(neighEl);
	      el1.erase(neighEl);
	      found = true;
	    }
	  }
	  
	  if (found)
	    break;
	}
      } while (found);
      
      if (!el1.empty()) {
680
681
682
683
	valid3dMesh = false;

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

685
686
687
	for (std::set<int>::iterator elIt0 = el0.begin(); elIt0 != el0.end(); ++elIt0) {
	  for (std::set<int>::iterator elIt1 = el1.begin(); elIt1 != el1.end(); ++elIt1) {
	    pair<Element*, int> edge0 = 
688
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
689
	    pair<Element*, int> edge1 = 
690
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
691

692
#if (DEBUG != 0)
693
694
695
696
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
697
698
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
699

700
701
702
703
704
705
706
707
708
	    MSG("FOUND EDGE %d/%d <-> %d/%d\n", 
		edge0.first->getIndex(), edge0.second,
		edge1.first->getIndex(), edge1.second);

	    MSG("FOUND EDGE: %d %d with coords %f %f %f   %f %f %f\n",
		dofEdge0.first, dofEdge0.second, 
		c0[0], c0[1], c0[2],
		c1[0], c1[1], c1[2]);
#endif
709
710
711
712
713
714
715
716
717

	    TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n");

	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge0, edge1));
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}
      }
    }
718
719

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
720
721
722
  }


723
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
724
					   int level,
725
					   DofContainer& dofs)
726
727
728
729
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
730
731
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
732
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
733
734
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
735
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
736
737
738
739
740
741

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


742
743
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
744
745
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

746
    // Remove periodic boundaries in boundary manager on matrices and vectors.
747
748
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
749
750
751

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
752
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
753
754
755
756
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
757
758
759

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
760
761
762
  }


763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
  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()));
    }
  }


785
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
786
  {
787
788
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
789
790
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
791
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
792
	boundaryMap.erase(it++);
793
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
794
795
796
797
798
	++it;      
    }    
  }


799
800
801
802
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

803
804
805
    if (mpiSize != 16)
      return;

806
    std::set<int> neighbours;
807
808
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
809
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
810
811
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
812
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
813
814
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
815
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
816
817
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
818
819
820
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
821
822
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
823
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
824
825
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
826
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
827
828
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
829
830
831
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
832
833
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
834
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);
835
836
      break;
    case 8:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
837
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(9);  neighbours.insert(10); neighbours.insert(11);
838
839
      break;
    case 9:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
840
841
842
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(6); 
      neighbours.insert(8); neighbours.insert(12);
      neighbours.insert(10); neighbours.insert(11); neighbours.insert(14); 
843
844
      break;
    case 10:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
845
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(11);
846
847
      break;
    case 11:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
848
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12);  neighbours.insert(10); neighbours.insert(14);
849
850
      break;
    case 12:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
851
852
853
      neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); 
      neighbours.insert(9); neighbours.insert(13);
      neighbours.insert(11); neighbours.insert(14); neighbours.insert(15); 
854
855
      break;
    case 13:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
856
      neighbours.insert(6); neighbours.insert(7); neighbours.insert(12);  neighbours.insert(14); neighbours.insert(15);
857
858
      break;
    case 14:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
859
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13);  neighbours.insert(11); neighbours.insert(15);
860
861
      break;
    case 15:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
862
      neighbours.insert(12); neighbours.insert(13); neighbours.insert(14);
863
864
      break;
    }    
865

866
    TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
867
868
869
870
871
872
873

    levelData.init(neighbours);

    bool multiLevelTest = false;
    Parameters::get("parallel->multi level test", multiLevelTest);
    if (multiLevelTest) {
      switch (mpiRank) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
874
      case 0: case 1: case 2: case 3:
875
876
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
877
	  m.insert(0); m.insert(1); m.insert(2); m.insert(3);
878
	  levelData.addLevel(m, 0);
Thomas Witkowski's avatar