MeshDistributor.cc 61 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
      mesh->dofCompress();
174
      updateLocalGlobalNumbering();
175
      
176
      // === Update periodic mapping, if there are periodic boundaries. ===
177
      
178
      createPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
179
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
180

181

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

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

186

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

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

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

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

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

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

263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278
      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;
279 280 281 282 283 284 285
    }

    probStat.push_back(probVec);
  }


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

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

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

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

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

312

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

    int nMacroElements = 0;

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

      elInfo = stack.traverseNext(elInfo);
    }

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

334

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

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

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

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

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

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


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

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

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

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

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

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

404

405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423
  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);
      }
    }
  }


424 425 426 427 428 429 430 431 432 433
  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())
434
	    removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455
	}
	
	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
456 457 458 459 460 461 462 463 464 465 466
  {
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
      if (it->second->isPeriodic())
	boundaryMap.erase(it++);
      else
	++it;      
    }    
  }


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

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

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

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

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

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

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

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

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

506

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

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

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

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

526 527

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

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

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

    createPeriodicMap();
542 543 544
  }

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

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

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

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

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

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

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

600
    return meshFitTogether;
601
  }
602 603


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

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

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

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

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

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

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

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

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

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

646 647 648
      meshChanged = true;

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

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

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

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

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

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

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

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

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

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

    return meshChanged;
691 692
  }

693

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

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

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

721 722
    if (!elInfo)
      return value;
723

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

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

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

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

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

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

765 766 767
    return value;
  }

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


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

    int vecSize = 0;
786
    SerUtil::deserialize(in, vecSize);
787
    data.clear();
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
  {