MeshDistributor.cc 67.6 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 <boost/filesystem.hpp>
6
7
#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_comparison.hpp>
8

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

33
34
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
35
  using boost::lexical_cast;
36
  using namespace boost::filesystem;
Thomas Witkowski's avatar
Thomas Witkowski committed
37

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

43

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
68

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

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

76
77
    TEST_EXIT(mesh)("No mesh has been defined for mesh distribution!\n");

78
79
80
    // 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).
81
    if (deserialized) {
82
      setRankDofs();
83
      removePeriodicBoundaryConditions();
84

85
      return;
86
    }
87
    
88

89
90
91
92
93
    // 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();

94

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

102

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

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

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

119

120
    // === Create interior boundary information. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
121

122
    createInteriorBoundaryInfo();
Thomas Witkowski's avatar
Thomas Witkowski committed
123

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

128

129
130
131
132
    // === Create new global and local DOF numbering. ===

    createLocalGlobalNumbering();

133

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

    removeMacroElements();
137
138
    macroElementStructureConsisten = true;

Thomas Witkowski's avatar
n  
Thomas Witkowski committed
139

140
141
142
143
    // === Reset all DOFAdmins of the mesh. ===

    updateDofAdmins();

144

145
146
    // === If in debug mode, make some tests. ===

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

154
    debug::writeMesh(feSpace, -1, "macro_mesh");   
155
#endif
156

157

158
159
    // === Create periodic dof mapping, if there are periodic boundaries. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
160
    createPeriodicMap();    
161

162

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

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

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

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

175
      updateLocalGlobalNumbering();
176
      
177
      // === Update periodic mapping, if there are periodic boundaries. ===
178
      
179
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
180
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
181

182

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

185
    setRankDofs();
Thomas Witkowski's avatar
Thomas Witkowski committed
186

187

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

190
    removePeriodicBoundaryConditions();
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
228
229
  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;
230
231
232
    GET_PARAMETER(0, probVec->getName() + "->output->write serialization", "%d", 
		  &writeSerialization);
    if (writeSerialization && !writeSerializationFile) {
233
234
235
236
237
      std::string filename = "";
      GET_PARAMETER(0, name + "->output->serialization filename", &filename);
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
238
239
240
241
242

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

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

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

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

264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
      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);
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      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
497
 	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
 	  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
      for (InteriorBoundary::iterator it(periodicBoundary); !it.end(); ++it)
504
505
 	if (it->rankObj.subObj == EDGE || it->rankObj.subObj == FACE)
 	  allBound[it.getRank()].push_back(*it);	
506

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
    data.clear();
789
790
791
    data.resize(vecSize);
    for (int i = 0; i < vecSize; i++) {
      int dofIndex = 0;
792
      SerUtil::deserialize(in, dofIndex);
793
794
795
796
797
798
799
800
801

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

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


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

  
814
  void MeshDistributor::deserialize(std::istream &in, RankToDofContainer &data,
815
				    std::map<int, const DegreeOfFreedom*> &dofMap)
816
  {
817
818
    data.clear();

819
    int mapSize = 0;
820
    SerUtil::deserialize(in, mapSize);
821
822
    for (int i = 0; i < mapSize; i++) {
      int rank = 0;
823
      SerUtil::deserialize(in, rank);
824
825
826
827
      deserialize(in, data[rank], dofMap);      
    }
  }

828

829
  double MeshDistributor::setElemWeights(AdaptInfo *adaptInfo) 
830
  {
831
    FUNCNAME("MeshDistributor::setElemWeights()");
832

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

850
851
852
	elemWeights[elNum] = elWeight;
	localWeightSum += elWeight;
      }
853

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

862
	elInfo = stack.traverseNext(elInfo);
863
864
865
866
867
868
      }
    }

    return localWeightSum;
  }

869

870
  void MeshDistributor::partitionMesh(AdaptInfo *