ParallelDomainBase.cc 65.2 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);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
38
39
40
41
42
43
44
45
46
47
  ParallelDomainBase::ParallelDomainBase(ProblemVec *problemStat,
					 ProblemInstatVec *problemInstat)
    : iterationIF(problemStat),
      timeIF(problemInstat),
      probStat(problemStat),
      name(problemStat->getName()),
      feSpace(problemStat->getFESpace(0)),
      mesh(feSpace->getMesh()),
      refineManager(problemStat->getRefinementManager(0)),
      info(problemStat->getInfo()),
48
      initialPartitionMesh(true),
49
      nRankDofs(0),
50
      rstart(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
51
      nComponents(problemStat->getNumComponents()),
52
53
      deserialized(false),
      lastMeshChangeIndex(0)
54
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
55
56
57
58
59
60
61
    FUNCNAME("ParallelDomainBase::ParalleDomainBase()");

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

62
63
64
65
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
    partitioner = new ParMetisPartitioner(mesh, &mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90

    // Test if all fe spaces are equal. Yet, different fe spaces are not supported for
    // domain parallelization.
    const FiniteElemSpace *fe = probStat->getFESpace(0);
    for (int i = 0; i < nComponents; i++)
      TEST_EXIT(fe == probStat->getFESpace(i))
	("Parallelization does not supported different FE spaces!\n");

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

    int readSerialization = 0;
    GET_PARAMETER(0, name + "->input->read serialization", "%d", &readSerialization);
    if (readSerialization) {
      std::string filename = "";
      GET_PARAMETER(0, name + "->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();
    }
91
92
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
93

94
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
95
  {
96
97
98
99
100
101
102
    FUNCNAME("ParallelDomainBase::initParallelization()");

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

    // If the problem has been already read from a file, we do not need to do anything.
    if (deserialized)
103
104
      return;

105
106
107
108
109
    // 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();

110
111
    MSG("START PARTITIONING!\n");

112
113
114
115
116
117
118
    // 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);   

Thomas Witkowski's avatar
Thomas Witkowski committed
119
#if (DEBUG != 0)
120
121
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
122
123
124
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
125

126
127
128
129
130
      if (writePartMesh > 0)
	writePartitioningMesh("part.vtu");
      else 
	MSG("Skip write part mesh!\n");
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
131
#endif
132

133
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
134

135
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
136

137
138
139
140
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
141
142
143
144
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();

145
146
147
148
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

149
150
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
151
#if (DEBUG != 0)
152
    MSG("AMDiS runs in debug mode, so make some test ...\n");
153
154
155
    debug::testSortedDofs(mesh, elMap);
    ParallelDomainDbg::testInteriorBoundary(*this);
    ParallelDomainDbg::testCommonDofs(*this, true);
156
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
157

158
159
    debug::writeMesh(feSpace, -1, "macromesh");   
#endif
160

161
162
163
    // === Create periodic dof mapping, if there are periodic boundaries. ===

    createPeriodicMap();
164

165
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
166

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

Thomas Witkowski's avatar
Thomas Witkowski committed
170
    if (globalRefinement > 0) {
171
      refineManager->globalRefine(mesh, globalRefinement);
172

173
      updateLocalGlobalNumbering();
174

175
      // === Update periodic mapping, if there are periodic boundaries. ===
176

177
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
178
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218


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

    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++)
 	if (probStat->getSystemMatrix(i, j))
 	  probStat->getSystemMatrix(i, j)->setRankDofs(isRankDof);      

      TEST_EXIT_DBG(probStat->getRHS()->getDOFVector(i))("No rhs vector!\n");
      TEST_EXIT_DBG(probStat->getSolution()->getDOFVector(i))("No solution vector!\n");

      probStat->getRHS()->getDOFVector(i)->setRankDofs(isRankDof);
      probStat->getSolution()->getDOFVector(i)->setRankDofs(isRankDof);
    }

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

    // Remove periodic boundaries in boundary manager on matrices and vectors.
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	DOFMatrix* mat = probStat->getSystemMatrix(i, j);
 	if (mat && mat->getBoundaryManager())
	  removeBoundaryCondition(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
      }

      if (probStat->getSolution()->getDOFVector(i)->getBoundaryManager())
	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat->getSolution()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));

      if (probStat->getRHS()->getDOFVector(i)->getBoundaryManager())
	removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probStat->getRHS()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap()));
    }

    // 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);
    }    
219
220
  }

221
222

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
223
  {}
224

225
226
227
  
  void ParallelDomainBase::updateDofAdmins()
  {
228
229
230
    FUNCNAME("ParallelDomainBase::updateDofAdmins()");

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
231
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
232
233
234

      // There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise,
      // it is not possible to define a first DOF hole.
235
      if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDofs.size())
236
	admin.enlargeDOFLists();
237
238
239
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
240
      for (unsigned int j = 0; j < mapLocalGlobalDofs.size(); j++)
241
242
 	admin.setDOFFree(j, false);

243
244
245
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
246
247
248
    }
  }

249

250
251
252
253
254
255
256
257
258
259
  void ParallelDomainBase::testForMacroMesh()
  {
    FUNCNAME("ParallelDomainBase::testForMacroMesh()");

    int nMacroElements = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      TEST_EXIT(elInfo->getLevel() == 0)
260
	("Mesh is already refined! This does not work with parallelization!\n");
261
262
263
264
265
266
267
268
269
270
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
  void ParallelDomainBase::synchVector(DOFVector<double> &vec)
273
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
274
    StdMpi<std::vector<double> > stdMpi(mpiComm);
275
276

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
277
	 sendIt != sendDofs.end(); ++sendIt) {
278
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
279
280
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
281
      
Thomas Witkowski's avatar
Thomas Witkowski committed
282
283
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
284
285
286
287

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

Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
290
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
291

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

Thomas Witkowski's avatar
Thomas Witkowski committed
294
295
296
297
298
    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];
  }
299
300


Thomas Witkowski's avatar
Thomas Witkowski committed
301
302
303
304
305
306
307
308
309
310
311
312
313
314
  void ParallelDomainBase::synchVector(SystemVector &vec)
  {
    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])]);
315
316
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
317
      stdMpi.send(sendIt->first, dofs);
318
319
320
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
324
    stdMpi.startCommunication<double>(MPI_DOUBLE);
325
326

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
329
330

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
331
332
333
334
335
      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++];
336
337
338
339
      }
    }
  }

340

Thomas Witkowski's avatar
Thomas Witkowski committed
341
342
343
344
345
346
347
348
349
350
351
352
  void ParallelDomainBase::removeBoundaryCondition(BoundaryIndexMap& boundaryMap)
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


353
354
  void ParallelDomainBase::checkMeshChange()
  {
355
356
    FUNCNAME("ParallelDomainBase::checkMeshChange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
357
358
    int dim = mesh->getDim();

359
360
361
362
363
    // === 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);
364

365
    if (recvAllValues == 0)
366
367
      return;

368
369
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
370

371
372
373
374
375
    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.
376
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
379
380
381
382
383
384
385

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
	if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
	  allBound[it.getRank()].push_back(it.getBound());

      for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
	if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
	  allBound[it.getRank()].push_back(it.getBound());

386

387
      // === Check the boundaries and adapt mesh if necessary. ===
388

389
390
391
392
393
394
395
396
      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);
    } while (recvAllValues != 0);
397

398
399

#if (DEBUG != 0)
400
    debug::writeMesh(feSpace, -1, "mesh");
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
#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. ===
   
    updateLocalGlobalNumbering();
  }

  
  bool ParallelDomainBase::checkAndAdaptBoundary(RankToBoundMap &allBound)
  {
    FUNCNAME("ParallelDomainBase::checkAndAdaptBoundary()");

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
421
422
423
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
424
	elCode.init(boundIt->rankObj);
425
426
427
428
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
429
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
430
    stdMpi.send(sendCodes);
431
432
433
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
434
    // === Compare received mesh structure codes. ===
435
    
436
437
    bool meshFitTogether = true;

438
439
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
440
441
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
442
      
443
444
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
445

446
447
448
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
449
450
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
451
452
453

	  int refCounter = 0;
	  bool b = fitElementToMeshCode(recvCodes[i], 
454
455
					boundIt->rankObj.el,
					boundIt->rankObj.ithObj, 
456
457
					boundIt->rankObj.elType, refCounter);  

458
459
	  if (b)
	    meshFitTogether = false;
460
	}
461
	
462
	i++;
463
464
465
      }
    }

466
    return meshFitTogether;
467
  }
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540


  bool ParallelDomainBase::fitElementToMeshCode(MeshStructure &code, 
						Element *el, 
						int ithSide, 
						int elType,
						int &refCounter)
  {
    FUNCNAME("fitElementToMeshCode()");

    if (code.empty())
      return false;
    
    int s1 = el->getSideOfChild(0, ithSide, elType);
    int s2 = el->getSideOfChild(1, ithSide, elType);

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

    if (s1 != -1 && s2 != -1)
      return fitElementToMeshCode2(code, el, ithSide, elType, refCounter);

    if (el->isLeaf()) {
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;
      
      el->setMark(1);
      refineManager->refineElement(el->getMesh(), el);
      refCounter++;
    }
   
    if (s1 != -1)
      return fitElementToMeshCode2(code, 
				   el->getFirstChild(), s1, el->getChildType(elType), refCounter);
    else
      return fitElementToMeshCode2(code, 
				   el->getSecondChild(), s2, el->getChildType(elType), refCounter);
  }

  bool ParallelDomainBase::fitElementToMeshCode2(MeshStructure &code, 
						 Element *el, 
						 int ithSide, 
						 int elType, int &refCounter)
  {
    if (code.isLeafElement())
      return false;

    bool value = false;

    if (el->isLeaf()) {
      el->setMark(1);
      refineManager->refineElement(el->getMesh(), el);     
      value = true;
      refCounter++;
    }

    int s1 = el->getSideOfChild(0, ithSide, elType);
    int s2 = el->getSideOfChild(1, ithSide, elType);
    
    if (s1 != -1) {
      code.nextElement();
      value |= fitElementToMeshCode2(code, el->getFirstChild(), 
				     s1, el->getChildType(elType), refCounter);
    }  

    if (s2 != -1) {
      code.nextElement();	
      value |= fitElementToMeshCode2(code, el->getSecondChild(), 
				     s2, el->getChildType(elType), refCounter);
    }

    return value;
  }

541
542
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
543
  {    
544
    int vecSize = data.size();
545
    SerUtil::serialize(out, vecSize);
546
547
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
548
      SerUtil::serialize(out, dofIndex);
549
550
551
552
553
554
555
556
557
558
    }
  }


  void ParallelDomainBase::deserialize(std::istream &in, DofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    FUNCNAME("ParallelDomainBase::deserialize()");

    int vecSize = 0;
559
    SerUtil::deserialize(in, vecSize);
560
561
562
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
563
      SerUtil::deserialize(in, dofIndex);
564
565
566
567
568
569
570
571
572
573
574
575

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

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


  void ParallelDomainBase::serialize(std::ostream &out, RankToDofContainer &data)
  {
    int mapSize = data.size();
576
    SerUtil::serialize(out, mapSize);
577
578
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
579
      SerUtil::serialize(out, rank);
580
581
582
583
584
585
586
587
588
      serialize(out, it->second);
    }
  }

  
  void ParallelDomainBase::deserialize(std::istream &in, RankToDofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    int mapSize = 0;
589
    SerUtil::deserialize(in, mapSize);
590
591
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
592
      SerUtil::deserialize(in, rank);
593
594
595
596
      deserialize(in, data[rank], dofMap);      
    }
  }

597

598
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
599
  {
600
    FUNCNAME("ParallelDomainBase::setElemWeights()");
601

602
    double localWeightSum = 0.0;
603
    elemWeights.clear();
604
605
606
607
608
609
610
611
612
613
614
615
616
617
      
    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;
618

619
620
621
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
622

623
624
625
626
627
628
629
630
631
632
633
      infile.close();
    } else {
           
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
      while (elInfo) {
	Element *element = elInfo->getElement();
	
	// get partition data
	PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	  (element->getElementData(PARTITION_ED));
634
	
635
636
637
638
639
640
641
642
643
644
	if (partitionData && partitionData->getPartitionStatus() == IN) {
	  int elNum = -1;
	  if (partitionData->getLevel() == 0)
	    elNum = element->getIndex();
	  
	  TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
	  if (element->isLeaf()) {
	    elemWeights[elNum] += 1.0;
	    localWeightSum += 1.0;
	  }
645
	}
646
647
	
	elInfo = stack.traverseNext(elInfo);
648
649
650
651
652
653
      }
    }

    return localWeightSum;
  }

654

655
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
656
  {
657
658
    FUNCNAME("ParallelDomainBase::partitionMesh()");

659
660
661
662
663
664
665
666
667
668
669
670
    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);
  }

671

Thomas Witkowski's avatar
Thomas Witkowski committed
672
673
674
675
676
677
678
679
680
  void ParallelDomainBase::solveInitialProblem(AdaptInfo *adaptInfo)
  {     
    if (timeIF)
      timeIF->solveInitialProblem(adaptInfo);
    
    synchVector(*(probStat->getSolution()));
  }


681
  void ParallelDomainBase::createInteriorBoundaryInfo()
682
  {
683
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
684

Thomas Witkowski's avatar
Thomas Witkowski committed
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
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
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
    // === 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);
   
    // === 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 ===
      // === indices. We now have to sort the corresponding list in this rank to    ===
      // === get the same order.                                                    ===
      
      for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
	// 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
  }


  void ParallelDomainBase::createMacroElementInteriorBoundaryInfo()
  {
    FUNCNAME("ParallelDomainBase::createMacroElementInteriorBoundaryInfo()");

796
797
798
799
800
801
802
803
804
805
806
807
808
    int nNeighbours = mesh->getGeo(NEIGH);
    int dim = mesh->getDim();
    GeoIndex subObj = CENTER;
    switch (dim) {
    case 2:
      subObj = EDGE;
      break;
    case 3:
      subObj = FACE;
      break;
    default:
      ERROR_EXIT("What is this?\n");
    }     
809
810
811

    // === 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
812

813
    TraverseStack stack;
814
815
816
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
817
818
819
820
    while (elInfo) {
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
821

822
      // Check, if the element is within rank's partition.
823
      if (partitionData->getPartitionStatus() == IN) {
824
	for (int i = 0; i < nNeighbours; i++) {
825
826
827
828
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
829
830
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->
						getElementData(PARTITION_ED));
831

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

834
835
836
837
838
	    // 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. ===

839
	    AtomicBoundary bound;	    	    
840
841
	    bound.rankObj.el = element;
	    bound.rankObj.elIndex = element->getIndex();
842
	    bound.rankObj.elType = elInfo->getType();
843
	    bound.rankObj.subObj = subObj;
844
	    bound.rankObj.ithObj = i;
845
846

	    bound.neighObj.el = elInfo->getNeighbour(i);
847
	    bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex();
848
	    bound.neighObj.elType = -1;
849
	    bound.neighObj.subObj = subObj;
850
	    bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i);
851
	    
852
853
854
855
856
	    if (dim == 2) {
	      // i == 2  =>  getSideOfNeighbour(i) == 2
	      TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
		("Should not happen!\n");
	    }
857

858
	    // Get rank number of the neighbouring element.
859
860
	    int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];

861
862
863
864
865
866
867
868
869

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

870
	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(subObj, i))) {	      
871
872
873
874
875
	      // 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.
876
877
878
879
880
881
882
883
884
	      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);
	      }	      
885
	    }
886
887
888
889
890
891
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
892
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
893
894


Thomas Witkowski's avatar
Thomas Witkowski committed
895
896
897
  void ParallelDomainBase::createSubstructureInteriorBoundaryInfo()
  {
    FUNCNAME("ParallelDomainBase::createSubstructureInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
898

Thomas Witkowski's avatar
Thomas Witkowski committed
899
900
901
902
903
    // === Seach for all vertices/edges, which are part of an interior boundary,   ===
    // === but are not a part of the interior boundaries that are created based on ===
    // === the information of macro elements.                                      ===

    int dim = mesh->getDim();
904
    const BasisFunction *basFcts = feSpace->getBasisFcts();
Thomas Witkowski's avatar
Thomas Witkowski committed
905
906
907
908
909
910
911
    std::vector<DegreeOfFreedom> localIndices(basFcts->getNumber(), 0);

    // Maps each DOF in the whole domain to the rank that ownes this DOF.
    std::map<DegreeOfFreedom, int> dofOwner;
    // Maps each DOF in ranks partition of the domain to the element object that 
    // contains this DOF.
    std::map<DegreeOfFreedom, BoundaryObject> rankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
912

913
914
    // Maps each edge in the whole domain to the rank that ownes this edge.
    std::map<GlobalEdge, int> edgeOwner;
Thomas Witkowski's avatar
Thomas Witkowski committed
915
916
    // Maps each edge in ranks partition of the domain to the element object that 
    // contains this edge.
917
    std::map<GlobalEdge, BoundaryObject> rankEdges;
918

919

Thomas Witkowski's avatar
Thomas Witkowski committed
920
    // === Traverse whole mesh and fill the maps defined above.                    ===
921

Thomas Witkowski's avatar
Thomas Witkowski committed
922
923
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
924
925
926
927
928
929
930
    while (elInfo) {
      Element *el = elInfo->getElement();
      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);

      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(el->getElementData(PARTITION_ED));
      
Thomas Witkowski's avatar
Thomas Witkowski committed
931
932
933
934
935
936
937
938
939
940
941
942
943
944
      // In 3d, traverse all edges of current element.
      if (dim == 3) {
	for (int i = 0; i < 6; i++) {
	  DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(i, 0)];
	  DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(i, 1)];
	  GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));

	  // Update the owner of the current edge.
	  edgeOwner[edge] = max(edgeOwner[edge], partitionVec[el->getIndex()]);
	  
	  // If the edge is part of an element that is part of rank's domain, add it
	  // to the set of all rank's edges.
	  if (partitionData && partitionData->getPartitionStatus() == IN)
	    rankEdges[edge] = BoundaryObject(el, elInfo->getType(), EDGE, i);
945
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
946
947
948
949
950
951
952
953
954
955
      }

      
      // In 2d and 3d, traverse all vertices of current element.
      for (int i = 0; i < 4; i++) {
	DegreeOfFreedom dof0 = localIndices[i];
	dofOwner[dof0] = max(dofOwner[dof0], partitionVec[el->getIndex()]);

	if (partitionData && partitionData->getPartitionStatus() == IN)
	  rankDofs[dof0] = BoundaryObject(el, elInfo->getType(), VERTEX, i);
956
957
958
959
960
      }      

      elInfo = stack.traverseNext(elInfo);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
961
962

    // === Create a set of all edges and vertices at rank's interior boundaries. ===
963
964
965

    // Stores all edges at rank's interior boundaries.
    std::set<GlobalEdge> rankBoundaryEdges;
Thomas Witkowski's avatar
Thomas Witkowski committed
966
    // Stores all vertices at rank's interior boundaries.
967
968
    std::set<DegreeOfFreedom> rankBoundaryDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
969
970
   
    // First, traverse the rank owned elements at the interior boundaries.
Thomas Witkowski's avatar
Thomas Witkowski committed
971
972
973
    for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
      Element *el = it.getBound().rankObj.el;
      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
974
      
Thomas Witkowski's avatar
Thomas Witkowski committed
975
      if (dim == 3) {
976
	for (int j = 0; j < 3; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
977
978
979
980
981
982
983
984
985
986
987
	  int edgeNo = el->getEdgeOfFace(it.getBound().rankObj.ithObj, j);
	  DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)];
	  DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)];
	  GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1));	
	  
	  // If the edge at rank's interior boundarie has a higher owner rank, than
	  // we have to remove this edge from the corresponding boundary element.
	  // Otherwise, it is part of the interior boundary and we add it to the set
	  // rankBoundaryEdges.
	  if (edgeOwner[edge] > mpiRank)
	    it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
988
	  else
Thomas Witkowski's avatar
Thomas Witkowski committed
989
	    rankBoundaryEdges.insert(edge);
990
	}
991
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
992
993
994
995
996
997
998
999
1000
      
      for (int j = 0; j < 3; j++) {
	int dofNo = el->getVertexOfPosition(FACE, it.getBound().rankObj.ithObj, j);
	DegreeOfFreedom dof = localIndices[dofNo];
	
	if (dofOwner[dof] > mpiRank)
	  it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
	else
	  rankBoundaryDofs.insert(dof);
For faster browsing, not all history is shown. View entire blame