InteriorBoundary.cc 22.2 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.


13
#include "parallel/InteriorBoundary.h"
14
#include "parallel/ElementObjectDatabase.h"
15
#include "parallel/MeshDistributor.h"
16
17
#include "FiniteElemSpace.h"
#include "BasisFunction.h"
18
#include "Serializer.h"
19
#include "VertexVector.h"
20
21

namespace AMDiS {
Thomas Witkowski's avatar
Thomas Witkowski committed
22

23
24
25
  using namespace std;


26
27
  void InteriorBoundary::create(MeshLevelData &levelData,
				int level,
28
29
30
31
32
33
34
35
36
37
				ElementObjectDatabase &elObjDb)
  { 
    FUNCNAME("InteriorBoundary::clear()");

    own.clear();
    other.clear();
    periodic.clear();

    Mesh *mesh = elObjDb.getMesh();
    TEST_EXIT_DBG(mesh)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
38
    TEST_EXIT_DBG(level < levelData.getLevelNumber())("Should not happen!\n");
39

Thomas Witkowski's avatar
Thomas Witkowski committed
40
41
42
    MPI::Intracomm mpiComm = levelData.getMpiComm(level);
    int levelMpiRank = mpiComm.Get_rank();
    int globalMpiRank = MPI::COMM_WORLD.Get_rank();
43
    std::set<int> levelRanks = levelData.getLevelRanks(level);
44
    bool allRanks = (levelRanks.size() == 1 && *(levelRanks.begin()) == -1);
45

46
47
48
49
    // === Create interior boundary data structure. ===
    
    for (int geoPos = 0; geoPos < mesh->getDim(); geoPos++) {
      GeoIndex geoIndex = INDEX_OF_DIM(geoPos, mesh->getDim());
50

51
52
      while (elObjDb.iterate(geoIndex)) {
	map<int, ElementObjectData>& objData = elObjDb.getIterateData();
53
54

	// Test, if this is a boundary object of this rank.
Thomas Witkowski's avatar
Thomas Witkowski committed
55
	if (!(objData.count(globalMpiRank) && objData.size() > 1))
56
57
	  continue;

58
59
60
	// Test, if the boundary object defines an interior boundary within the
	// ranks of the MPI group. If not, go to next element.
	bool boundaryWithinMpiGroup = false;
61
	if (allRanks) {
62
63
64
65
	  boundaryWithinMpiGroup = true;
	} else {
	  for (map<int, ElementObjectData>::iterator it = objData.begin();
	       it != objData.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
66
	    if (it->first != globalMpiRank && levelRanks.count(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
67
	      boundaryWithinMpiGroup = true;
68
69
70
71
72
73
74
75
	      break;
	    }
	  }
	}
	if (!boundaryWithinMpiGroup)
	  continue;

	int owner = elObjDb.getIterateOwner(level);
Thomas Witkowski's avatar
Thomas Witkowski committed
76
	ElementObjectData& rankBoundEl = objData[globalMpiRank];
77
78
79
80
81
82
83
84
	
	AtomicBoundary bound;
	bound.maxLevel = elObjDb.getIterateMaxLevel();
	bound.rankObj.el = elObjDb.getElementPtr(rankBoundEl.elIndex);
	bound.rankObj.elIndex = rankBoundEl.elIndex;
	bound.rankObj.elType = elObjDb.getElementType(rankBoundEl.elIndex);
	bound.rankObj.subObj = geoIndex;
	bound.rankObj.ithObj = rankBoundEl.ithObject;
Thomas Witkowski's avatar
Thomas Witkowski committed
85
86


87
88
89
90
91
92
93
94
	if (geoIndex == FACE) {
	  for (int edgeNo = 0; edgeNo < 3; edgeNo++) {
	    int edgeOfFace = 
	      bound.rankObj.el->getEdgeOfFace(bound.rankObj.ithObj, edgeNo);
	    
	    bound.rankObj.excludedSubstructures.push_back(make_pair(EDGE, edgeOfFace));
	  }
	}
95
	
Thomas Witkowski's avatar
Thomas Witkowski committed
96

Thomas Witkowski's avatar
Thomas Witkowski committed
97
	if (owner == globalMpiRank) {
98
99
	  for (map<int, ElementObjectData>::iterator it2 = objData.begin();
	       it2 != objData.end(); ++it2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
100
	    if (it2->first == globalMpiRank)
101
	      continue;
102
103
104
105

	    if (!allRanks && levelRanks.count(it2->first) == 0)
	      continue;

106
107
108
109
110
111
112
113
	    bound.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
	    bound.neighObj.elIndex = it2->second.elIndex;
	    bound.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
	    bound.neighObj.subObj = geoIndex;
	    bound.neighObj.ithObj = it2->second.ithObject;
	    
	    bound.type = INTERIOR;
	    
114
	    AtomicBoundary& b = getNewOwn(it2->first);
115
116
117
118
119
120
121
122
123
124
125
126
	    b = bound;
	    if (geoIndex == EDGE)
	      b.neighObj.reverseMode = 
		elObjDb.getEdgeReverseMode(rankBoundEl, it2->second);
	    if (geoIndex == FACE)
	      b.neighObj.reverseMode = 
		elObjDb.getFaceReverseMode(rankBoundEl, it2->second);
	  }
	  
	} else {
	  TEST_EXIT_DBG(objData.count(owner) == 1)
	    ("Should not happen!\n");
127
	  
128
	  ElementObjectData& ownerBoundEl = objData[owner];
129
  
130
131
132
133
134
135
136
137
	  bound.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
	  bound.neighObj.elIndex = ownerBoundEl.elIndex;
	  bound.neighObj.elType = -1;
	  bound.neighObj.subObj = geoIndex;
	  bound.neighObj.ithObj = ownerBoundEl.ithObject;
	  
	  bound.type = INTERIOR;
	  
138
	  AtomicBoundary& b = getNewOther(owner);
139
140
141
142
143
144
145
146
	  b = bound;	    
	  if (geoIndex == EDGE)
	    b.rankObj.reverseMode =
	      elObjDb.getEdgeReverseMode(rankBoundEl, ownerBoundEl);
	  if (geoIndex == FACE)
	    b.rankObj.reverseMode = 
	      elObjDb.getFaceReverseMode(rankBoundEl, ownerBoundEl);
	}
147
      }
148
149
150
    }


151
    // === Create periodic boundary data structure. ===
152

153
154
155
156
    int removePeriodicBoundary = 0;
    Parameters::get("parallel->remove periodic boundary", removePeriodicBoundary);


157
158
    for (PerBoundMap<DegreeOfFreedom>::iterator it = elObjDb.getPeriodicVertices().begin();
	 it != elObjDb.getPeriodicVertices().end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
159
      if (elObjDb.isInRank(it->first.first, globalMpiRank) == false)
160
161
162
	continue;

      ElementObjectData& perDofEl0 = 
Thomas Witkowski's avatar
Thomas Witkowski committed
163
	elObjDb.getElementsInRank(it->first.first)[globalMpiRank];
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183

      for (map<int, ElementObjectData>::iterator elIt = elObjDb.getElementsInRank(it->first.second).begin();
	   elIt != elObjDb.getElementsInRank(it->first.second).end(); ++elIt) {

	int otherElementRank = elIt->first;
	ElementObjectData& perDofEl1 = elIt->second;

	AtomicBoundary bound;
	bound.rankObj.el = elObjDb.getElementPtr(perDofEl0.elIndex);
	bound.rankObj.elIndex = perDofEl0.elIndex;
	bound.rankObj.elType = elObjDb.getElementType(perDofEl0.elIndex);
	bound.rankObj.subObj = VERTEX;
	bound.rankObj.ithObj = perDofEl0.ithObject;

	bound.neighObj.el = elObjDb.getElementPtr(perDofEl1.elIndex);
	bound.neighObj.elIndex = perDofEl1.elIndex;
	bound.neighObj.elType = elObjDb.getElementType(perDofEl1.elIndex);
	bound.neighObj.subObj = VERTEX;
	bound.neighObj.ithObj = perDofEl1.ithObject;

184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
	if (removePeriodicBoundary) {
	  bound.type = INTERIOR;
	  
	  map<int, ElementObjectData> objData = 
	    elObjDb.getElementsInRank(it->first.first);
	  objData.insert(elObjDb.getElementsInRank(it->first.second).begin(),
			 elObjDb.getElementsInRank(it->first.second).end());
	  
	  int owner = std::max(elObjDb.getOwner(it->first.first, level),
			       elObjDb.getOwner(it->first.second, level));
	  if (owner == globalMpiRank) {
	    for (map<int, ElementObjectData>::iterator it2 = objData.begin();
	     	 it2 != objData.end(); ++it2) {
	      if (it2->first == globalMpiRank)
	     	continue;
	      
	      if (!allRanks && levelRanks.count(it2->first) == 0)
	     	continue;
	      
	      AtomicBoundary& b = getNewOwn(it2->first);	  
	      b = bound;	  
	      b.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
	      b.neighObj.elIndex = it2->second.elIndex;
	      b.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
	      b.neighObj.ithObj = it2->second.ithObject;
	    }
	  } else {
	    AtomicBoundary& b = getNewOther(owner);
	    b = bound;
	    
	    ElementObjectData& ownerBoundEl = objData[owner];    
	    b.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
	    b.neighObj.elIndex = ownerBoundEl.elIndex;
	    b.neighObj.elType = -1;
	    b.neighObj.ithObj = ownerBoundEl.ithObject;
219
220
	  }
	} else {
221
222
223
224
225
226
227
228
229
230
231
232
233
234
	  bound.type = it->second;
	  
	  if (MeshDistributor::sebastianMode) {
	    if (bound.rankObj.elIndex == 3 &&
		bound.rankObj.ithObj == 1 ||
		bound.rankObj.elIndex == 78 &&
		bound.rankObj.ithObj == 2) {
	      AtomicBoundary& b = getNewPeriodic(otherElementRank);
	      b = bound;	    	    
	    }
	  } else {
	    AtomicBoundary& b = getNewPeriodic(otherElementRank);
	    b = bound;	    	  
	  }
235
	}
236
237
238
239
240
241
      }
    }


    for (PerBoundMap<DofEdge>::iterator it = elObjDb.getPeriodicEdges().begin();
	 it != elObjDb.getPeriodicEdges().end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
242
      if (elObjDb.isInRank(it->first.first, globalMpiRank) == false)
243
244
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
      ElementObjectData& perEdgeEl0 = 
	elObjDb.getElementsInRank(it->first.first)[globalMpiRank];
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265

      for (map<int, ElementObjectData>::iterator elIt = elObjDb.getElementsInRank(it->first.second).begin();
 	   elIt != elObjDb.getElementsInRank(it->first.second).end(); ++elIt) {
      
	int otherElementRank = elIt->first;
	ElementObjectData& perEdgeEl1 = elIt->second;

	AtomicBoundary bound;	    	    
	bound.rankObj.el = elObjDb.getElementPtr(perEdgeEl0.elIndex);
	bound.rankObj.elIndex = perEdgeEl0.elIndex;
	bound.rankObj.elType = elObjDb.getElementType(perEdgeEl0.elIndex);
	bound.rankObj.subObj = EDGE;
	bound.rankObj.ithObj = perEdgeEl0.ithObject;
	
	bound.neighObj.el = elObjDb.getElementPtr(perEdgeEl1.elIndex);
	bound.neighObj.elIndex = perEdgeEl1.elIndex;
	bound.neighObj.elType = elObjDb.getElementType(perEdgeEl1.elIndex);
	bound.neighObj.subObj = EDGE;
	bound.neighObj.ithObj = perEdgeEl1.ithObject;
266

267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
	if (removePeriodicBoundary) {
	  bound.type = INTERIOR;

	  map<int, ElementObjectData> objData = 
	    elObjDb.getElementsInRank(it->first.first);
	  objData.insert(elObjDb.getElementsInRank(it->first.second).begin(),
			 elObjDb.getElementsInRank(it->first.second).end());
	  
	  int owner = std::max(elObjDb.getOwner(it->first.first, level),
			       elObjDb.getOwner(it->first.second, level));
	  ElementObjectData& rankBoundEl = objData[globalMpiRank];

	  if (owner == globalMpiRank) {
	    for (map<int, ElementObjectData>::iterator it2 = objData.begin();
	     	 it2 != objData.end(); ++it2) {
	      if (it2->first == globalMpiRank)
	     	continue;
	      
	      if (!allRanks && levelRanks.count(it2->first) == 0)
	     	continue;
	      
	      AtomicBoundary& b = getNewOwn(it2->first);	  
	      b = bound;	  
	      b.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
	      b.neighObj.elIndex = it2->second.elIndex;
	      b.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
	      b.neighObj.ithObj = it2->second.ithObject;
	      b.neighObj.reverseMode = 
		elObjDb.getEdgeReverseMode(rankBoundEl, it2->second);
	    }
	  } else {
	    AtomicBoundary& b = getNewOther(owner);
	    b = bound;
	    
	    ElementObjectData& ownerBoundEl = objData[owner];    
	    b.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
	    b.neighObj.elIndex = ownerBoundEl.elIndex;
	    b.neighObj.elType = -1;
	    b.neighObj.ithObj = ownerBoundEl.ithObject;
	    b.rankObj.reverseMode = 
	      elObjDb.getEdgeReverseMode(rankBoundEl, ownerBoundEl);
	  }
	} else {
	  bound.type = it->second;
	  
	  AtomicBoundary& b = getNewPeriodic(otherElementRank);
	  b = bound;
	  
	  if (globalMpiRank > otherElementRank)
	    b.neighObj.reverseMode = 
	      elObjDb.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
	  else
	    b.rankObj.reverseMode = 
	      elObjDb.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
	}
322
323
324
325
326
327
      }
    }


    for (PerBoundMap<DofFace>::iterator it = elObjDb.getPeriodicFaces().begin();
	 it != elObjDb.getPeriodicFaces().end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
328
      if (elObjDb.isInRank(it->first.first, globalMpiRank) == false)
329
330
331
332
333
334
335
	continue;

      TEST_EXIT_DBG(elObjDb.getElements(it->first.first).size() == 1)
 	("Should not happen!\n");
      TEST_EXIT_DBG(elObjDb.getElements(it->first.second).size() == 1)
 	("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
      ElementObjectData& perFaceEl0 = 
	elObjDb.getElementsInRank(it->first.first)[globalMpiRank];
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357

      for (map<int, ElementObjectData>::iterator elIt = elObjDb.getElementsInRank(it->first.second).begin();
 	   elIt != elObjDb.getElementsInRank(it->first.second).end(); ++elIt) {
      
	int otherElementRank = elIt->first;
	ElementObjectData& perFaceEl1 = elIt->second;

	AtomicBoundary bound;	    	    
	bound.rankObj.el = elObjDb.getElementPtr(perFaceEl0.elIndex);
	bound.rankObj.elIndex = perFaceEl0.elIndex;
	bound.rankObj.elType = elObjDb.getElementType(perFaceEl0.elIndex);
	bound.rankObj.subObj = FACE;
	bound.rankObj.ithObj = perFaceEl0.ithObject;
	
	bound.neighObj.el = elObjDb.getElementPtr(perFaceEl1.elIndex);
	bound.neighObj.elIndex = perFaceEl1.elIndex;
	bound.neighObj.elType = elObjDb.getElementType(perFaceEl1.elIndex);
	bound.neighObj.subObj = FACE;
	bound.neighObj.ithObj = perFaceEl1.ithObject;
	
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
	if (removePeriodicBoundary) {
	  bound.type = INTERIOR;

	  map<int, ElementObjectData> objData = 
	    elObjDb.getElementsInRank(it->first.first);
	  objData.insert(elObjDb.getElementsInRank(it->first.second).begin(),
			 elObjDb.getElementsInRank(it->first.second).end());
	  
	  int owner = std::max(elObjDb.getOwner(it->first.first, level),
			       elObjDb.getOwner(it->first.second, level));
	  ElementObjectData& rankBoundEl = objData[globalMpiRank];

	  if (owner == globalMpiRank) {
	    for (map<int, ElementObjectData>::iterator it2 = objData.begin();
	     	 it2 != objData.end(); ++it2) {
	      if (it2->first == globalMpiRank)
	     	continue;
	      
	      if (!allRanks && levelRanks.count(it2->first) == 0)
	     	continue;
	      
	      AtomicBoundary& b = getNewOwn(it2->first);	  
	      b = bound;	  
	      b.neighObj.el = elObjDb.getElementPtr(it2->second.elIndex);
	      b.neighObj.elIndex = it2->second.elIndex;
	      b.neighObj.elType = elObjDb.getElementType(it2->second.elIndex);
	      b.neighObj.ithObj = it2->second.ithObject;
	      b.neighObj.reverseMode = 
		elObjDb.getFaceReverseMode(rankBoundEl, it2->second);
	    }
	  } else {
	    AtomicBoundary& b = getNewOther(owner);
	    b = bound;
	    
	    ElementObjectData& ownerBoundEl = objData[owner];    
	    b.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
	    b.neighObj.elIndex = ownerBoundEl.elIndex;
	    b.neighObj.elType = -1;
	    b.neighObj.ithObj = ownerBoundEl.ithObject;
	    b.rankObj.reverseMode = 
	      elObjDb.getFaceReverseMode(rankBoundEl, ownerBoundEl);
	  }
	} else {
	  bound.type = it->second;
	  
	  AtomicBoundary& b = getNewPeriodic(otherElementRank);
	  b = bound;
	  
	  if (globalMpiRank > otherElementRank)
	    b.neighObj.reverseMode = 
	      elObjDb.getFaceReverseMode(perFaceEl0, perFaceEl1);
	  else
	    b.rankObj.reverseMode = 
	      elObjDb.getFaceReverseMode(perFaceEl0, perFaceEl1);
	}
413
414
415
416
417
418
419
420
421
      }
    }

    // === Once we have this information, we must care about the order of the ===
    // === atomic bounds in the three boundary handling object. Eventually    ===
    // === all the boundaries have to be in the same order on both ranks that ===
    // === share the bounday.                                                 ===

    StdMpi<vector<AtomicBoundary> > stdMpi(mpiComm);
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
424
425
426
    if (level == 0) {
      stdMpi.send(own);
      stdMpi.recv(other);
    } else {
      for (RankToBoundMap::iterator rankIt = own.begin();	 
427
428
	   rankIt != own.end(); ++rankIt)
	stdMpi.send(levelData.mapRank(rankIt->first, 0, level), rankIt->second);
Thomas Witkowski's avatar
Thomas Witkowski committed
429
      for (RankToBoundMap::iterator rankIt = other.begin();
430
431
	   rankIt != other.end(); ++rankIt)
	stdMpi.recv(levelData.mapRank(rankIt->first, 0, level));
Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
    }

434
435
436
437
438
439
440
441
442
443
444
    stdMpi.startCommunication();


    // === The information about all neighbouring boundaries has been         ===
    // === received. So the rank tests if its own atomic boundaries are in    ===
    // === the same order. If not, the atomic boundaries are swaped to the    ===
    // === correct order.                                                     ===

    for (RankToBoundMap::iterator rankIt = other.begin();
	 rankIt != other.end(); ++rankIt) {

445
446
447
448
449
450
451
      int rank = rankIt->first;
      if (level > 0)
	rank = levelData.mapRank(rank, 0, level);

      // === We have received from "rank" the ordered list of element       ===
      // === indices. Now, we have to sort the corresponding list in this   ===
      // === rank to get the same order.                                    ===
452
453
454
455
456
     
      for (unsigned int j = 0; j < rankIt->second.size(); j++) {

	// If the expected object is not at place, search for it.

Thomas Witkowski's avatar
Thomas Witkowski committed
457
	BoundaryObject &recvedBound = 
458
	  stdMpi.getRecvData()[rank][j].rankObj;
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488

	if ((rankIt->second)[j].neighObj != recvedBound) {
	  unsigned int k = j + 1;

	  for (; k < rankIt->second.size(); k++)
 	    if ((rankIt->second)[k].neighObj == recvedBound)
	      break;

	  // The element must always be found, because the list is just in
	  // another order.
	  TEST_EXIT_DBG(k < rankIt->second.size())("Should never happen!\n");

	  // Swap the current with the found element.
	  AtomicBoundary tmpBound = (rankIt->second)[k];
	  (rankIt->second)[k] = (rankIt->second)[j];
	  (rankIt->second)[j] = tmpBound;	
	}
      }
    }


    // === Do the same for the periodic boundaries. ===

    if (periodic.size() > 0) {
      stdMpi.clear();

      RankToBoundMap sendBounds, recvBounds;
      for (RankToBoundMap::iterator rankIt = periodic.begin();
	   rankIt != periodic.end(); ++rankIt) {

Thomas Witkowski's avatar
Thomas Witkowski committed
489
	if (rankIt->first == globalMpiRank)
490
491
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
492
	if (rankIt->first < globalMpiRank)
493
494
495
496
497
498
499
500
501
502
503
504
	  sendBounds[rankIt->first] = rankIt->second;
	else
	  recvBounds[rankIt->first] = rankIt->second;	
      }

      stdMpi.send(sendBounds);
      stdMpi.recv(recvBounds);
      stdMpi.startCommunication();

      for (RankToBoundMap::iterator rankIt = periodic.begin();
	   rankIt != periodic.end(); ++rankIt) {

Thomas Witkowski's avatar
Thomas Witkowski committed
505
 	if (rankIt->first <= globalMpiRank)
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
 	  continue;
  
	for (unsigned int j = 0; j < rankIt->second.size(); j++) {
	  BoundaryObject &recvRankObj = 
	    stdMpi.getRecvData()[rankIt->first][j].rankObj;
	  BoundaryObject &recvNeighObj = 
	    stdMpi.getRecvData()[rankIt->first][j].neighObj;

	  if (periodic[rankIt->first][j].neighObj != recvRankObj ||
	      periodic[rankIt->first][j].rankObj != recvNeighObj) {
	    unsigned int k = j + 1;	    
	    for (; k < rankIt->second.size(); k++)
	      if (periodic[rankIt->first][k].neighObj == recvRankObj &&
		  periodic[rankIt->first][k].rankObj == recvNeighObj)
		break;
	    
	    // The element must always be found, because the list is just in 
	    // another order.
	    TEST_EXIT_DBG(k < rankIt->second.size())("Should never happen!\n");
525

526
527
528
529
530
531
532
533
	    // Swap the current with the found element.
	    AtomicBoundary tmpBound = (rankIt->second)[k];
	    (rankIt->second)[k] = (rankIt->second)[j];
	    (rankIt->second)[j] = tmpBound;	
	  } 
	}
      }     
    } // periodicBoundary.boundary.size() > 0
534
535
536
  }


537
  void InteriorBoundary::serialize(ostream &out)
538
  {
539
540
541
542
    serialize(out, own);
    serialize(out, other);
    serialize(out, periodic);
  }
543

544

545
546
547
548
549
  void InteriorBoundary::serialize(ostream &out,
				   RankToBoundMap& boundary)
  {
    FUNCNAME("InteriorBoundary::serialize()");

550
    int mSize = boundary.size();
551
    SerUtil::serialize(out, mSize);
552
553
    for (RankToBoundMap::iterator it = boundary.begin(); 
	 it != boundary.end(); ++it) {
554
555
      int rank = it->first;
      int boundSize = it->second.size();
556
557
      SerUtil::serialize(out, rank);
      SerUtil::serialize(out, boundSize);
558
559
560
      for (int i = 0; i < boundSize; i++) {
	AtomicBoundary &bound = (it->second)[i];

561
	SerUtil::serialize(out, bound.rankObj.elIndex);
562
	SerUtil::serialize(out, bound.rankObj.elType);
563
564
	SerUtil::serialize(out, bound.rankObj.subObj);
	SerUtil::serialize(out, bound.rankObj.ithObj);
565
	SerUtil::serialize(out, bound.rankObj.reverseMode);
566
	serializeExcludeList(out, bound.rankObj.excludedSubstructures);
567

568
	SerUtil::serialize(out, bound.neighObj.elIndex);
569
	SerUtil::serialize(out, bound.neighObj.elType);
570
571
	SerUtil::serialize(out, bound.neighObj.subObj);
	SerUtil::serialize(out, bound.neighObj.ithObj);
572
	SerUtil::serialize(out, bound.neighObj.reverseMode);
573
	serializeExcludeList(out, bound.neighObj.excludedSubstructures);
574
575

	SerUtil::serialize(out, bound.type);
576
577
      }
    }
578
579
  }

580

581
  void InteriorBoundary::deserialize(istream &in, Mesh *mesh)				     
582
  {
583
584
585
586
587
588
589
    map<int, Element*> elIndexMap;
    mesh->getElementIndexMap(elIndexMap);

    deserialize(in, own, elIndexMap);
    deserialize(in, other, elIndexMap);
    deserialize(in, periodic, elIndexMap);
  }
590

591

592
593
594
595
596
597
  void InteriorBoundary::deserialize(istream &in, 
				     RankToBoundMap& boundary,
				     map<int, Element*> &elIndexMap)
  {
    FUNCNAME("InteriorBoundary::deserialize()");

598
    int mSize = 0;
599
    SerUtil::deserialize(in, mSize);
600
601
602
    for (int i = 0; i < mSize; i++) {
      int rank = 0;
      int boundSize = 0;
603
604
      SerUtil::deserialize(in, rank);
      SerUtil::deserialize(in, boundSize);
605
606
607
608
609

      boundary[rank].resize(boundSize);
      for (int i = 0; i < boundSize; i++) {
	AtomicBoundary &bound = boundary[rank][i];

610
	SerUtil::deserialize(in, bound.rankObj.elIndex);
611
	SerUtil::deserialize(in, bound.rankObj.elType);
612
613
	SerUtil::deserialize(in, bound.rankObj.subObj);
	SerUtil::deserialize(in, bound.rankObj.ithObj);
614
	SerUtil::deserialize(in, bound.rankObj.reverseMode);
615
	deserializeExcludeList(in, bound.rankObj.excludedSubstructures);
616

617
	SerUtil::deserialize(in, bound.neighObj.elIndex);
618
	SerUtil::deserialize(in, bound.neighObj.elType);
619
620
	SerUtil::deserialize(in, bound.neighObj.subObj);
	SerUtil::deserialize(in, bound.neighObj.ithObj);
621
	SerUtil::deserialize(in, bound.neighObj.reverseMode);
622
	deserializeExcludeList(in, bound.neighObj.excludedSubstructures);
623

624
625
	SerUtil::deserialize(in, bound.type);

626
627
628
629
	TEST_EXIT_DBG(elIndexMap.count(bound.rankObj.elIndex) == 1)
	  ("Cannot find element with index %d for deserialization!\n", 
	   bound.rankObj.elIndex);

630
631
632
	TEST_EXIT_DBG(elIndexMap[bound.rankObj.elIndex]->getIndex() == 
		      bound.rankObj.elIndex)("Should not happen!\n");

633
	bound.rankObj.el = elIndexMap[bound.rankObj.elIndex];
634

635
636
	// For the case of periodic interior boundaries, a rank may have an
	// boundary with itself. In this case, also the pointer to the neighbour
637
	// object must be set correctly.
638
639
640
641
	if (elIndexMap.count(bound.neighObj.elIndex))
	  bound.neighObj.el = elIndexMap[bound.neighObj.elIndex];
	else
	  bound.neighObj.el = NULL;
642
643
      }
    }
644
  }
645
646


647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
  AtomicBoundary& InteriorBoundary::getNewOwn(int rank)
  {
    int size = own[rank].size();
    own[rank].resize(size + 1);
    return own[rank][size];
  }


  AtomicBoundary& InteriorBoundary::getNewOther(int rank)
  {
    int size = other[rank].size();
    other[rank].resize(size + 1);
    return other[rank][size];
  }


  AtomicBoundary& InteriorBoundary::getNewPeriodic(int rank)
  {
    int size = periodic[rank].size();
    periodic[rank].resize(size + 1);
    return periodic[rank][size];
  }


671
  void InteriorBoundary::serializeExcludeList(ostream &out, 
672
					      ExcludeList &list)
673
674
675
676
677
678
679
680
681
682
  {
    int size = list.size();
    SerUtil::serialize(out, size);
    for (int i = 0; i < size; i++) {
      SerUtil::serialize(out, list[i].first);
      SerUtil::serialize(out, list[i].second);
    }
  }


683
  void InteriorBoundary::deserializeExcludeList(istream &in, 
684
						ExcludeList &list)
685
686
687
688
689
690
691
692
693
694
695
696
  {
    int size = 0;
    SerUtil::deserialize(in, size);
    list.resize(0);
    list.reserve(size);

    for (int i = 0; i < size; i++) {
      GeoIndex a;
      int b;

      SerUtil::deserialize(in, a);
      SerUtil::deserialize(in, b);
697
      list.push_back(make_pair(a, b));
698
699
700
    }
  }

701
}