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


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

21
#include "parallel/MeshDistributor.h"
22
#include "parallel/MeshManipulation.h"
23
#include "parallel/ParallelDebug.h"
24
#include "parallel/StdMpi.h"
25
#include "parallel/MeshPartitioner.h"
26
#include "parallel/ParMetisPartitioner.h"
27
#include "parallel/ZoltanPartitioner.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28
#include "parallel/SimplePartitioner.h"
29
#include "parallel/CheckerPartitioner.h"
30
#include "parallel/MpiHelper.h"
31
#include "parallel/DofComm.h"
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
      feSpaces(0),
78
79
80
81
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
82
      deserialized(false),
83
      writeSerializationFile(false),
84
      repartitioningAllowed(false),
85
      repartitionIthChange(20),
86
87
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
88
      repartitioningFailed(0),
89
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
90
      lastMeshChangeIndex(0),
91
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
92
93
      boundaryDofInfo(1),
      meshAdaptivity(true)
94
  {
95
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
96

97
    mpiComm = MPI::COMM_WORLD;
98
99
    mpiRank = mpiComm.Get_rank();
    mpiSize = mpiComm.Get_size();
100

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

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

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

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

120
121
122
    if (partStr == "checker")
      partitioner = new CheckerPartitioner(&mpiComm);

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

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

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
150
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");

151
152
    TEST_EXIT(feSpaces.size() > 0)
      ("No FE space has been defined for the mesh distributor!\n");
153
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
154

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

    bool doNext = false;
    do {
160
161
      doNext = false;
      for (unsigned int i = 0; i < feSpaces.size() - 1; i++) {
162
163
164
165
166
167
168
	if (feSpaces[i]->getBasisFcts()->getDegree() >
	    feSpaces[i + 1]->getBasisFcts()->getDegree()) {
	  const FiniteElemSpace *tmp = feSpaces[i + 1];
	  feSpaces[i + 1] = feSpaces[i];
	  feSpaces[i] = tmp;
	  doNext = true;
	}
169
      }
170
171
    } while (doNext);

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

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

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

185
      removePeriodicBoundaryConditions();
186

187
      elObjDb.createMacroElementInfo(allMacroElements);
188

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
193
#if (DEBUG != 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
194
195
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], dofMap, 
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
196
#endif    
197
198

      initialized = true;
199
      return;
200
    }
201

202

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

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

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

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

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

223

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

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

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

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

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

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

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

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

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

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

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

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

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

281
282
283
	}
      }
    }
284

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

    removeMacroElements();
288

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

Thomas Witkowski's avatar
Thomas Witkowski committed
291

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

299
    updateLocalGlobalNumbering();
300

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

305
306
    check3dValidMesh();

307

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

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

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

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

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

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

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

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

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

341
      updateLocalGlobalNumbering();
342

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

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

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

355
    initialized = true;
356
357
  }

358

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

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

366
367
368
369
370
371
372
373
374
375
    // === Add FE spaces from stationary problem to mesh distributor. ===

    for (unsigned int i = 0; i < probStat->getFeSpaces().size(); i++) {
      bool foundFeSpace = false;
      for (unsigned int j = 0; j < feSpaces.size(); j++) {
	if (feSpaces[j] == probStat->getFeSpaces()[i])
	  foundFeSpace = true;

	TEST_EXIT(feSpaces[j]->getMesh() == probStat->getFeSpaces()[i]->getMesh())
	  ("MeshDistributor does not yet support usage of different meshes!\n");
376
377
      }

378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
      if (!foundFeSpace)
	feSpaces.push_back(probStat->getFeSpaces()[i]);
    }

    TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n");

    mesh = feSpaces[0]->getMesh();
    info = probStat->getInfo();
    
    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());
396
    }
397
398
399
    
    partitioner->setMesh(mesh);
    
400
401
402

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

      int tsModulo = -1;
413
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
414
		      tsModulo);
415
      
416
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
417
418
      writeSerializationFile = true;
    }    
419
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 set rank
    // DOFs object to the matrices and vectors of the added stationary problem,
    // and 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
  void MeshDistributor::testForMacroMesh()
493
  {
494
    FUNCNAME("MeshDistributor::testForMacroMesh()");
495
496
497
498
499
500
501

    int nMacroElements = 0;

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

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
506
507
508
509
510
511
512
513
514
515
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

516

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
517
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
518
  {
519
520
    FUNCNAME("MeshDistributor::synchVector()");

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
523
524
525
526
527
528
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

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

    for (DofComm::Iterator it(dc.getSendDofs()); 
529
	 !it.end(); it.nextRank()) {
530
      vector<double> dofs;
531

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

Thomas Witkowski's avatar
Thomas Witkowski committed
535
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
536
537
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
538

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
539
540
541
542
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
543
    }
544
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
545
546
547
548
549
550
    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);
    }
551

552
    stdMpi.startCommunication();
553

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
554
555
556
557
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
558
559

      int counter = 0;
560
      vector<double> &recvData =  stdMpi.getRecvData(rank);
561
562
563
564

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

565
566
567
568
569
	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++];
	}
570
571
572
573
      }
    }
  }

574

575
576
577
578
579
580
581
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

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

582
583
584
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

585
586
587
588
589
590
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
591
      for (int i = 0; i < 6; i++) {
592
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
593
594
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
595
596
597
598
599
600
601
      }

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

      elInfo = stack.traverseNext(elInfo);
    }

602

603
604
605
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

606
    FixRefinementPatch::connectedEdges.clear();
607
    bool valid3dMesh = true;
608

609
    // Traverse all edges
610
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
611
      // For each edge get all macro elements that contain this edge.
612
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
613
614
615
616

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

617
618
619
620
621
622
623

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

624
      std::set<int> el0, el1;
625
      map<int, int> edgeNoInEl;
626
627
628

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
629
630
631
632
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
633
634
635
636

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
637
           
638
639
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

640
641
      // 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.
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
      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()) {
664
665
666
667
	valid3dMesh = false;

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

669
670
671
	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 = 
672
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
673
	    pair<Element*, int> edge1 = 
674
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
675

676
#if (DEBUG != 0)
677
678
679
680
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
681
682
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
683

684
685
686
687
688
689
690
691
692
	    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
693
694
695
696
697
698
699
700
701

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

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
704
705
706
  }


707
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
708
					   int level,
709
					   DofContainer& dofs)
710
711
712
713
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
714
715
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
716
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
717
718
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
719
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
720
721
722
723
724
725

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


726
727
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
728
729
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

730
    // Remove periodic boundaries in boundary manager on matrices and vectors.
731
732
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
733
734
735

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
736
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
737
738
739
740
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
741
742
743

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
744
745
746
  }


747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
  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()));
    }
  }


769
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
770
  {
771
772
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
773
774
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
775
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
776
	boundaryMap.erase(it++);
777
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
778
779
780
781
782
	++it;      
    }    
  }


783
784
785
786
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

787
788
789
    if (mpiSize != 16)
      return;

790
    std::set<int> neighbours;
791
792
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
793
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
794
795
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
796
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
797
798
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
799
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
800
801
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
802
803
804
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
805
806
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
807
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
808
809
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
810
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
811
812
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
813
814
815
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
816
817
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
818
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);
819
820
      break;
    case 8:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
821
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(9);  neighbours.insert(10); neighbours.insert(11);
822
823
      break;
    case 9:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
824
825
826
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(6); 
      neighbours.insert(8); neighbours.insert(12);
      neighbours.insert(10); neighbours.insert(11); neighbours.insert(14); 
827
828
      break;
    case 10:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
829
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(11);
830
831
      break;
    case 11:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
832
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12);  neighbours.insert(10); neighbours.insert(14);
833
834
      break;
    case 12:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
835
836
837
      neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); 
      neighbours.insert(9); neighbours.insert(13);
      neighbours.insert(11); neighbours.insert(14); neighbours.insert(15); 
838
839
      break;
    case 13:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
840
      neighbours.insert(6); neighbours.insert(7); neighbours.insert(12);  neighbours.insert(14); neighbours.insert(15);
841
842
      break;
    case 14:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
843
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13);  neighbours.insert(11); neighbours.insert(15);
844
845
      break;
    case 15:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
846
      neighbours.insert(12); neighbours.insert(13); neighbours.insert(14);
847
848
      break;
    }    
849

850
    TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
851
852
853
854
855
856
857

    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
858
      case 0: case 1: case 2: case 3:
859
860
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
861
	  m.insert(0); m.insert(1); m.insert(2); m.insert(3);
862
	  levelData.addLevel(m, 0);
863
864
	}
	break;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
865
      case 4: case 5: case 6: case 7:
866
867
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
868
	  m.insert(4); m.insert(5); m.insert(6); m.insert(7);
869
	  levelData.addLevel(m, 1);
870
871
	}
	break;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
872
      case 8: case 9: case 10: case 11:
873
874
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
875
	  m.insert(8); m.insert(9); m.insert(10); m.insert(11);
876
	  levelData.addLevel(m, 2);