ParallelDomainBase.cc 76.8 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
6
#include "parallel/ParallelDomainBase.h"
#include "parallel/StdMpi.h"
7
8
9
10
11
12
13
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
14
15
#include "DOFMatrix.h"
#include "DOFVector.h"
16
#include "SystemVector.h"
17
#include "VtkWriter.h"
18
#include "ElementDofIterator.h"
19
20
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
21
#include "ElementFileWriter.h"
22
#include "VertexVector.h"
23
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24
25
#include "ProblemVec.h"
#include "ProblemInstat.h"
26
#include "Debug.h"
27

28
29
namespace AMDiS {

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
37
38
39
40
41
42
43
44
45
46
  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()),
47
      initialPartitionMesh(true),
48
      nRankDofs(0),
49
      rstart(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
50
      nComponents(problemStat->getNumComponents()),
51
52
      deserialized(false),
      lastMeshChangeIndex(0)
53
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
54
55
56
57
58
59
60
    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");

61
62
63
64
    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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

    // 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();
    }
90
91
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
92

93
  void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo)
94
  {
95
96
97
98
99
100
101
    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)
102
103
      return;

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

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

111
112
113
114
115
116
117
    // 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
118
119
#if (DEBUG != 0)
    ElementIdxToDofs elMap;
120
    dbgCreateElementMap(elMap);
121
122
123
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
124

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

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

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

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

    createLocalGlobalNumbering();

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

    removeMacroElements();

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

    updateDofAdmins();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
150
#if (DEBUG != 0)
151
    MSG("AMDiS runs in debug mode, so make some test ...\n");
152
153
    dbgTestElementMap(elMap);
    dbgTestInteriorBoundary();
154
    dbgTestCommonDofs(true);
155
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
156

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

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

    createPeriodicMap();
163

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

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

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

172
      updateLocalGlobalNumbering();
173

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

176
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
177
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
178
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


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

220
221

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

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

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

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

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

248

249
250
251
252
253
254
255
256
257
258
  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)
259
	("Mesh is already refined! This does not work with parallelization!\n");
260
261
262
263
264
265
266
267
268
269
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

270

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

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

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

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

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

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


Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
302
303
304
305
306
307
308
309
310
311
312
313
  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])]);
314
315
      }

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

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

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

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

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

339

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


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

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

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

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

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

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

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

385

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

388
389
390
391
392
393
394
395
      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);
396

397
398

#if (DEBUG != 0)
399
    debug::writeMesh(feSpace, -1, "mesh");
400
401
402
403
404
405
406
407
#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();
408
409
410
411

    // === If there is a non zero pattern computed for Petsc matrix, delete it. So ===
    // === it will be recomputed after next assembling.                            ===

412
#if 0
413
414
415
416
417
418
419
420
421
    if (d_nnz) {
      delete [] d_nnz;
      d_nnz = NULL;
    }

    if (o_nnz) {
      delete [] o_nnz;
      o_nnz = NULL;
    }
422
#endif
423
424
425
426
427
428
429
430
431
432
433
434
  }

  
  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) {
435
436
437
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
438
	elCode.init(boundIt->rankObj);
439
440
441
442
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
443
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
444
    stdMpi.send(sendCodes);
445
446
447
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
448
    // === Compare received mesh structure codes. ===
449
    
450
451
    bool meshFitTogether = true;

452
453
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
454
455
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
456
      
457
458
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
459

460
461
462
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
463
464
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
465
466
467

	  int refCounter = 0;
	  bool b = fitElementToMeshCode(recvCodes[i], 
468
469
					boundIt->rankObj.el,
					boundIt->rankObj.ithObj, 
470
471
					boundIt->rankObj.elType, refCounter);  

472
473
	  if (b)
	    meshFitTogether = false;
474
	}
475
	
476
	i++;
477
478
479
      }
    }

480
    return meshFitTogether;
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
541
542
543
544
545
546
547
548
549
550
551
552
553
554


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

555
556
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
557
  {    
558
    int vecSize = data.size();
559
    SerUtil::serialize(out, vecSize);
560
561
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
562
      SerUtil::serialize(out, dofIndex);
563
564
565
566
567
568
569
570
571
572
    }
  }


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

    int vecSize = 0;
573
    SerUtil::deserialize(in, vecSize);
574
575
576
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
577
      SerUtil::deserialize(in, dofIndex);
578
579
580
581
582
583
584
585
586
587
588
589

      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();
590
    SerUtil::serialize(out, mapSize);
591
592
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
593
      SerUtil::serialize(out, rank);
594
595
596
597
598
599
600
601
602
      serialize(out, it->second);
    }
  }

  
  void ParallelDomainBase::deserialize(std::istream &in, RankToDofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    int mapSize = 0;
603
    SerUtil::deserialize(in, mapSize);
604
605
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
606
      SerUtil::deserialize(in, rank);
607
608
609
610
      deserialize(in, data[rank], dofMap);      
    }
  }

611

612
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
613
  {
614
    FUNCNAME("ParallelDomainBase::setElemWeights()");
615

616
    double localWeightSum = 0.0;
617
    elemWeights.clear();
618
619
620
621
622
623
624
625
626
627
628
629
630
631
      
    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;
632

633
634
635
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
636

637
638
639
640
641
642
643
644
645
646
647
      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));
648
	
649
650
651
652
653
654
655
656
657
658
	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;
	  }
659
	}
660
661
	
	elInfo = stack.traverseNext(elInfo);
662
663
664
665
666
667
      }
    }

    return localWeightSum;
  }

668

669
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
670
  {
671
672
    FUNCNAME("ParallelDomainBase::partitionMesh()");

673
674
675
676
677
678
679
680
681
682
683
684
    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);
  }

685

Thomas Witkowski's avatar
Thomas Witkowski committed
686
687
688
689
690
691
692
693
694
  void ParallelDomainBase::solveInitialProblem(AdaptInfo *adaptInfo)
  {     
    if (timeIF)
      timeIF->solveInitialProblem(adaptInfo);
    
    synchVector(*(probStat->getSolution()));
  }


695
  void ParallelDomainBase::createInteriorBoundaryInfo()
696
  {
697
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
698

Thomas Witkowski's avatar
Thomas Witkowski committed
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
796
797
798
799
800
801
802
803
804
805
806
807
808
809
    // === 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()");

810
811
812
813
814
815
816
817
818
819
820
821
822
    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");
    }     
823
824
825

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

827
    TraverseStack stack;
828
829
830
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
831
832
833
834
    while (elInfo) {
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
835

836
      // Check, if the element is within rank's partition.
837
      if (partitionData->getPartitionStatus() == IN) {
838
	for (int i = 0; i < nNeighbours; i++) {
839
840
841
842
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
843
844
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->
						getElementData(PARTITION_ED));
845

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

848
849
850
851
852
	    // 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. ===

853
	    AtomicBoundary bound;	    	    
854
855
	    bound.rankObj.el = element;
	    bound.rankObj.elIndex = element->getIndex();
856
	    bound.rankObj.elType = elInfo->getType();
857
	    bound.rankObj.subObj = subObj;
858
	    bound.rankObj.ithObj = i;
859
860

	    bound.neighObj.el = elInfo->getNeighbour(i);
861
	    bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex();
862
	    bound.neighObj.elType = -1;
863
	    bound.neighObj.subObj = subObj;
864
	    bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i);
865
	    
866
867
868
869
870
	    if (dim == 2) {
	      // i == 2  =>  getSideOfNeighbour(i) == 2
	      TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
		("Should not happen!\n");
	    }
871

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

875
876
877
878
879
880
881
882
883

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

884
	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(subObj, i))) {	      
885
886
887
888
889
	      // 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.
890
891
892
893
894
895
896
897
898
	      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);
	      }	      
899
	    }
900
901
902
903
904
905
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
906
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
907
908


Thomas Witkowski's avatar
Thomas Witkowski committed
909
910
911
  void ParallelDomainBase::createSubstructureInteriorBoundaryInfo()
  {
    FUNCNAME("ParallelDomainBase::createSubstructureInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
912

Thomas Witkowski's avatar
Thomas Witkowski committed
913
914
915
916
917
    // === 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();
918
    const BasisFunction *basFcts = feSpace->getBasisFcts();
Thomas Witkowski's avatar
Thomas Witkowski committed
919
920
921
922
923
924
925
    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
926

927
928
    // 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
929
930
    // Maps each edge in ranks partition of the domain to the element object that 
    // contains this edge.
931
    std::map<GlobalEdge, BoundaryObject> rankEdges;
932

933

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

Thomas Witkowski's avatar
Thomas Witkowski committed
936
937
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
938
939
940
941
942
943
944
    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
945
946
947
948
949
950
951
952
953
954
955
956
957
958
      // 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);
959
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
960
961
962
963
964
965
966
967
968
969
      }

      
      // 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);
970
971
972
973
974
      }      

      elInfo = stack.traverseNext(elInfo);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
975
976

    // === Create a set of all edges and vertices at rank's interior boundaries. ===
977
978
979

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

Thomas Witkowski's avatar
Thomas Witkowski committed
983
984
   
    // First, traverse the rank owned elements at the interior boundaries.
Thomas Witkowski's avatar
Thomas Witkowski committed
985
986
987
    for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
      Element *el = it.getBound().rankObj.el;
      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
988
      
Thomas Witkowski's avatar
Thomas Witkowski committed
989
      if (dim == 3) {
990
	for (int j = 0; j < 3; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
991
992
993
994
995
996
997
998
999
1000
	  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)
For faster browsing, not all history is shown. View entire blame