ParallelDomainBase.cc 66.3 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

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

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

    // 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();
    }
92
93
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
94

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

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

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
#if (DEBUG != 0)
119
120
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, 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");
    }
130
131

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

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

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

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

    createLocalGlobalNumbering();

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

    removeMacroElements();

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

    updateDofAdmins();

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

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

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

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

    createPeriodicMap();
165

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

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

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

174
175
176
177
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "grmesh");
#endif

178
      updateLocalGlobalNumbering();
179

180
      // === Update periodic mapping, if there are periodic boundaries. ===
181

182
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
183
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
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
219
220
221
222

    /// === 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);
    }    
223
224
  }

225
226

  void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo)
227
  {}
228

229
230
231
  
  void ParallelDomainBase::updateDofAdmins()
  {
232
233
234
    FUNCNAME("ParallelDomainBase::updateDofAdmins()");

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
235
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
236
237
238

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

247
248
249
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
250
251
252
    }
  }

253

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

      elInfo = stack.traverseNext(elInfo);
    }

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

275

Thomas Witkowski's avatar
Thomas Witkowski committed
276
  void ParallelDomainBase::synchVector(DOFVector<double> &vec)
277
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    StdMpi<std::vector<double> > stdMpi(mpiComm);
279
280

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
292
293
294
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
295

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

Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
300
301
302
    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];
  }
303
304


Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
307
308
309
310
311
312
313
314
315
316
317
318
  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])]);
319
320
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
321
      stdMpi.send(sendIt->first, dofs);
322
323
324
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
328
    stdMpi.startCommunication<double>(MPI_DOUBLE);
329
330

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
331
332
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
333
334

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
337
338
339
      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++];
340
341
342
343
      }
    }
  }

344

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


357
358
  void ParallelDomainBase::checkMeshChange()
  {
359
360
    FUNCNAME("ParallelDomainBase::checkMeshChange()");

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

363
364
365
366
367
    // === 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);
368

369
    if (recvAllValues == 0)
370
371
      return;

372
373
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
374

375
376
377
378
379
    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.
380
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
381
382

      for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it)
383
384
	if (it->rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
385
386

      for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it)
387
388
	if (it->rankObj.subObj == INDEX_OF_DIM(dim - 1, dim))
	  allBound[it.getRank()].push_back(*it);
Thomas Witkowski's avatar
Thomas Witkowski committed
389

390

391
      // === Check the boundaries and adapt mesh if necessary. ===
392

393
394
395
396
397
398
399
400
      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);
401

402
403

#if (DEBUG != 0)
404
    debug::writeMesh(feSpace, -1, "mesh");
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
#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) {
425
426
427
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
428
	elCode.init(boundIt->rankObj);
429
430
431
432
	sendCodes[it->first].push_back(elCode);
      }
    }

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

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

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

	  bool b = fitElementToMeshCode(recvCodes[i], 
457
458
					boundIt->rankObj.el,
					boundIt->rankObj.ithObj, 
459
					boundIt->rankObj.elType);  
460

461
462
	  if (b)
	    meshFitTogether = false;
463
	}
464
	
465
	i++;
466
467
468
      }
    }

469
    return meshFitTogether;
470
  }
471
472
473
474
475


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

480
481
    TEST_EXIT_DBG(el)("No element given!\n");

482
483
484
485
486
487
488
489
    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");

490
491
492
    bool meshChanged = false;

    if (s1 != -1 && s2 != -1) {
493
494
495
496
497
498
499
      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);
      }

500
      meshChanged = fitElementToMeshCode2(code, stack, ithSide, elType);
501
502
      return meshChanged;
    }
503
504

    if (el->isLeaf()) {
505
506
507
508
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

      ERROR_EXIT("NOT YET WORKING!\n");
509
      
510
511
      ElInfo *elInfo = el->getMesh()->createNewElInfo();    
      //      el->getMesh()->fillElInfo(elInfo, el, Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
512
      el->setMark(1);
513
      refineManager->refineFunction(elInfo);
514
515
      delete elInfo;

516
      meshChanged = true;
517
    }
518
519

    if (s1 != -1) {
520
521
522
523
524
525
526
527
      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);
      }

528
      meshChanged |= fitElementToMeshCode2(code, stack, s1, el->getChildType(elType));
529
    } else {
530
531
532
533
534
535
536
537
      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);
      }

538
      meshChanged |= fitElementToMeshCode2(code, stack, s2, el->getChildType(elType));
539
540
541
    }

    return meshChanged;
542
543
  }

544

545
  bool ParallelDomainBase::fitElementToMeshCode2(MeshStructure &code, 
546
						 TraverseStack &stack,
547
						 int ithSide, 
548
						 int elType)
549
  {
550
551
    FUNCNAME("ParallelDomainBase::fitElementToMeshCode2()");

552
553
    ElInfo *elInfo = stack.getElInfo();

554
    bool value = false;
555
556
557
558
559
560
561
562
563
564
565
    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);
566
     
567
568
      return value;
    }
569

570
571
    if (!elInfo)
      return value;
572

573
    if (el->isLeaf()) {
574
      el->setMark(1);
575
576
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
577
      refineManager->refineFunction(elInfo);
578
579
580
581
582
      value = true;
    }

    int s1 = el->getSideOfChild(0, ithSide, elType);
    int s2 = el->getSideOfChild(1, ithSide, elType);
583

584
    if (s1 != -1) {
585
      stack.traverseNext(elInfo);
586
      code.nextElement();
587
      value |= fitElementToMeshCode2(code, stack, s1, el->getChildType(elType));
588
589
590
591
592
      elInfo = stack.getElInfo();
    } else {
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getElement() != el->getSecondChild());      
593
594
    }  

595
596
597
    TEST_EXIT_DBG(elInfo->getElement() == el->getSecondChild())
      ("This should not happen!\n");

598
    if (s2 != -1) {
599
      code.nextElement();
600
      value |= fitElementToMeshCode2(code, stack, s2, el->getChildType(elType));
601
602
    } else {
      int level = elInfo->getLevel();
603

604
605
606
607
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
608

609
610
611
    return value;
  }

612
613
  
  void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
614
  {    
615
    int vecSize = data.size();
616
    SerUtil::serialize(out, vecSize);
617
618
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
619
      SerUtil::serialize(out, dofIndex);
620
621
622
623
624
625
626
627
628
629
    }
  }


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

    int vecSize = 0;
630
    SerUtil::deserialize(in, vecSize);
631
632
633
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
634
      SerUtil::deserialize(in, dofIndex);
635
636
637
638
639
640
641
642
643
644
645
646

      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();
647
    SerUtil::serialize(out, mapSize);
648
649
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
650
      SerUtil::serialize(out, rank);
651
652
653
654
655
656
657
658
659
      serialize(out, it->second);
    }
  }

  
  void ParallelDomainBase::deserialize(std::istream &in, RankToDofContainer &data,
				       std::map<int, const DegreeOfFreedom*> &dofMap)
  {
    int mapSize = 0;
660
    SerUtil::deserialize(in, mapSize);
661
662
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
663
      SerUtil::deserialize(in, rank);
664
665
666
667
      deserialize(in, data[rank], dofMap);      
    }
  }

668

669
  double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) 
670
  {
671
    FUNCNAME("ParallelDomainBase::setElemWeights()");
672

673
    double localWeightSum = 0.0;
674
    elemWeights.clear();
675
676
677
678
679
680
681
682
683
684
685
686
687
688
      
    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;
689

690
691
692
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
693

694
      infile.close();
695
    } else {           
696
      TraverseStack stack;
697
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
698
      while (elInfo) {
699
700
701
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
	localWeightSum++;

702
	elInfo = stack.traverseNext(elInfo);
703
704
705
706
707
708
      }
    }

    return localWeightSum;
  }

709

710
  void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo)
711
  {
712
713
    FUNCNAME("ParallelDomainBase::partitionMesh()");

714
715
716
717
718
719
720
721
722
723
724
725
    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);
  }

726

Thomas Witkowski's avatar
Thomas Witkowski committed
727
728
729
730
731
732
733
734
735
  void ParallelDomainBase::solveInitialProblem(AdaptInfo *adaptInfo)
  {     
    if (timeIF)
      timeIF->solveInitialProblem(adaptInfo);
    
    synchVector(*(probStat->getSolution()));
  }


736
  void ParallelDomainBase::createInteriorBoundaryInfo()
737
  {
738
    FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
739

Thomas Witkowski's avatar
Thomas Witkowski committed
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
    // === 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);
759

Thomas Witkowski's avatar
Thomas Witkowski committed
760
761
762
763
764
765
766
767
    // === 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 ===
768
      // === indices. Now, we have to sort the corresponding list in this rank to   ===
Thomas Witkowski's avatar
Thomas Witkowski committed
769
      // === get the same order.                                                    ===
770
     
Thomas Witkowski's avatar
Thomas Witkowski committed
771
      for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
772

Thomas Witkowski's avatar
Thomas Witkowski committed
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
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
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
	// 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()");

852
853
    int nNeighbours = mesh->getGeo(NEIGH);
    int dim = mesh->getDim();
854
    GeoIndex subObj = INDEX_OF_DIM(dim - 1, dim);
855
856
857

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

859
    TraverseStack stack;
860
861
862
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
863
864
865
    while (elInfo) {
      Element *element = elInfo->getElement();
      PartitionElementData *partitionData = 
866
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));
867

868
      // Check, if the element is within rank's partition.
869
      if (partitionData->getPartitionStatus() == IN) {
870
	for (int i = 0; i < nNeighbours; i++) {
871
872
873
874
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
875
876
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->
						getElementData(PARTITION_ED));
877

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

880
881
882
883
884
	    // 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. ===

885
	    AtomicBoundary bound;	    	    
886
887
	    bound.rankObj.el = element;
	    bound.rankObj.elIndex = element->getIndex();
888
	    bound.rankObj.elType = elInfo->getType();
889
	    bound.rankObj.subObj = subObj;
890
	    bound.rankObj.ithObj = i;
891
892

	    bound.neighObj.el = elInfo->getNeighbour(i);
893
	    bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex();
894
	    bound.neighObj.elType = -1;
895
	    bound.neighObj.subObj = subObj;
896
	    bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i);
897
	    
898
	    // Get rank number of the neighbouring element.
899
900
	    int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];

901
902
903
904
905
906
907
908
	    // === 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.                                         ===

909
	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(subObj, i))) {	      
910
911
912
913
914
	      // 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.
915
916
917
918
919
920
921
922
923
	      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);
	      }	      
924
	    }
925
926
927
928
929
930
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
931
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
932
933


Thomas Witkowski's avatar
Thomas Witkowski committed
934
935
936
  void ParallelDomainBase::createSubstructureInteriorBoundaryInfo()
  {
    FUNCNAME("ParallelDomainBase::createSubstructureInteriorBoundaryInfo()");
Thomas Witkowski's avatar
Thomas Witkowski committed
937

938
939
940
    // === 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
941
942

    int dim = mesh->getDim();
943
    const BasisFunction *basFcts = feSpace->getBasisFcts();
944
    std::vector<DegreeOfFreedom> localIndices(basFcts->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
945
946
947

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

951
952
    // Maps each edge in the whole domain to the rank that ownes this edge.
    std::map<GlobalEdge, int> edgeOwner;
953
    // Maps each edge in ranks partition to the element object that contains this edge.
954
    std::map<GlobalEdge, BoundaryObject> rankEdges;
955

956

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

Thomas Witkowski's avatar
Thomas Witkowski committed
959
960
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
961
962
963
964
965
966
    while (elInfo) {
      Element *el = elInfo->getElement();
      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);

      PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	(el->getElementData(PARTITION_ED));
967
968
          
      // In 2d and 3d, traverse all vertices of current element.
969
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
970
971
972
973
974
975
976
977
978
979
	DegreeOfFreedom dof0 = localIndices[i];
	// Update the owner of the current vertex DOF.
	dofOwner[dof0] = max(dofOwner[dof0], partitionVec[el->getIndex()]);

	// If the DOF is part of an element that is within ranks partition, add it
	// to the set of all ranks vertex DOFs.
	if (partitionData && partitionData->getPartitionStatus() == IN)
	  rankDofs[dof0] = BoundaryObject(el, elInfo->getType(), VERTEX, i);
      }      

Thomas Witkowski's avatar
Thomas Witkowski committed
980
981
982
983
984
985
986
987
988
989
      // 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()]);
	  
990
991
	  // If the edge is part of an element that is within ranks partition, add it
	  // to the set of all ranks edges.
Thomas Witkowski's avatar
Thomas Witkowski committed
992
993
	  if (partitionData && partitionData->getPartitionStatus() == IN)
	    rankEdges[edge] = BoundaryObject(el, elInfo->getType(), EDGE, i);
994
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
995
996
      }

997
998
999
      elInfo = stack.traverseNext(elInfo);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1000

1001
    // === Create a set of all edges and vertices at ranks interior boundaries. ===
1002
1003
1004

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1008
1009
   
    // First, traverse the rank owned elements at the interior boundaries.
Thomas Witkowski's avatar
Thomas Witkowski committed
1010
    for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) {
1011
      Element *el = it->rankObj.el;
Thomas Witkowski's avatar
Thomas Witkowski committed
1012
      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
1013
      
1014
1015
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
	int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, i);
1016
1017
1018
1019
1020
1021
1022
1023
	DegreeOfFreedom dof = localIndices[dofNo];
	
	if (dofOwner[dof] > mpiRank)
	  it->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
	else
	  rankBoundaryDofs.insert(dof);
      }    

Thomas Witkowski's avatar
Thomas Witkowski committed
1024
      if (dim == 3) {
1025
1026
	for (int i = 0; i < 3; i++) {
	  int edgeNo = el->getEdgeOfFace(it->rankObj.ithObj, i);
Thomas Witkowski's avatar
Thomas Witkowski committed
1027
1028
1029
1030
1031
1032
1033
1034
1035
	  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)
1036
	    it->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo));
1037
	  else
Thomas Witkowski's avatar
Thomas Witkowski committed
1038
	    rankBoundaryEdges.insert(edge);
1039
	}
1040
      }     
Thomas Witkowski's avatar
Thomas Witkowski committed
1041
1042
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1043
1044
    // Now, do the same with all other elements at the interior boundaries.
    for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) {
1045
      Element *el = it->rankObj.el;
Thomas Witkowski's avatar
Thomas Witkowski committed
1046
1047
      basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
      
1048
1049
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
	int dofNo = el->getVertexOfPosition(it->rankObj.subObj, it->rankObj.ithObj, i);
1050
1051
1052
1053
1054
1055
1056
1057
1058

	DegreeOfFreedom dof = localIndices[dofNo];
	
	if (dofOwner[dof] > it.getRank())
	  it->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo));
	else
	  rankBoundaryDofs.insert(dof);
      }	      

Thomas Witkowski's avatar
Thomas Witkowski committed
1059
      if (dim == 3) {