MeshDistributor.cc 64.7 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
  MeshDistributor::MeshDistributor()
72
    : problemStat(0),
73
      initialized(false),
74
      name("parallel"),
75
76
77
      mesh(NULL),
      refineManager(NULL),
      partitioner(NULL),
78
      deserialized(false),
79
      writeSerializationFile(false),
80
      repartitioningAllowed(false),
81
      repartitionIthChange(20),
82
83
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
84
      repartitioningFailed(0),
85
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
86
      lastMeshChangeIndex(0),
87
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
88
      boundaryDofInfo(1),
89
      meshAdaptivity(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
90
91
      hasPeriodicBoundary(false),
      printTimings(false)
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);
Thomas Witkowski's avatar
Thomas Witkowski committed
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));

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
133

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


141
  void MeshDistributor::initParallelization()
142
  {
143
    FUNCNAME("MeshDistributor::initParallelization()");
144

145
146
147
    if (initialized)
      return;

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

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

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

171
    elObjDb.setFeSpace(feSpaces[0]);
172

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

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

184
      removePeriodicBoundaryConditions();
185

186
      elObjDb.createMacroElementInfo(allMacroElements);
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188
      updateLocalGlobalNumbering();
189

Thomas Witkowski's avatar
Thomas Witkowski committed
190
191
      elObjDb.setData(partitionMap, levelData);

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

      initialized = true;
200
      return;
201
    }
202

203

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

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

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

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

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

224

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

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

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

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

247
248
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

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

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

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

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

    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))) {
272
273
274

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

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

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

282
283
284
	}
      }
    }
285

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

    removeMacroElements();
289

290
291
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
292

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

300
    updateLocalGlobalNumbering();
301

302
303
304
    // === 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.                         ===
305

306
307
    check3dValidMesh();

308

309
310
    // === If in debug mode, make some tests. ===

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

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

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

    MSG("Debug mode tests finished!\n");
322
#endif
323

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

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

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

339
      updateLocalGlobalNumbering();
340

341
342
343
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
344
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
345

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
346
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
347
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
348
349
      elObjDb.clear();

350
    initialized = true;
351
352
  }

353

354
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
355
  {
356
    FUNCNAME("MeshDistributor::addProblemStat()");
357

358
359
360
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

361

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

395
396
397
    // === Check whether the stationary problem should be serialized. ===


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

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

417
418
419

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

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

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

434
435
436
437
438
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

439
      probStat->deserialize(in);
440
      in.close();
441
442
      MSG("Deserialization from file: %s\n", filename.c_str());

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

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

466
    problemStat.push_back(probStat);
467

468
469
470

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

472
    if (initialized)
473
      removePeriodicBoundaryConditions(probStat);
474
475
476
  }


477
478
479
480
481
482
483
484
485
486
487
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

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

    globalMeshDistributor->addProblemStat(probStat);
  }


488
  void MeshDistributor::exitParallelization()
489
  {}
490

491
492
493
494
495
496
497
498
499
500
501
502

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

503
  
504
  void MeshDistributor::testForMacroMesh()
505
  {
506
    FUNCNAME("MeshDistributor::testForMacroMesh()");
507
508
509
510
511
512
513

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
518
519
520
521
522
523
524
525
526
527
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

528

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
529
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
530
  {
531
532
    FUNCNAME("MeshDistributor::synchVector()");

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
535
536
537
538
539
540
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

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

    for (DofComm::Iterator it(dc.getSendDofs()); 
541
	 !it.end(); it.nextRank()) {
542
      vector<double> dofs;
543

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

Thomas Witkowski's avatar
Thomas Witkowski committed
547
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
548
549
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
550

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
551
552
553
554
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
555
    }
556
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
557
558
559
560
561
562
    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);
    }
563

564
    stdMpi.startCommunication();
565

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

      int counter = 0;
572
      vector<double> &recvData =  stdMpi.getRecvData(rank);
573
574
575
576

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

577
578
579
580
581
	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++];
	}
582
583
584
585
      }
    }
  }

586

587
588
589
590
591
592
593
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

594
595
596
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

597
598
599
600
601
602
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

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

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

      elInfo = stack.traverseNext(elInfo);
    }

614

615
616
617
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

618
    FixRefinementPatch::connectedEdges.clear();
619
    bool valid3dMesh = true;
620

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

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

629
630
631
632
633
634
635

      // 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". 

636
      std::set<int> el0, el1;
637
      map<int, int> edgeNoInEl;
638
639
640

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
641
642
643
644
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
645
646
647
648

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
649
           
650
651
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

652
653
      // 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.
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
      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()) {
676
677
678
679
	valid3dMesh = false;

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

681
682
683
	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 = 
684
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
685
	    pair<Element*, int> edge1 = 
686
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
687

688
#if (DEBUG != 0)
689
690
691
692
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
693
694
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
695

696
697
698
699
700
701
702
703
704
	    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
705
706
707
708
709
710
711
712
713

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

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
716
717
718
  }


719
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
720
					   int level,
721
					   DofContainer& dofs)
722
723
724
725
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
726
727
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
728
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
729
730
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
731
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
732
733
734
735
736
737

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


738
739
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
740
741
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

742
    // Remove periodic boundaries in boundary manager on matrices and vectors.
743
744
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
745
746
747

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

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
756
757
758
  }


759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
  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()));
    }
  }


781
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
782
  {
783
784
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

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


795
796
797
798
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

799
800
801
    if (mpiSize != 16)
      return;

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

862
    TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
863
864
865
866
867
868
869

    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
870
      case 0: case 1: case 2: case 3:
871
872
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
873
	  m.insert(0); m.insert(1); m.insert(2); m.insert(3);
874
	  levelData.addLevel(m, 0);
875
876
	}
	break;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
877
      case 4: case 5: case 6: case 7:
878
879
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
880
	  m.insert(4); m.