MeshDistributor.cc 66 KB
Newer Older
1
//
2
3
4
5
6
7
8
9
10
11
12
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include <algorithm>
14
15
#include <iostream>
#include <fstream>
16
17
#include <limits>
#include <stdint.h>
Thomas Witkowski's avatar
Thomas Witkowski committed
18
#include <boost/lexical_cast.hpp>
19
20
#include <boost/filesystem.hpp>

21
#include "parallel/MeshDistributor.h"
22
#include "parallel/MeshManipulation.h"
23
#include "parallel/ParallelDebug.h"
24
#include "parallel/StdMpi.h"
25
#include "parallel/MeshPartitioner.h"
26
#include "parallel/ParMetisPartitioner.h"
27
#include "parallel/ZoltanPartitioner.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28
#include "parallel/SimplePartitioner.h"
29
#include "parallel/CheckerPartitioner.h"
30
#include "parallel/MpiHelper.h"
31
#include "parallel/DofComm.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
32
#include "parallel/ParallelProblemStatBase.h"
33
34
#include "io/ElementFileWriter.h"
#include "io/MacroInfo.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
35
#include "io/MacroWriter.h"
36
#include "io/VtkWriter.h"
37
38
39
40
41
#include "Mesh.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "Element.h"
#include "MacroElement.h"
42
43
#include "DOFMatrix.h"
#include "DOFVector.h"
44
#include "SystemVector.h"
45
#include "ElementDofIterator.h"
46
47
#include "ProblemStatBase.h"
#include "StandardProblemIteration.h"
48
#include "VertexVector.h"
49
#include "MeshStructure.h"
50
#include "ProblemStat.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
51
#include "ProblemInstat.h"
52
#include "RefinementManager3d.h"
53
#include "Debug.h"
54

55
56
namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
57
  using boost::lexical_cast;
58
  using namespace boost::filesystem;
59
  using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
60

61
62
  MeshDistributor* MeshDistributor::globalMeshDistributor = NULL;

63
64
65
  const Flag MeshDistributor::BOUNDARY_SUBOBJ_SORTED              = 0X01L;
  const Flag MeshDistributor::BOUNDARY_FILL_INFO_SEND_DOFS        = 0X02L;
  const Flag MeshDistributor::BOUNDARY_FILL_INFO_RECV_DOFS        = 0X04L;
Thomas Witkowski's avatar
Thomas Witkowski committed
66

67
68
69
70
71
  inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2)
  {
    return (*dof1 < *dof2);
  }

72
  MeshDistributor::MeshDistributor()
73
    : problemStat(0),
74
      initialized(false),
75
      name("parallel"),
76
77
78
      mesh(NULL),
      refineManager(NULL),
      partitioner(NULL),
79
      deserialized(false),
80
      writeSerializationFile(false),
81
      repartitioningAllowed(false),
82
      repartitionIthChange(20),
83
84
      nMeshChangesAfterLastRepartitioning(0),
      repartitioningCounter(0),
85
      repartitioningFailed(0),
86
      debugOutputDir(""),
Thomas Witkowski's avatar
Thomas Witkowski committed
87
      lastMeshChangeIndex(0),
88
      createBoundaryDofFlag(0),
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
89
      boundaryDofInfo(1),
90
      meshAdaptivity(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
91
92
      hasPeriodicBoundary(false),
      printTimings(false)
93
  {
94
    FUNCNAME("MeshDistributor::ParalleDomainBase()");
Thomas Witkowski's avatar
Thomas Witkowski committed
95

96
    mpiComm = MPI::COMM_WORLD;
97
98
    mpiRank = mpiComm.Get_rank();
    mpiSize = mpiComm.Get_size();
99

100
    Parameters::get(name + "->repartitioning", repartitioningAllowed);
101
102
    Parameters::get(name + "->debug output dir", debugOutputDir);
    Parameters::get(name + "->repartition ith change", repartitionIthChange);
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
103
    Parameters::get(name + "->mesh adaptivity", meshAdaptivity);
Thomas Witkowski's avatar
Thomas Witkowski committed
104
    
105
    string partStr = "parmetis";
106
    Parameters::get(name + "->partitioner", partStr);
107
108
109
110

    if (partStr == "parmetis") 
      partitioner = new ParMetisPartitioner(&mpiComm);

Thomas Witkowski's avatar
Thomas Witkowski committed
111
112
    if (partStr == "zoltan") {
#ifdef HAVE_ZOLTAN
113
      partitioner = new ZoltanPartitioner(&mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
114
#else
Thomas Witkowski's avatar
Thomas Witkowski committed
115
      ERROR_EXIT("AMDiS was compiled without Zoltan support. Therefore you cannot make use of it!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
116
117
#endif
    }
118

119
120
121
    if (partStr == "checker")
      partitioner = new CheckerPartitioner(&mpiComm);

Thomas Witkowski's avatar
Thomas Witkowski committed
122
123
124
    if (partStr == "simple")
      partitioner = new SimplePartitioner(&mpiComm);

125
    int tmp = 0;
126
    Parameters::get(name + "->box partitioning", tmp);
Thomas Witkowski's avatar
Thomas Witkowski committed
127
128
    partitioner->setBoxPartitioning(static_cast<bool>(tmp));

Thomas Witkowski's avatar
Thomas Witkowski committed
129
130
    Parameters::get(name + "->print timings", printTimings);

131
    TEST_EXIT(partitioner)("Could not create partitioner \"%s\"!\n", partStr.c_str());
132
133
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
134

Thomas Witkowski's avatar
Thomas Witkowski committed
135
136
137
138
139
140
141
  MeshDistributor::~MeshDistributor()
  {
    if (partitioner)
      delete partitioner;
  }


142
  void MeshDistributor::initParallelization()
143
  {
144
    FUNCNAME("MeshDistributor::initParallelization()");
145

146
147
148
    if (initialized)
      return;

Thomas Witkowski's avatar
Thomas Witkowski committed
149
150
151
152
    double first = MPI::Wtime();
    MSG("Initialization phase 1 needed %.5f seconds\n", 
	first - ParallelProblemStatBase::initTimeStamp);

153
154
    TEST_EXIT(mpiSize > 1)
      ("Parallelization does not work with only one process!\n");
155
    TEST_EXIT(feSpaces.size() > 0)
156
      ("No FE space has been defined for the mesh distributor!\n");
157
    TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n");
158

159
160
161
162
163
    // === Sort FE spaces with respect to the degree of the basis ===
    // === functions. Use stuiped bubble sort for this.           ===

    bool doNext = false;
    do {
164
      doNext = false;
165
166
167
168
169
170
      for (unsigned int i = 0; i < feSpaces.size() - 1; i++) {
	if (feSpaces[i]->getBasisFcts()->getDegree() >
	    feSpaces[i + 1]->getBasisFcts()->getDegree()) {
	  const FiniteElemSpace *tmp = feSpaces[i + 1];
	  feSpaces[i + 1] = feSpaces[i];
	  feSpaces[i] = tmp;
171
172
	  doNext = true;
	}
173
      }
174
175
    } while (doNext);

176
    elObjDb.setFeSpace(feSpaces[0]);
177

Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
180
    // If required, create hierarchical mesh level structure.
    createMeshLevelStructure();

181
182
183
    // 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).
184
    if (deserialized) {
Thomas Witkowski's avatar
Thomas Witkowski committed
185
186
      createMeshLevelStructure();
      
187
188
      updateMacroElementInfo();

189
      removePeriodicBoundaryConditions();
190

191
      elObjDb.createMacroElementInfo(allMacroElements);
192

Thomas Witkowski's avatar
Thomas Witkowski committed
193
      updateLocalGlobalNumbering();
194

Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
      elObjDb.setData(partitionMap, levelData);

Thomas Witkowski's avatar
Thomas Witkowski committed
197
#if (DEBUG != 0)
198
      TEST_EXIT_DBG(dofMaps.size())("No DOF mapping defined!\n");
199
      ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], 
200
				    *(dofMaps[0]), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
201
				    debugOutputDir + "mpi-dbg", "dat");
Thomas Witkowski's avatar
Thomas Witkowski committed
202
#endif    
203
204

      initialized = true;
205
      return;
206
    }
207

208

209
210
211
    // 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.
212
213
    testForMacroMesh();

214
215
    // For later mesh repartitioning, we need to store some information about
    // the macro mesh.
216
217
    createMacroElementInfo();

218
    // create an initial partitioning of the mesh
219
    partitioner->createInitialPartitioning();
220

221
    // set the element weights, which are 1 at the very first begin
222
    setInitialElementWeights();
223
224

    // and now partition the mesh    
225
226
    bool partitioningSucceed = partitioner->partition(elemWeights, INITIAL);
    TEST_EXIT(partitioningSucceed)("Initial partitioning does not work!\n");
227
    partitioner->createPartitionMap(partitionMap);
228

229

230
#if (DEBUG != 0)    
231
232
    debug::ElementIdxToDofs elMap;
    debug::createSortedDofs(mesh, elMap);
233
234
#endif

235
    if (mpiRank == 0) {
236
#if (DEBUG != 0)    
237
      int writePartMesh = 1;
238
239
240
#else
      int writePartMesh = 0;
#endif
241
      Parameters::get("dbg->write part mesh", writePartMesh);
242

243
      if (writePartMesh > 0) {
244
	debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex");
245
	ParallelDebug::writePartitioning(*this, debugOutputDir + "part.vtu");
246
      }
247
    }
248

249
    // Create interior boundary information.
250
    createInteriorBoundary(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
251

252
253
    // === Remove neighbourhood relations due to periodic bounday conditions. ===

254
255
256
257
258
    for (deque<MacroElement*>::iterator it = mesh->firstMacroElement();
	 it != mesh->endOfMacroElements(); ++it) {
      for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
	if ((*it)->getNeighbour(i) && 
	    mesh->isPeriodicAssociation((*it)->getBoundary(i))) {
259
260
261

	  int neighIndex = (*it)->getNeighbour(i)->getIndex();

262
263
264
265
266
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
267
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
268
269
270
271
272
273
274
275
276
	}
      }
    }

    for (vector<MacroElement*>::iterator it = allMacroElements.begin();
	 it != allMacroElements.end(); ++it) {
      for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
	if ((*it)->getNeighbour(i) && 
	    mesh->isPeriodicAssociation((*it)->getBoundary(i))) {
277
278
279

	  int neighIndex = (*it)->getNeighbour(i)->getIndex();

280
281
282
283
284
	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
	  (*it)->setNeighbour(i, NULL);
	  (*it)->setBoundary(i, 0);

	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
285
	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
286

287
288
289
	}
      }
    }
290

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
291
292
293
    // === Remove all macro elements that are not part of the rank partition. ===

    removeMacroElements();
294

295
296
    // === Create new global and local DOF numbering. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
297

298
    // We have to remove the VertexVectors, which contain periodic assoiciations, 
299
300
    // because they are not valid anymore after some macro elements have been
    // removed and the corresponding DOFs were deleted.
301
    for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin();
302
303
	 it != mesh->getPeriodicAssociations().end(); ++it)
      const_cast<DOFAdmin&>(mesh->getDofAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second));
304

305
    updateLocalGlobalNumbering();
306

307
308
309
    // === In 3D we have to make some test, if the resulting mesh is valid.  ===
    // === If it is not valid, there is no possiblity yet to fix this        ===
    // === problem, just exit with an error message.                         ===
310

311
312
    check3dValidMesh();

313

314
315
    // === If in debug mode, make some tests. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
316
#if (DEBUG != 0)
317
    MSG("AMDiS runs in debug mode, so make some test ...\n");
318

319
    ParallelDebug::testAllElements(*this);
320
    debug::testSortedDofs(mesh, elMap);
321
    ParallelDebug::testInteriorBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
322
    ParallelDebug::followBoundary(*this);
Thomas Witkowski's avatar
Thomas Witkowski committed
323

324
    debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh");   
325
326

    MSG("Debug mode tests finished!\n");
327
#endif
328

329
330
331
    // Remove periodic boundary conditions in sequential problem definition. 
    removePeriodicBoundaryConditions();

332
333
334
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
335

336
    // === Global refinements. ===
337
    
Thomas Witkowski's avatar
Thomas Witkowski committed
338
    int globalRefinement = 0;
339
    Parameters::get(mesh->getName() + "->global refinements", globalRefinement);
340
    
Thomas Witkowski's avatar
Thomas Witkowski committed
341
    if (globalRefinement > 0) {
342
      refineManager->globalRefine(mesh, globalRefinement);
343

344
      updateLocalGlobalNumbering();
345

346
347
348
#if (DEBUG != 0)
    ParallelDebug::testPeriodicBoundary(*this);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
349
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
350

Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
351
    // And delete some data, we there is no mesh adaptivty.
Thomas Witkowski's avatar
Thomas Witkowski committed
352
    if (!meshAdaptivity)
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
353
354
      elObjDb.clear();

355
    initialized = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
356
    MSG("Init parallelization needed %.5f seconds\n", MPI::Wtime() - first);
357
358
  }

359

360
  void MeshDistributor::addProblemStat(ProblemStatSeq *probStat)
361
  {
362
    FUNCNAME("MeshDistributor::addProblemStat()");
363

364
365
366
    TEST_EXIT_DBG(probStat->getFeSpaces().size())
      ("No FE spaces in stationary problem!\n");

367

368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
    // === Add all FE spaces from stationary problem. ===

    vector<const FiniteElemSpace*> newFeSpaces = probStat->getFeSpaces();
    for (int i = 0; i < static_cast<int>(newFeSpaces.size()); i++)
      if (find(feSpaces.begin(), feSpaces.end(), newFeSpaces[i]) == 
	  feSpaces.end())
	feSpaces.push_back(newFeSpaces[i]);
    

    // === Add mesh of stationary problem and create a corresponding  ===
    // === refinement manager object.                                 ===
    
    if (mesh != NULL) {
      TEST_EXIT(mesh == probStat->getMesh())
	("Does not yet support for different meshes!\n");
    } else {
      mesh = probStat->getMesh();
    }
386
387
388
389
390
391
392
393
394
395
    
    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());
396
    }
397
398
399
    
    partitioner->setMesh(mesh);
    
400

401
402
403
    // === Check whether the stationary problem should be serialized. ===


404
405
    // Create parallel serialization file writer, if needed.
    int writeSerialization = 0;
406
    Parameters::get(probStat->getName() + "->output->write serialization",  
407
		    writeSerialization);
408
    if (writeSerialization && !writeSerializationFile) {
409
      string filename = "";
410
      Parameters::get(name + "->output->serialization filename", filename);
411
412
413
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel serialization file!\n");
414
415

      int tsModulo = -1;
416
      Parameters::get(probStat->getName() + "->output->write every i-th timestep", 
417
		      tsModulo);
418
      
419
      probStat->getFileWriterList().push_back(new Serializer<MeshDistributor>(this, filename, tsModulo));
420
421
      writeSerializationFile = true;
    }    
422

423
424
425

    // === Check whether the stationary problem should be deserialized. ===

426
    int readSerialization = 0;
427
    Parameters::get(probStat->getName() + "->input->read serialization", 
428
		    readSerialization);
429
    if (readSerialization) {
430
      string filename = "";
431
      Parameters::get(probStat->getName() + "->input->serialization filename", 
432
		      filename);
433
      filename += ".p" + lexical_cast<string>(mpiRank);
434
      MSG("Start deserialization with %s\n", filename.c_str());
435
      ifstream in(filename.c_str());
436
437
438
439

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

440
441
442
443
444
      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      int revNumber = -1;
      SerUtil::deserialize(in, revNumber);

445
      probStat->deserialize(in);
446
      in.close();
447
448
      MSG("Deserialization from file: %s\n", filename.c_str());

449
      filename = "";
450
      Parameters::get(name + "->input->serialization filename", filename);
451
452
453
454
      
      TEST_EXIT(filename != "")
	("No filename defined for parallel deserialization file!\n");
      
455
      string rankFilename = filename + ".p" + lexical_cast<string>(mpiRank);
456
457
458
459
      in.open(rankFilename.c_str());
      
      TEST_EXIT(!in.fail())("Could not open parallel deserialization file: %s\n",
			    filename.c_str());
460
461
462
463
464

      // Read the revision number of the AMDiS version which was used to create 
      // the serialization file.
      revNumber = -1;
      SerUtil::deserialize(in, revNumber);
465
466
467
468
469
      
      deserialize(in);
      in.close();
      MSG("Deserializtion of mesh distributor from file: %s\n", rankFilename.c_str());
      deserialized = true;
470
471
    }

472
    problemStat.push_back(probStat);
473

474
475
476

    // === If the mesh distributor is already initialized, don't forget to  ===
    // === remove the periodic boundary conditions on these objects.        ===
477

478
    if (initialized)
479
      removePeriodicBoundaryConditions(probStat);
480
481
482
  }


483
484
485
486
487
488
489
490
491
492
493
  void MeshDistributor::addProblemStatGlobal(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::addProblemStatGlobal()");

    if (globalMeshDistributor == NULL)
      globalMeshDistributor = new MeshDistributor();

    globalMeshDistributor->addProblemStat(probStat);
  }


494
  void MeshDistributor::exitParallelization()
495
  {}
496

497
498
499
500
501
502
503
504
505
506
507
508

  void MeshDistributor::registerDofMap(ParallelDofMapping &dofMap)
  {
    FUNCNAME("MeshDistributor::registerDofMap()");

    TEST_EXIT(find(dofMaps.begin(), dofMaps.end(), &dofMap) ==
	      dofMaps.end())
      ("Parallel DOF mapping already registerd in mesh distributor object!\n");

    dofMaps.push_back(&dofMap);
  }

509
  
510
  void MeshDistributor::testForMacroMesh()
511
  {
512
    FUNCNAME("MeshDistributor::testForMacroMesh()");
513
514
515
516
517
518
519

    int nMacroElements = 0;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      TEST_EXIT(elInfo->getLevel() == 0)
520
	("Mesh is already refined! This does not work with parallelization!\n");
521
522
523

      TEST_EXIT(elInfo->getType() == 0)
	("Only macro elements with level 0 are supported!\n");
524
525
526
527
528
529
530
531
532
533
      
      nMacroElements++;

      elInfo = stack.traverseNext(elInfo);
    }

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

534

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
535
  void MeshDistributor::synchVector(SystemVector &vec, int level)
Thomas Witkowski's avatar
Thomas Witkowski committed
536
  {
537
538
    FUNCNAME("MeshDistributor::synchVector()");

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
539
    TEST_EXIT_DBG(level >= 0 && level <= 1)("Wrong level number!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
540

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
541
542
543
544
545
546
    MPI::Intracomm &levelComm = levelData.getMpiComm(level);
    DofComm &dc = (level == 0 ? dofComm : dofCommSd);

    StdMpi<vector<double> > stdMpi(levelComm);

    for (DofComm::Iterator it(dc.getSendDofs()); 
547
	 !it.end(); it.nextRank()) {
548
      vector<double> dofs;
549

550
551
      for (int i = 0; i < vec.getSize(); i++) {
	DOFVector<double> &dofVec = *(vec.getDOFVector(i));
552

Thomas Witkowski's avatar
Thomas Witkowski committed
553
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof())
554
555
	  dofs.push_back(dofVec[it.getDofIndex()]);
      }
556

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
557
558
559
560
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.send(rank, dofs);
561
    }
562
	   
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
563
564
565
566
567
568
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
      stdMpi.recv(rank);
    }
569

570
    stdMpi.startCommunication();
571

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
572
573
574
575
    for (DofComm::Iterator it(dc.getRecvDofs()); !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);
576
577

      int counter = 0;
578
      vector<double> &recvData =  stdMpi.getRecvData(rank);
579
580
581
582

      for (int i = 0; i < vec.getSize(); i++) {
	DOFVector<double> &dofVec = *(vec.getDOFVector(i));

583
584
585
586
587
	for (it.beginDofIter(vec.getFeSpace(i)); !it.endDofIter(); it.nextDof()) {
	  TEST_EXIT_DBG(counter < recvData.size())
	    ("Recv data from rank %d has only %d entries!\n", rank, recvData.size());
	  dofVec[it.getDofIndex()] = recvData[counter++];
	}
588
589
590
591
      }
    }
  }

592

593
594
595
596
597
598
599
  void MeshDistributor::check3dValidMesh()
  {
    FUNCNAME("MeshDistributor::check3dValidMesh()");

    if (mesh->getDim() != 3)
      return;

600
601
602
    // === Create set of all edges and all macro element indices in rank's ===
    // === subdomain.                                                      ===

603
604
605
606
607
608
    std::set<DofEdge> allEdges;
    std::set<int> rankMacroEls;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
609
      for (int i = 0; i < 6; i++) {
610
	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
611
612
	DofEdge e = elObjDb.getEdgeLocalMap(elData);
	allEdges.insert(e);
613
614
615
616
617
618
619
      }

      rankMacroEls.insert(elInfo->getElement()->getIndex());

      elInfo = stack.traverseNext(elInfo);
    }

620

621
622
623
    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
    // === for the opposite.                                                 ===

624
    FixRefinementPatch::connectedEdges.clear();
625
    bool valid3dMesh = true;
626

627
    // Traverse all edges
628
    for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
629
      // For each edge get all macro elements that contain this edge.
630
      vector<ElementObjectData>& edgeEls = elObjDb.getElements(*it);
631
632
633
634

      TEST_EXIT_DBG(edgeEls.size() > 0)
	("No edge %d/%d in elObjDB!\n", it->first, it->second);

635
636
637
638
639
640
641

      // We have now to check, if the current edge connects two macro elements,
      // which are otherwise not connected. The basic idea to check this is very
      // simple: We take the first macro element in rank's subdomain that contain
      // this edge and add it to the set variable "el0". All other macro elements
      // which share this edge are added to "el1". 

642
      std::set<int> el0, el1;
643
      map<int, int> edgeNoInEl;
644
645
646

      for (unsigned int i = 0; i < edgeEls.size(); i++) {
	if (rankMacroEls.count(edgeEls[i].elIndex)) {
647
648
649
650
	  if (el0.empty())
	    el0.insert(edgeEls[i].elIndex);
	  else
	    el1.insert(edgeEls[i].elIndex);
651
652
653
654

	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
	}
      }
655
           
656
657
      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");

658
659
      // If el1 is empty, there is only on macro element in the mesh which
      // contains this edge. Hence, we can continue with checking another edge.
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
      if (el1.empty())
	continue;
      
      bool found = false;      
      do {
	found = false;
	for (std::set<int>::iterator elIt = el0.begin(); elIt != el0.end(); ++elIt) {
	  for (int i = 0; i < 4; i++) {
	    int neighEl = macroElementNeighbours[*elIt][i];
	    if (neighEl != -1 && el1.count(neighEl)) {
	      el0.insert(neighEl);
	      el1.erase(neighEl);
	      found = true;
	    }
	  }
	  
	  if (found)
	    break;
	}
      } while (found);
      
      if (!el1.empty()) {
682
683
684
685
	valid3dMesh = false;

	MSG_DBG("Unvalid 3D mesh with %d (%d - %d) elements on this edge!\n", 
		edgeEls.size(), el0.size(), el1.size());
686

687
688
689
	for (std::set<int>::iterator elIt0 = el0.begin(); elIt0 != el0.end(); ++elIt0) {
	  for (std::set<int>::iterator elIt1 = el1.begin(); elIt1 != el1.end(); ++elIt1) {
	    pair<Element*, int> edge0 = 
690
	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
691
	    pair<Element*, int> edge1 = 
692
	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
693

694
#if (DEBUG != 0)
695
696
697
698
	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
	    DofEdge dofEdge1 = edge1.first->getEdge(edge1.second);

	    WorldVector<double> c0, c1;
699
700
	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
701

702
703
704
705
706
707
708
709
710
	    MSG("FOUND EDGE %d/%d <-> %d/%d\n", 
		edge0.first->getIndex(), edge0.second,
		edge1.first->getIndex(), edge1.second);

	    MSG("FOUND EDGE: %d %d with coords %f %f %f   %f %f %f\n",
		dofEdge0.first, dofEdge0.second, 
		c0[0], c0[1], c0[2],
		c1[0], c1[1], c1[2]);
#endif
711
712
713
714
715
716
717
718
719

	    TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n");

	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge0, edge1));
	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
	  }
	}
      }
    }
720
721

    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
722
723
724
  }


725
  void MeshDistributor::getAllBoundaryDofs(const FiniteElemSpace *feSpace,
726
					   int level,
727
					   DofContainer& dofs)
728
729
730
731
  {
    FUNCNAME("MeshDistributor::getAllBoundaryDofs()");

    DofContainerSet dofSet;
732
733
    for (DofComm::Iterator it(dofComm.getSendDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
734
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
735
736
    for (DofComm::Iterator it(dofComm.getRecvDofs(), level, feSpace); 
	 !it.end(); it.nextRank())
737
      dofSet.insert(it.getDofs().begin(), it.getDofs().end());
738
739
740
741
742
743

    dofs.clear();
    dofs.insert(dofs.begin(), dofSet.begin(), dofSet.end());
  }


744
745
  void MeshDistributor::removePeriodicBoundaryConditions()
  {
746
747
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

748
    // Remove periodic boundaries in boundary manager on matrices and vectors.
749
750
    for (unsigned int i = 0; i < problemStat.size(); i++)
      removePeriodicBoundaryConditions(problemStat[i]);
751
752
753

    // Remove periodic boundaries on elements in mesh.
    TraverseStack stack;
754
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
755
756
757
758
    while (elInfo) {
      elInfo->getElement()->deleteElementData(PERIODIC);
      elInfo = stack.traverseNext(elInfo);
    }    
759
760
761

    // Remove periodic vertex associations
    mesh->getPeriodicAssociations().clear();
762
763
764
  }


765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
  void MeshDistributor::removePeriodicBoundaryConditions(ProblemStatSeq *probStat)
  {
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

    int nComponents = probStat->getNumComponents();

    for (int j = 0; j < nComponents; j++) {
      for (int k = 0; k < nComponents; k++) {
	DOFMatrix* mat = probStat->getSystemMatrix(j, k);
	if (mat && mat->getBoundaryManager())
	  removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(mat->getBoundaryManager()->getBoundaryConditionMap()));
      }
      
      if (probStat->getSolution()->getDOFVector(j)->getBoundaryManager())
	removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat->getSolution()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
      
      if (probStat->getRhs()->getDOFVector(j)->getBoundaryManager())
	removePeriodicBoundaryConditions(const_cast<BoundaryIndexMap&>(probStat->getRhs()->getDOFVector(j)->getBoundaryManager()->getBoundaryConditionMap()));
    }
  }


787
  void MeshDistributor::removePeriodicBoundaryConditions(BoundaryIndexMap& boundaryMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
788
  {
789
790
    FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()");

Thomas Witkowski's avatar
Thomas Witkowski committed
791
792
    BoundaryIndexMap::iterator it = boundaryMap.begin();
    while (it != boundaryMap.end()) {
793
      if (it->second->isPeriodic())
Thomas Witkowski's avatar
Thomas Witkowski committed
794
	boundaryMap.erase(it++);
795
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
796
797
798
799
800
	++it;      
    }    
  }


801
802
803
804
  void MeshDistributor::createMeshLevelStructure()
  {
    FUNCNAME("MeshDistributor::createMeshLevelStructure()");

805
806
807
    if (mpiSize != 16)
      return;

808
    std::set<int> neighbours;
809
810
    switch (mpiRank) {
    case 0:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
811
      neighbours.insert(1); neighbours.insert(2); neighbours.insert(3);
812
813
      break;
    case 1:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
814
      neighbours.insert(0); neighbours.insert(4); neighbours.insert(2); neighbours.insert(3); neighbours.insert(6);
815
816
      break;
    case 2:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
817
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(3); neighbours.insert(8); neighbours.insert(9);
818
819
      break;
    case 3:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
820
821
822
      neighbours.insert(0); neighbours.insert(1); neighbours.insert(4); 
      neighbours.insert(2); neighbours.insert(6);
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12); 
823
824
      break;
    case 4:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
825
      neighbours.insert(1); neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); neighbours.insert(5);
826
827
      break;
    case 5:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
828
      neighbours.insert(4); neighbours.insert(6); neighbours.insert(7);
829
830
      break;
    case 6:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
831
832
833
      neighbours.insert(1); neighbours.insert(4); neighbours.insert(5); 
      neighbours.insert(3); neighbours.insert(7);
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13); 
834
835
      break;
    case 7:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
836
      neighbours.insert(4); neighbours.insert(5); neighbours.insert(6);  neighbours.insert(12); neighbours.insert(13);
837
838
      break;
    case 8:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
839
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(9);  neighbours.insert(10); neighbours.insert(11);
840
841
      break;
    case 9:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
842
843
844
      neighbours.insert(2); neighbours.insert(3); neighbours.insert(6); 
      neighbours.insert(8); neighbours.insert(12);
      neighbours.insert(10); neighbours.insert(11); neighbours.insert(14); 
845
846
      break;
    case 10:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
847
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(11);
848
849
      break;
    case 11:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
850
      neighbours.insert(8); neighbours.insert(9); neighbours.insert(12);  neighbours.insert(10); neighbours.insert(14);
851
852
      break;
    case 12:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
853
854
855
      neighbours.insert(3); neighbours.insert(6); neighbours.insert(7); 
      neighbours.insert(9); neighbours.insert(13);
      neighbours.insert(11); neighbours.insert(14); neighbours.insert(15); 
856
857
      break;
    case 13:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
858
      neighbours.insert(6); neighbours.insert(7); neighbours.insert(12);  neighbours.insert(14); neighbours.insert(15);
859
860
      break;
    case 14:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
861
      neighbours.insert(9); neighbours.insert(12); neighbours.insert(13);  neighbours.insert(11); neighbours.insert(15);
862
863
      break;
    case 15:
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
864
      neighbours.insert(12); neighbours.insert(13); neighbours.insert(14);
865
866
      break;
    }    
867

868
    TEST_EXIT(neighbours.size() > 0)("Should not happen!\n");
869
870
871
872
873
874
875

    levelData.init(neighbours);

    bool multiLevelTest = false;
    Parameters::get("parallel->multi level test", multiLevelTest);
    if (multiLevelTest) {
      switch (mpiRank) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
876
      case 0: case 1: case 2: case 3:
877
878
	{
	  std::set<int> m;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
879
	  m.insert(0); m.insert(1); m.insert(2); m.insert(3);
880
	  levelData.addLevel(m, 0);