MeshDistributor.cc 73.4 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 <boost/filesystem.hpp>

7
8
#include "parallel/MeshDistributor.h"
#include "parallel/ParallelDebug.h"
9
#include "parallel/StdMpi.h"
10
11
12
13
14
15
16
#include "ParMetisPartitioner.h"
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
#include "PartitionElementData.h"
17
18
#include "DOFMatrix.h"
#include "DOFVector.h"
19
#include "SystemVector.h"
20
#include "VtkWriter.h"
21
#include "ElementDofIterator.h"
22
23
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
24
#include "ElementFileWriter.h"
25
#include "VertexVector.h"
26
#include "MeshStructure.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
27
28
#include "ProblemVec.h"
#include "ProblemInstat.h"
29
#include "Debug.h"
30

31
32
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
33
  using boost::lexical_cast;
34
  using namespace boost::filesystem;
Thomas Witkowski's avatar
Thomas Witkowski committed
35

36
37
38
39
40
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

41

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

61
62
63
64
65
    mpiRank = MPI::COMM_WORLD.Get_rank();
    mpiSize = MPI::COMM_WORLD.Get_size();
    mpiComm = MPI::COMM_WORLD;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
66

67
  void MeshDistributor::initParallelization(AdaptInfo *adaptInfo)
68
  {
69
    FUNCNAME("MeshDistributor::initParallelization()");
70
71
72
73

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

74
75
    TEST_EXIT(mesh)("No mesh has been defined for mesh distribution!\n");

76
77
78
    // If the problem has been already read from a file, we need only to set
    // isRankDofs to all matrices and rhs vector and to remove periodic 
    // boundary conditions (if there are some).
79
    if (deserialized) {
80
      setRankDofs();
81
      removePeriodicBoundaryConditions();
82

83
      return;
84
    }
85
    
86

87
88
89
90
91
    // 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();

92

93
94
95
96
97
98
99
    // 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);   

100

Thomas Witkowski's avatar
Thomas Witkowski committed
101
#if (DEBUG != 0)
102
103
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
104
105
106
    if (mpiRank == 0) {
      int writePartMesh = 1;
      GET_PARAMETER(0, "dbg->write part mesh", "%d", &writePartMesh);
107

108
109
110
111
112
      if (writePartMesh > 0)
	writePartitioningMesh("part.vtu");
      else 
	MSG("Skip write part mesh!\n");
    }
113

114
    ParallelDebug::testAllElements(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
115
#endif
116

117

118
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
119

120
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
121

122
#if (DEBUG != 0)
123
    ParallelDebug::printBoundaryInfo(*this);
124
125
#endif

126

127
128
129
130
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

131

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
132
133
134
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
135
136
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
137

138
139
140
141
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

142

143
144
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
145
#if (DEBUG != 0)
146
    MSG("AMDiS runs in debug mode, so make some test ...\n");
147
    debug::testSortedDofs(mesh, elMap);
148
149
    ParallelDebug::testInteriorBoundary(*this);
    ParallelDebug::testCommonDofs(*this, true);
150
    MSG("Debug mode tests finished!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
151

152
    debug::writeMesh(feSpace, -1, "macro_mesh");   
153
#endif
154

155

156
157
    // === Create periodic dof mapping, if there are periodic boundaries. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
158
    createPeriodicMap();    
159

160

161
    // === Global refinements. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
162

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

Thomas Witkowski's avatar
Thomas Witkowski committed
166
    if (globalRefinement > 0) {
167
      refineManager->globalRefine(mesh, globalRefinement);
168

169
#if (DEBUG != 0)
170
      debug::writeMesh(feSpace, -1, "gr_mesh");
171
172
#endif

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

Thomas Witkowski's avatar
Thomas Witkowski committed
181
182
    /// === Set DOF rank information to all matrices and vectors. ===

183
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
184

185

Thomas Witkowski's avatar
Thomas Witkowski committed
186
187
    // === Remove periodic boundary conditions in sequential problem definition. ===

188
    removePeriodicBoundaryConditions();
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
223
224
225
226
227
  void MeshDistributor::addProblemStat(ProblemVec *probVec)
  {
    FUNCNAME("MeshDistributor::addProblemVec()");

    if (feSpace != NULL) {
      std::vector<FiniteElemSpace*> feSpaces = probVec->getFeSpaces();
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	TEST_EXIT(feSpace == feSpaces[i])
	  ("Parallelizaton is not supported for multiple FE spaces!\n");
      }
    } else {
      feSpace = probVec->getFeSpace(0);
      mesh = feSpace->getMesh();
      info = probVec->getInfo();
      
      TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1)
	("Only meshes with one DOFAdmin are supported!\n");
      TEST_EXIT(mesh->getDOFAdmin(0).getNumberOfPreDOFs(0) == 0)
	("Wrong pre dof number for DOFAdmin!\n");
      
      switch (mesh->getDim()) {
      case 2:
	refineManager = new RefinementManager2d();
	break;
      case 3:
	refineManager = new RefinementManager3d();
	break;
      default:
	ERROR_EXIT("This should not happen for dim = %d!\n", mesh->getDim());
      }

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

    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
228
229
230
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
231
232
233
234
235
      std::string filename = "";
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
236
237
238
239
240

      int tsModulo = -1;
      GET_PARAMETER(0, probVec->getName() + "->output->write every i-th timestep", 
		    "%d", &tsModulo);
      
241
      probVec->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
242
243
      writeSerializationFile = true;
    }    
244
245

    int readSerialization = 0;
246
247
    GET_PARAMETER(0, probVec->getName() + "->input->read serialization", "%d", 
		  &readSerialization);
248
249
250
251
    if (readSerialization) {
      std::string filename = "";
      GET_PARAMETER(0, probVec->getName() + "->input->serialization filename", &filename);
      filename += ".p" + lexical_cast<std::string>(mpiRank);
252
      MSG("Start deserialization with %s\n", filename.c_str());
253
      std::ifstream in(filename.c_str());
254
255
256
257

      TEST_EXIT(!in.fail())("Could not open deserialization file: %s\n",
			    filename.c_str());

258
      probVec->deserialize(in);
259
      in.close();
260
261
262
      MSG("Deserialization from file: %s\n", filename.c_str());

      if (!deserialized) {
263
264
265
266
267
268
269
270
271
272
273
274
	filename = "";
	GET_PARAMETER(0, name + "->input->serialization filename", &filename);

	TEST_EXIT(filename != "")
	  ("No filename defined for parallel deserialization file!\n");

	std::string rankFilename = filename + ".p" + lexical_cast<std::string>(mpiRank);
	std::ifstream in(rankFilename.c_str());

	TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			      filename.c_str());

275
276
	deserialize(in);
	in.close();
277
	MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
278
279
	deserialized = true;
      }
280
281
282
283
284
285
286
    }

    probStat.push_back(probVec);
  }


  void MeshDistributor::exitParallelization(AdaptInfo *adaptInfo)
287
  {}
288

289
  
290
  void MeshDistributor::updateDofAdmins()
291
  {
292
    FUNCNAME("MeshDistributor::updateDofAdmins()");
293
294

    for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
295
      DOFAdmin& admin = const_cast<DOFAdmin&>(mesh->getDOFAdmin(i));
296
297
298

      // There must be always more allocated DOFs than used DOFs in DOFAdmin. Otherwise,
      // it is not possible to define a first DOF hole.
299
      if (static_cast<unsigned int>(admin.getSize()) == mapLocalGlobalDofs.size())
300
	admin.enlargeDOFLists();
301
302
303
      
      for (int j = 0; j < admin.getSize(); j++)
	admin.setDOFFree(j, true);
304
      for (unsigned int j = 0; j < mapLocalGlobalDofs.size(); j++)
305
306
 	admin.setDOFFree(j, false);

307
308
309
      admin.setUsedSize(mapLocalGlobalDofs.size());
      admin.setUsedCount(mapLocalGlobalDofs.size());
      admin.setFirstHole(mapLocalGlobalDofs.size());
310
311
312
    }
  }

313

314
  void MeshDistributor::testForMacroMesh()
315
  {
316
    FUNCNAME("MeshDistributor::testForMacroMesh()");
317
318
319
320
321
322
323

    int nMacroElements = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      TEST_EXIT(elInfo->getLevel() == 0)
324
	("Mesh is already refined! This does not work with parallelization!\n");
325
326
327
328
329
330
331
332
333
334
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

335

336
  void MeshDistributor::synchVector(DOFVector<double> &vec)
337
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
338
    StdMpi<std::vector<double> > stdMpi(mpiComm);
339
340

    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
341
	 sendIt != sendDofs.end(); ++sendIt) {
342
      std::vector<double> dofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
      int nSendDofs = sendIt->second.size();
      dofs.reserve(nSendDofs);
345
      
Thomas Witkowski's avatar
Thomas Witkowski committed
346
347
      for (int i = 0; i < nSendDofs; i++)
	dofs.push_back(vec[*((sendIt->second)[i])]);
348
349
350
351

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

Thomas Witkowski's avatar
Thomas Witkowski committed
352
353
354
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt)
      stdMpi.recv(recvIt->first, recvIt->second.size());
355

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

Thomas Witkowski's avatar
Thomas Witkowski committed
358
359
360
361
362
    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];
  }
363
364


365
  void MeshDistributor::synchVector(SystemVector &vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
366
  {
367
    int nComponents = vec.getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
368
369
370
371
372
373
374
375
376
377
378
379
    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])]);
380
381
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
382
      stdMpi.send(sendIt->first, dofs);
383
384
385
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
389
    stdMpi.startCommunication<double>(MPI_DOUBLE);
390
391

    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
	 recvIt != recvDofs.end(); ++recvIt) {
      int nRecvDofs = recvIt->second.size();
394
395

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
396
397
398
399
400
      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++];
401
402
403
404
      }
    }
  }

405

406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
  void MeshDistributor::setRankDofs()
  {
    for (unsigned int i = 0; i < probStat.size(); i++) {
      int nComponents = probStat[i]->getNumComponents();
      for (int j = 0; j < nComponents; j++) {
	for (int k = 0; k < nComponents; k++)
	  if (probStat[i]->getSystemMatrix(j, k))
	    probStat[i]->getSystemMatrix(j, k)->setRankDofs(isRankDof);

	TEST_EXIT_DBG(probStat[i]->getRhs()->getDOFVector(j))("No RHS vector!\n");
	TEST_EXIT_DBG(probStat[i]->getSolution()->getDOFVector(j))("No solution vector!\n");
	
	probStat[i]->getRhs()->getDOFVector(j)->setRankDofs(isRankDof);
	probStat[i]->getSolution()->getDOFVector(j)->setRankDofs(isRankDof);
      }
    }
  }


425
426
427
428
429
430
431
432
433
434
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
    // Remove periodic boundaries in boundary manager on matrices and vectors.
    for (unsigned int i = 0; i < probStat.size(); i++) {
      int nComponents = probStat[i]->getNumComponents();

      for (int j = 0; j < nComponents; j++) {
	for (int k = 0; k < nComponents; k++) {
	  DOFMatrix* mat = probStat[i]->getSystemMatrix(j, k);
	  if (mat && mat->getBoundaryManager())
435
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
	}
	
	if (probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat[i]->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
	
	if (probStat[i]->getRhs()->getDOFVector(j)->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat[i]->getRhs()->getDOFVector(j)->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);
    }    
  }


  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
457
458
459
460
461
462
463
464
465
466
467
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


468
  void MeshDistributor::checkMeshChange()
469
  {
470
    FUNCNAME("MeshDistributor::checkMeshChange()");
471

472
473
474
475
#if (DEBUG != 0)
    debug::writeMesh(feSpace, -1, "before_check_mesh");
#endif

476
477
478
479
480
    // === 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);
481

482
    if (recvAllValues == 0)
483
484
      return;

485
486
    // === At least one rank mesh has been changed, so the boundaries must be ===
    // === adapted to the new mesh structure.                                 ===
487

488
489
490
491
492
    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.
493
      RankToBoundMap allBound;
Thomas Witkowski's avatar
Thomas Witkowski committed
494
495

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

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

503
504
505
506
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
	  allBound[it.getRank()].push_back(*it);	

507

508
      // === Check the boundaries and adapt mesh if necessary. ===
509

510
511
512
513
#if (DEBUG != 0)
      MSG("Run checkAndAdaptBoundary ...\n");
#endif

514
515
516
517
518
519
520
      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);
521
522
523
524

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

527
528

#if (DEBUG != 0)
529
    debug::writeMesh(feSpace, -1, "mesh");
530
531
532
533
534
535
#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. ===
536

537
    mesh->dofCompress();
538
    updateLocalGlobalNumbering();
539
540
541
542

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

    createPeriodicMap();
543
544
545
  }

  
546
  bool MeshDistributor::checkAndAdaptBoundary(RankToBoundMap &allBound)
547
  {
548
    FUNCNAME("MeshDistributor::checkAndAdaptBoundary()");
549
550
551
552
553
554

    // === Create mesh structure codes for all ranks boundary elements. ===
       
    std::map<int, MeshCodeVec> sendCodes;
   
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
555
556
557
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	MeshStructure elCode;
558
	elCode.init(boundIt->rankObj);
559
560
561
562
	sendCodes[it->first].push_back(elCode);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
563
    StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
564
    stdMpi.send(sendCodes);
565
566
567
    stdMpi.recv(allBound);
    stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG);
 
568
    // === Compare received mesh structure codes. ===
569
    
570
571
    bool meshFitTogether = true;

572
573
    for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) {
      
574
575
      MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first];
      int i = 0;
576
      
577
578
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
579

580
581
582
	MeshStructure elCode;	
	elCode.init(boundIt->rankObj);
	
583
584
	if (elCode.getCode() != recvCodes[i].getCode()) {
	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
585
586

	  bool b = fitElementToMeshCode(recvCodes[i], 
587
					boundIt->rankObj.el,
588
					boundIt->rankObj.subObj,
589
					boundIt->rankObj.ithObj, 
590
591
					boundIt->rankObj.elType,
					boundIt->rankObj.reverseMode);  
592

593
594
	  if (b)
	    meshFitTogether = false;
595
	}
596
	
597
	i++;
598
599
600
      }
    }

601
    return meshFitTogether;
602
  }
603
604


605
  bool MeshDistributor::fitElementToMeshCode(MeshStructure &code, 
606
607
608
609
610
					     Element *el, 
					     GeoIndex subObj,
					     int ithObj, 
					     int elType,
					     bool reverseMode)
611
  {
612
    FUNCNAME("MeshDistributor::fitElementToMeshCode()");
613

614
615
    TEST_EXIT_DBG(el)("No element given!\n");

616
617
    if (code.empty())
      return false;
618
619
620

    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
621
622
623

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

624
    bool meshChanged = false;
625
626
627
628
629
630
631
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

    if (reverseMode) {
      std::swap(s1, s2);
      traverseFlag |= Mesh::CALL_REVERSE_MODE;
    }
632
633

    if (s1 != -1 && s2 != -1) {
634
      TraverseStack stack;
635
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
636
637
      while (elInfo && elInfo->getElement() != el)
	elInfo = stack.traverseNext(elInfo);      
638

639
      meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode);
640
641
      return meshChanged;
    }
642
643

    if (el->isLeaf()) {
644
645
646
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

647
648
649
      meshChanged = true;

      TraverseStack stack;
650
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
651
652
653
654
655
656
657

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

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

658
      el->setMark(1);
659
660
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
661
      refineManager->refineFunction(elInfo);
662
    }
663

664
665
666
667
668
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode)
      std::swap(child0, child1);

669
    if (s1 != -1) {
670
      TraverseStack stack;
671
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
672

673
      while (elInfo && elInfo->getElement() != child0) {
674
675
676
	elInfo = stack.traverseNext(elInfo);
      }

677
      meshChanged |= 
678
	fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
679
    } else {
680
      TraverseStack stack;
681
      ElInfo *elInfo = stack.traverseFirst(el->getMesh(), -1, traverseFlag);
682

683
      while (elInfo && elInfo->getElement() != child1) {
684
685
686
	elInfo = stack.traverseNext(elInfo);
      }

687
      meshChanged |= 
688
	fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
689
690
691
    }

    return meshChanged;
692
693
  }

694

695
  bool MeshDistributor::fitElementToMeshCode2(MeshStructure &code, 
696
697
698
699
700
					      TraverseStack &stack,
					      GeoIndex subObj,
					      int ithObj, 
					      int elType,
					      bool reverseMode)
701
  {
702
    FUNCNAME("MeshDistributor::fitElementToMeshCode2()");
703

704
705
    ElInfo *elInfo = stack.getElInfo();

706
    bool value = false;
707
708
709
710
711
712
713
714
715
716
717
    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);
718
     
719
720
      return value;
    }
721

722
723
    if (!elInfo)
      return value;
724

725
    if (el->isLeaf()) {
726
      el->setMark(1);
727
728
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
729
      refineManager->refineFunction(elInfo);
730
731
732
      value = true;
    }

733
734
    int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
    int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);
735
736
737
738
739
740
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
      std::swap(s1, s2);
      std::swap(child0, child1);
    }
741

742
    if (s1 != -1) {
743
      stack.traverseNext(elInfo);
744
      code.nextElement();
745
      value |= fitElementToMeshCode2(code, stack, subObj, s1, el->getChildType(elType), reverseMode);
746
747
748
749
      elInfo = stack.getElInfo();
    } else {
      do {
	elInfo = stack.traverseNext(elInfo);
750
      } while (elInfo && elInfo->getElement() != child1);      
751
752
    }  

753
    TEST_EXIT_DBG(elInfo->getElement() == child1)("This should not happen!\n");
754

755
    if (s2 != -1) {
756
      code.nextElement();
757
      value |= fitElementToMeshCode2(code, stack, subObj, s2, el->getChildType(elType), reverseMode);
758
759
    } else {
      int level = elInfo->getLevel();
760

761
762
763
764
      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }
765

766
767
768
    return value;
  }

769
  
770
  void MeshDistributor::serialize(std::ostream &out, DofContainer &data)
Thomas Witkowski's avatar
Thomas Witkowski committed
771
  {    
772
    int vecSize = data.size();
773
    SerUtil::serialize(out, vecSize);
774
775
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = (*data[i]);
776
      SerUtil::serialize(out, dofIndex);
777
778
779
780
    }
  }


781
  void MeshDistributor::deserialize(std::istream &in, DofContainer &data,
782
				    std::map<int, const DegreeOfFreedom*> &dofMap)
783
  {
784
    FUNCNAME("MeshDistributor::deserialize()");
785
786

    int vecSize = 0;
787
    SerUtil::deserialize(in, vecSize);
788
789
790
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
791
      SerUtil::deserialize(in, dofIndex);
792
793
794
795
796
797
798
799
800

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

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


801
  void MeshDistributor::serialize(std::ostream &out, RankToDofContainer &data)
802
803
  {
    int mapSize = data.size();
804
    SerUtil::serialize(out, mapSize);
805
806
    for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) {
      int rank = it->first;
807
      SerUtil::serialize(out, rank);
808
809
810
811
812
      serialize(out, it->second);
    }
  }

  
813
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
814
				    std::map<int, const DegreeOfFreedom*> &dofMap)
815
816
  {
    int mapSize = 0;
817
    SerUtil::deserialize(in, mapSize);
818
819
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
820
      SerUtil::deserialize(in, rank);
821
822
823
824
      deserialize(in, data[rank], dofMap);      
    }
  }

825

826
  double MeshDistributor::setElemWeights(AdaptInfo *adaptInfo) 
827
  {
828
    FUNCNAME("MeshDistributor::setElemWeights()");
829

830
    double localWeightSum = 0.0;
831
    elemWeights.clear();
832
833
834
835
836
837
838
839
840
841
842
843
844
845
      
    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;
846

847
848
849
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
850

851
      infile.close();
852
    } else {           
853
      TraverseStack stack;
854
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
855
      while (elInfo) {
856
857
858
	elemWeights[elInfo->getElement()->getIndex()] = 1.0;
	localWeightSum++;

859
	elInfo = stack.traverseNext(elInfo);
860
861
862
863
864
865
      }
    }

    return localWeightSum;
  }

866

867
  void MeshDistributor::partitionMesh(AdaptInfo *adaptInfo)
868
  {
869
    FUNCNAME("MeshDistributor::partitionMesh()");
870