ParallelDomainBase.cc 66.5 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
#include <algorithm>
2
3
#include <iostream>
#include <fstream>
Thomas Witkowski's avatar
Thomas Witkowski committed
4
#include <boost/lexical_cast.hpp>
5
#include "parallel/ParallelDomainBase.h"
6
#include "parallel/ParallelDomainDbg.h"
7
#include "parallel/StdMpi.h"
8
9
10
11
12
13
14
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
15
16
#include "DOFMatrix.h"
#include "DOFVector.h"
17
#include "SystemVector.h"
18
#include "VtkWriter.h"
19
#include "ElementDofIterator.h"
20
21
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
22
#include "ElementFileWriter.h"
23
#include "VertexVector.h"
24
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
25
26
#include "ProblemVec.h"
#include "ProblemInstat.h"
27
#include "Debug.h"
28

29
30
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
31
32
  using boost::lexical_cast;

33
34
35
36
37
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

38

39
40
41
42
43
44
45
46
  MeshDistributor::MeshDistributor(std::string str)
    : probStat(0),
      name(str),
      feSpace(NULL),
      mesh(NULL),
      refineManager(NULL),
      info(10),
      partitioner(NULL),
47
      initialPartitionMesh(true),
48
      nRankDofs(0),
49
      nOverallDofs(0),
50
      rstart(0),
51
52
      deserialized(false),
      lastMeshChangeIndex(0)
53
  {
54
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
55

56
57
58
59
60
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
61

62
  void MeshDistributor::initParallelization(AdaptInfo *adaptInfo)
63
  {
64
    FUNCNAME("MeshDistributor::initParallelization()");
65
66
67
68

    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");

69
70
    TEST_EXIT(mesh)("No mesh has been defined for mesh distribution!\n");

71
72
    // If the problem has been already read from a file, we do not need to do anything.
    if (deserialized)
73
      return;
74
    
75

76
77
78
79
80
    // Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported
    // only for macro meshes, so it will not work yet if the mesh is already refined
    // in some way.
    testForMacroMesh();

81

82
83
84
85
86
87
88
    // create an initial partitioning of the mesh
    partitioner->createPartitionData();
    // set the element weights, which are 1 at the very first begin
    setElemWeights(adaptInfo);
    // and now partition the mesh
    partitionMesh(adaptInfo);   

89

Thomas Witkowski's avatar
Thomas Witkowski committed
90
#if (DEBUG != 0)
91
92
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
93
94
95
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
96

97
98
99
100
101
      if (writePartMesh > 0)
	writePartitioningMesh("part.vtu");
      else 
	MSG("Skip write part mesh!\n");
    }
102
103

    ParallelDomainDbg::testAllElements(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
104
#endif
105

106

107
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
108

109
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
110

111
112
113
114
#if (DEBUG != 0)
    ParallelDomainDbg::printBoundaryInfo(*this);
#endif

115

116
117
118
119
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
120
121
122
123
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

124
125
126
127
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

128
129
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
130
#if (DEBUG != 0)
131
    MSG("AMDiS runs in debug mode, so make some test ...\n");
132
133
134
    debug::testSortedDofs(mesh, elMap);
    ParallelDomainDbg::testInteriorBoundary(*this);
    ParallelDomainDbg::testCommonDofs(*this, true);
135
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
136

137
    debug::writeMesh(feSpace, -1, "macro_mesh");   
138
#endif
139

140
141
142
    // === Create periodic dof mapping, if there are periodic boundaries. ===

    createPeriodicMap();
143

144
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
145

Thomas Witkowski's avatar
Thomas Witkowski committed
146
    int globalRefinement = 0;
147
    GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinement);
Thomas Witkowski's avatar
Thomas Witkowski committed
148

Thomas Witkowski's avatar
Thomas Witkowski committed
149
    if (globalRefinement > 0) {
150
      refineManager->globalRefine(mesh, globalRefinement);
151

152
#if (DEBUG != 0)
153
      debug::writeMesh(feSpace, -1, "gr_mesh");
154
155
#endif

156
      updateLocalGlobalNumbering();
157
      
158
      // === Update periodic mapping, if there are periodic boundaries. ===
159
      
160
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
161
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
162
163
164

    /// === Set DOF rank information to all matrices and vectors. ===

165
166
167
168
169
170
    for (unsigned int i = 0; i < probStat.size(); i++) {
      int nComponents = probStat[i]->getNumComponents();
      for (int j = 0; j < nComponents; j++) {
	for (int k = 0; k < nComponents; k++)
	  if (probStat[i]->getSystemMatrix(j, k))
	    probStat[i]->getSystemMatrix(j, k)->setRankDofs(isRankDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
171

172
173
174
175
176
177
	TEST_EXIT_DBG(probStat[i]->getRhs()->getDOFVector(j))("No RHS vector!\n");
	TEST_EXIT_DBG(probStat[i]->getSolution()->getDOFVector(j))("No solution vector!\n");
	
	probStat[i]->getRhs()->getDOFVector(j)->setRankDofs(isRankDof);
	probStat[i]->getSolution()->getDOFVector(j)->setRankDofs(isRankDof);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
180
181
182
    }

    // === Remove periodic boundary conditions in sequential problem definition. ===

    // Remove periodic boundaries in boundary manager on matrices and vectors.
183
184
185
    for (unsigned int i = 0; i < probStat.size(); i++) {
      int nComponents = probStat[i]->getNumComponents();

Thomas Witkowski's avatar
Thomas Witkowski committed
186
      for (int j = 0; j < nComponents; j++) {
187
188
189
190
191
192
193
194
195
	for (int k = 0; k < nComponents; k++) {
	  DOFMatrix* mat = probStat[i]->getSystemMatrix(j, k);
	  if (mat && mat->getBoundaryManager())
	    removeBoundaryCondition(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
	}
	
	if (probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager())
	  removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
	
196
	if (probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager())
197
	  removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
Thomas Witkowski's avatar
Thomas Witkowski committed
198
199
200
201
202
203
204
205
206
207
      }
    }

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh,  -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
208
209
  }

210

211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

    if (feSpace != NULL) {
      std::vector<FiniteElemSpace*> feSpaces = probVec->getFeSpaces();
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	TEST_EXIT(feSpace == feSpaces[i])
	  ("Parallelizaton is not supported for multiple FE spaces!\n");
      }
    } else {
      feSpace = probVec->getFeSpace(0);
      mesh = feSpace->getMesh();
      info = probVec->getInfo();
      
      TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
	("Only meshes with one DOFAdmin are supported!\n");
      TEST_EXIT(mesh->getDOFAdmin(0).getNumberOfPreDOFs(0) == 0)
	("Wrong pre dof number for DOFAdmin!\n");
      
      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());
      }

      partitioner = new ParMetisPartitioner(mesh, &mpiComm);
    }

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", &writeSerialization);
    if (writeSerialization)
      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this));

    int readSerialization = 0;
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", &readSerialization);
    if (readSerialization) {
      ERROR_EXIT("Must be reimplemented!\n");
#if 0      
      std::string filename = "";
      GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
      filename += ".p" + lexical_cast<std::string>(mpiRank);
      MSG("Start serialization with %s\n", filename.c_str());
      std::ifstream in(filename.c_str());
      deserialize(in);
      in.close();
#endif
    }

    probStat.push_back(probVec);
  }


  void MeshDistributor::exitParallelization(AdaptInfo *adaptInfo)
271
  {}
272

273
  
274
  void MeshDistributor::updateDofAdmins()
275
  {
276
    FUNCNAME("MeshDistributor::updateDofAdmins()");
277
278

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
279
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
280
281
282

      // There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise,
      // it is not possible to define a first DOF hole.
283
      if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDofs.size())
284
	admin.enlargeDOFLists();
285
286
287
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
288
      for (unsigned int j = 0; j < mapLocalGlobalDofs.size(); j++)
289
290
 	admin.setDOFFree(j, false);

291
292
293
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
294
295
296
    }
  }

297

298
  void MeshDistributor::testForMacroMesh()
299
  {
300
    FUNCNAME("MeshDistributor::testForMacroMesh()");
301
302
303
304
305
306
307

    int nMacroElements = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      TEST_EXIT(elInfo->getLevel() == 0)
308
	("Mesh is already refined! This does not work with parallelization!\n");
309
310
311
312
313
314
315
316
317
318
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

319

320
  void MeshDistributor::synchVector(DOFVector<double> &vec)
321
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
322
    StdMpi<std::vector<double> > stdMpi(mpiComm);
323
324

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
325
	 sendIt != sendDofs.end(); ++sendIt) {
326
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
329
      
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
332
333
334
335

      stdMpi.send(sendIt->first, dofs);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
338
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
339

Thomas Witkowski's avatar
Thomas Witkowski committed
340
    stdMpi.startCommunication<double>(MPI_DOUBLE);
341

Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
344
345
346
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      for (unsigned int i = 0; i < recvIt->second.size(); i++)
	vec[*(recvIt->second)[i]] = stdMpi.getRecvData(recvIt->first)[i];
  }
347
348


349
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
350
  {
351
    int nComponents = vec.getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
354
355
356
357
358
359
360
361
362
363
    StdMpi<std::vector<double> > stdMpi(mpiComm);

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt) {
      std::vector<double> dofs;
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nComponents * nSendDofs);
      
      for (int i = 0; i < nComponents; i++) {
	DOFVector<double> *dofvec = vec.getDOFVector(i);
	for (int j = 0; j < nSendDofs; j++)
	  dofs.push_back((*dofvec)[*((sendIt->second)[j])]);
364
365
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
366
      stdMpi.send(sendIt->first, dofs);
367
368
369
    }

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size() * nComponents);
372

Thomas Witkowski's avatar
Thomas Witkowski committed
373
    stdMpi.startCommunication<double>(MPI_DOUBLE);
374
375

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
376
377
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
378
379

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
380
381
382
383
384
      for (int i = 0; i < nComponents; i++) {
	DOFVector<double> *dofvec = vec.getDOFVector(i);
 	for (int j = 0; j < nRecvDofs; j++)
	  (*dofvec)[*(recvIt->second)[j]] = 
	    stdMpi.getRecvData(recvIt->first)[counter++];
385
386
387
388
      }
    }
  }

389

390
  void MeshDistributor::removeBoundaryCondition(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
393
394
395
396
397
398
399
400
401
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


402
  void MeshDistributor::checkMeshChange()
403
  {
404
    FUNCNAME("MeshDistributor::checkMeshChange()");
405

406
407
408
409
410
    // === If mesh has not been changed on all ranks, return. ===

    int recvAllValues = 0;
    int sendValue = static_cast<int>(mesh->getChangeIndex() != lastMeshChangeIndex);
    mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
411

412
    if (recvAllValues == 0)
413
414
      return;

415
416
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
417

418
419
420
421
422
    clock_t first = clock();
    
    do {
      // To check the interior boundaries, the ownership of the boundaries is not 
      // important. Therefore, we add all boundaries to one boundary container.
423
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
424
425

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
426
	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
427
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
428
429

      for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
430
	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
431
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
432

433

434
      // === Check the boundaries and adapt mesh if necessary. ===
435

436
437
438
439
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

440
441
442
443
444
445
446
      bool meshChanged = checkAndAdaptBoundary(allBound);

      // === Check on all ranks if at least one rank's mesh has changed. ===

      int sendValue = static_cast<int>(!meshChanged);
      recvAllValues = 0;
      mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
447
448
449
450

#if (DEBUG != 0)
      MSG("Mesh changed on %d ranks!\n", recvAllValues);
#endif
451
    } while (recvAllValues != 0);
452

453
454

#if (DEBUG != 0)
455
    debug::writeMesh(feSpace, -1, "mesh");
456
457
458
459
460
461
#endif

    INFO(info, 8)("Parallel mesh adaption needed %.5f seconds\n", 
		  TIME_USED(first, clock()));

    // === Because the mesh has been changed, update the DOF numbering and mappings. ===
462

463
464
465
466
    updateLocalGlobalNumbering();
  }

  
467
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
468
  {
469
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
470
471
472
473
474
475

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
476
477
478
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
479
	elCode.init(boundIt->rankObj);
480
481
482
483
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
484
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
485
    stdMpi.send(sendCodes);
486
487
488
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
489
    // === Compare received mesh structure codes. ===
490
    
491
492
    bool meshFitTogether = true;

493
494
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
495
496
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
497
      
498
499
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
500

501
502
503
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
504
505
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
506
507

	  bool b = fitElementToMeshCode(recvCodes[i], 
508
					boundIt->rankObj.el,
509
					boundIt->rankObj.subObj,
510
					boundIt->rankObj.ithObj, 
511
					boundIt->rankObj.elType);  
512

513
514
	  if (b)
	    meshFitTogether = false;
515
	}
516
	
517
	i++;
518
519
520
      }
    }

521
    return meshFitTogether;
522
  }
523
524


525
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
526
						Element *el, 
527
528
						GeoIndex subObj,
						int ithObj, 
529
						int elType)
530
  {
531
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");
532

533
534
    TEST_EXIT_DBG(el)("No element given!\n");

535
536
    if (code.empty())
      return false;
537
538
539

    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
540
541
542

    TEST_EXIT_DBG(s1 != -1 || s2 != -1)("This should not happen!\n");

543
544
545
    bool meshChanged = false;

    if (s1 != -1 && s2 != -1) {
546
547
548
549
550
551
552
      TraverseStack stack;
      ElInfo *elInfo = 
	stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
      while (elInfo && elInfo->getElement() != el) {
	elInfo = stack.traverseNext(elInfo);
      }

553
      meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType);
554
555
      return meshChanged;
    }
556
557

    if (el->isLeaf()) {
558
559
560
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

561
562
563
564
565
566
567
568
569
570
571
572
      meshChanged = true;

      TraverseStack stack;
      ElInfo *elInfo = 
	stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);

      while (elInfo && elInfo->getElement() != el) {
	elInfo = stack.traverseNext(elInfo);
      }

      TEST_EXIT_DBG(elInfo)("This should not happen!\n");

573
      el->setMark(1);
574
575
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
576
      refineManager->refineFunction(elInfo);
577
    }
578
579

    if (s1 != -1) {
580
581
582
583
584
585
586
587
      TraverseStack stack;
      ElInfo *elInfo = 
	stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);

      while (elInfo && elInfo->getElement() != el->getFirstChild()) {
	elInfo = stack.traverseNext(elInfo);
      }

588
589
      meshChanged |= 
	fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType));
590
    } else {
591
592
593
594
595
596
597
598
      TraverseStack stack;
      ElInfo *elInfo = 
	stack.traverseFirst(el->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);

      while (elInfo && elInfo->getElement() != el->getSecondChild()) {
	elInfo = stack.traverseNext(elInfo);
      }

599
600
      meshChanged |= 
	fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType));
601
602
603
    }

    return meshChanged;
604
605
  }

606

607
  bool MeshDistributor::fitElementToMeshCode2(MeshStructure &code, 
608
						 TraverseStack &stack,
609
610
						 GeoIndex subObj,
						 int ithObj, 
611
						 int elType)
612
  {
613
    FUNCNAME("MeshDistributor::fitElementToMeshCode2()");
614

615
616
    ElInfo *elInfo = stack.getElInfo();

617
    bool value = false;
618
619
620
621
622
623
624
625
626
627
628
    if (!elInfo)
      return value;

    Element *el = elInfo->getElement();

    if (code.isLeafElement()) {
      int level = elInfo->getLevel();

      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
629
     
630
631
      return value;
    }
632

633
634
    if (!elInfo)
      return value;
635

636
    if (el->isLeaf()) {
637
      el->setMark(1);
638
639
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
640
      refineManager->refineFunction(elInfo);
641
642
643
      value = true;
    }

644
645
    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
646

647
    if (s1 != -1) {
648
      stack.traverseNext(elInfo);
649
      code.nextElement();
650
      value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType));
651
652
653
654
655
      elInfo = stack.getElInfo();
    } else {
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getElement() != el->getSecondChild());      
656
657
    }  

658
659
660
    TEST_EXIT_DBG(elInfo->getElement() == el->getSecondChild())
      ("This should not happen!\n");

661
    if (s2 != -1) {
662
      code.nextElement();
663
      value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType));
664
665
    } else {
      int level = elInfo->getLevel();
666

667
668
669
670
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
671

672
673
674
    return value;
  }

675
  
676
  void MeshDistributor::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
677
  {    
678
    int vecSize = data.size();
679
    SerUtil::serialize(out, vecSize);
680
681
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
682
      SerUtil::serialize(out, dofIndex);
683
684
685
686
    }
  }


687
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
688
689
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
690
    FUNCNAME("MeshDistributor::deserialize()");
691
692

    int vecSize = 0;
693
    SerUtil::deserialize(in, vecSize);
694
695
696
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
697
      SerUtil::deserialize(in, dofIndex);
698
699
700
701
702
703
704
705
706

      TEST_EXIT_DBG(dofMap.count(dofIndex) != 0)
	("Dof index could not be deserialized correctly!\n");

      data[i] = dofMap[dofIndex];
    }
  }


707
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)
708
709
  {
    int mapSize = data.size();
710
    SerUtil::serialize(out, mapSize);
711
712
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
713
      SerUtil::serialize(out, rank);
714
715
716
717
718
      serialize(out, it->second);
    }
  }

  
719
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
720
721
722
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    int mapSize = 0;
723
    SerUtil::deserialize(in, mapSize);
724
725
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
726
      SerUtil::deserialize(in, rank);
727
728
729
730
      deserialize(in, data[rank], dofMap);      
    }
  }

731

732
  double MeshDistributor::setElemWeights(AdaptInfo *adaptInfo) 
733
  {
734
    FUNCNAME("MeshDistributor::setElemWeights()");
735

736
    double localWeightSum = 0.0;
737
    elemWeights.clear();
738
739
740
741
742
743
744
745
746
747
748
749
750
751
      
    std::string filename = "";
    GET_PARAMETER(0, mesh->getName() + "->macro weights", &filename);
    if (filename != "") {
      MSG("Read macro weights from %s\n", filename.c_str());

      std::ifstream infile;
      infile.open(filename.c_str(), std::ifstream::in);
      while (!infile.eof()) {
	int elNum, elWeight;
	infile >> elNum;
	if (infile.eof())
	  break;
	infile >> elWeight;
752

753
754
755
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
756

757
      infile.close();
758
    } else {           
759
      TraverseStack stack;
760
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
761
      while (elInfo) {
762
763
764
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
	localWeightSum++;

765
	elInfo = stack.traverseNext(elInfo);
766
767
768
769
770
771
      }
    }

    return localWeightSum;
  }

772

773
  void MeshDistributor::partitionMesh(AdaptInfo *adaptInfo)
774
  {
775
    FUNCNAME("MeshDistributor::partitionMesh()");
776

777
778
779
780
781
782
783
784
785
786
787
788
    if (initialPartitionMesh) {
      initialPartitionMesh = false;
      partitioner->fillCoarsePartitionVec(&oldPartitionVec);
      partitioner->partition(&elemWeights, INITIAL);
    } else {
      oldPartitionVec = partitionVec;
      partitioner->partition(&elemWeights, ADAPTIVE_REPART, 100.0 /*0.000001*/);
    }    

    partitioner->fillCoarsePartitionVec(&partitionVec);
  }

789

790
  void MeshDistributor::createInteriorBoundaryInfo()
791
  {
792
    FUNCNAME("MeshDistributor::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
793

Thomas Witkowski's avatar
Thomas Witkowski committed
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
    // === First, create the interior boundaries based on macro element's  ===
    // === neighbour informations.                                         ===

    createMacroElementInteriorBoundaryInfo();

    // === Second, search the whole mesh for interior boundaries that consists of ===
    // === substructures of the macro elements.                                   ===

    createSubstructureInteriorBoundaryInfo();


    // === Once we have this information, we must care about the order of the atomic ===
    // === bounds in the three boundary handling object. Eventually all the bound-   ===
    // === aries have to be in the same order on both ranks that share the bounday.  ===

    StdMpi<std::vector<AtomicBoundary> > stdMpi(mpiComm);
    stdMpi.send(myIntBoundary.boundary);
    stdMpi.recv(otherIntBoundary.boundary);
    stdMpi.startCommunication<int>(MPI_INT);
813

Thomas Witkowski's avatar
Thomas Witkowski committed
814
815
816
817
818
819
820
821
    // === The information about all neighbouring boundaries has been received. So ===
    // === the rank tests if its own atomic boundaries are in the same order. If   ===
    // === not, the atomic boundaries are swaped to the correct order.             ===

    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {

      // === We have received from rank "rankIt->first" the ordered list of element ===
822
      // === indices. Now, we have to sort the corresponding list in this rank to   ===
Thomas Witkowski's avatar
Thomas Witkowski committed
823
      // === get the same order.                                                    ===
824
     
Thomas Witkowski's avatar
Thomas Witkowski committed
825
      for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
826

Thomas Witkowski's avatar
Thomas Witkowski committed
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
	// If the expected object is not at place, search for it.

	BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj;

	if ((rankIt->second)[j].neighObj != recvedBound) {
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
 	    if ((rankIt->second)[k].neighObj == recvedBound)
	      break;

	  // The element must always be found, because the list is just in another order.
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
	    ("Should never happen!\n");

	  // Swap the current with the found element.
	  AtomicBoundary tmpBound = (rankIt->second)[k];
	  (rankIt->second)[k] = (rankIt->second)[j];
	  (rankIt->second)[j] = tmpBound;	
	}
      }
    }

    // === Do the same for the periodic boundaries. ===

    if (periodicBoundary.boundary.size() > 0) {
      stdMpi.clear();

      InteriorBoundary sendBounds, recvBounds;     
      for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
	   rankIt != periodicBoundary.boundary.end(); ++rankIt) {

	TEST_EXIT_DBG(rankIt->first != mpiRank)
	  ("It is no possible to have an interior boundary within a rank partition!\n");

	if (rankIt->first < mpiRank)
	  sendBounds.boundary[rankIt->first] = rankIt->second;
	else
	  recvBounds.boundary[rankIt->first] = rankIt->second;
      }

      stdMpi.send(sendBounds.boundary);
      stdMpi.recv(recvBounds.boundary);
      stdMpi.startCommunication<int>(MPI_INT);

      for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
	   rankIt != periodicBoundary.boundary.end(); ++rankIt) {
	if (rankIt->first <= mpiRank)
	  continue;
	  
	for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
	  
	  BoundaryObject &recvedBound = stdMpi.getRecvData()[rankIt->first][j].rankObj;
	  
	  if (periodicBoundary.boundary[rankIt->first][j].neighObj != recvedBound) {    
	    int k = j + 1;	    
	    for (; k < static_cast<int>(rankIt->second.size()); k++)
	      if (periodicBoundary.boundary[rankIt->first][k].neighObj == recvedBound)
		break;
	    
	    // The element must always be found, because the list is just in 
	    // another order.
	    TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
	      ("Should never happen!\n");
	    
	    // Swap the current with the found element.
	    AtomicBoundary tmpBound = (rankIt->second)[k];
	    (rankIt->second)[k] = (rankIt->second)[j];
	    (rankIt->second)[j] = tmpBound;	
	  }  	  
	}
      }     
    } // periodicBoundary.boundary.size() > 0
  }


902
  void MeshDistributor::createMacroElementInteriorBoundaryInfo()
Thomas Witkowski's avatar
Thomas Witkowski committed
903
  {
904
    FUNCNAME("MeshDistributor::createMacroElementInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
905

906
907
    int nNeighbours = mesh->getGeo(NEIGH);
    int dim = mesh->getDim();
908
    GeoIndex subObj = INDEX_OF_DIM(dim - 1, dim);
909
910
911

    // === First, traverse the mesh and search for all elements that have an  ===
    // === boundary with an element within another partition.                 ===
Thomas Witkowski's avatar
Thomas Witkowski committed
912

913
    TraverseStack stack;
914
915
916
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
917
918
919
    while (elInfo) {
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
920
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
921

922
      // Check, if the element is within rank's partition.
923
      if (partitionData->getPartitionStatus() == IN) {
924
	for (int i = 0; i < nNeighbours; i++) {
925
926
927
928
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
929
930
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->
						getElementData(PARTITION_ED));
931

932
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
Thomas Witkowski's avatar
Thomas Witkowski committed
933

934
935
936
937
938
	    // We have found an element that is in rank's partition, but has a 
	    // neighbour outside of the rank's partition.

	    // === Create information about the boundary between the two elements. ===

939
	    AtomicBoundary bound;	    	    
940
941
	    bound.rankObj.el = element;
	    bound.rankObj.elIndex = element->getIndex();
942
	    bound.rankObj.elType = elInfo->getType();
943
	    bound.rankObj.subObj = subObj;
944
	    bound.rankObj.ithObj = i;
945
946

	    bound.neighObj.el = elInfo->getNeighbour(i);
947
	    bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex();
948
	    bound.neighObj.elType = -1;
949
	    bound.neighObj.subObj = subObj;
950
	    bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i);
951
	    
952
	    // Get rank number of the neighbouring element.
953
954
	    int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];

955
956
957
958
959
960
961
962
	    // === Add the boundary information object to the corresponding overall ===
	    // === boundary. There are three possibilities: if the boundary is a    ===
	    // === periodic boundary, just add it to \ref periodicBounadry. Here it ===
	    // === does not matter which rank is responsible for this boundary.     ===
	    // === Otherwise, the boundary is added either to \ref myIntBoundary or ===
	    // === to \ref otherIntBoundary. It dependes on which rank is respon-   ===
	    // === sible for this boundary.                                         ===

963
	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(subObj, i))) {	      
964
965
966
967
968
	      // We have found an element that is at an interior, periodic boundary.
	      AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
	      b = bound;
	    } else {
	      // We have found an element that is at an interior, non-periodic boundary.
969
970
971
972
973
974
975
976
977
	      if (mpiRank > otherElementRank) {
		AtomicBoundary& b = myIntBoundary.getNewAtomic(otherElementRank);
		b = bound;
		b.rankObj.setReverseMode(b.neighObj, feSpace);
	      } else {
		AtomicBoundary& b = otherIntBoundary.getNewAtomic(otherElementRank);
		b = bound;	 
		b.neighObj.setReverseMode(b.rankObj, feSpace);
	      }	      
978
	    }
979
980
981
982
983
984
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
985
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
986
987


988
  void MeshDistributor::createSubstructureInteriorBoundaryInfo()
Thomas Witkowski's avatar
Thomas Witkowski committed
989
  {
990
    FUNCNAME("MeshDistributor::createSubstructureInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
991

992
993
994
    // === Seach for all vertices/edges, which are part of an interior boundary,  ===
    // === but are not part of the interior boundaries that are created based on  ===
    // === the information of macro elements.                                     ===
Thomas Witkowski's avatar
Thomas Witkowski committed
995
996

    int dim = mesh->getDim();
997
    const BasisFunction *basFcts = feSpace->getBasisFcts();
998
    std::vector<DegreeOfFreedom> localIndices(basFcts->getNumber());