InteriorBoundary.cc 25.9 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");
38
39
    TEST_EXIT_DBG(level < levelData.getNumberOfLevels())
      ("Should not happen!\n");
40

Thomas Witkowski's avatar
Thomas Witkowski committed
41
    MPI::Intracomm mpiComm = levelData.getMpiComm(level);
42
43
44
45

    if (mpiComm == MPI::COMM_SELF)
      return;

Thomas Witkowski's avatar
Thomas Witkowski committed
46
47
    int levelMpiRank = mpiComm.Get_rank();
    int globalMpiRank = MPI::COMM_WORLD.Get_rank();
48
    std::set<int> levelRanks = levelData.getLevelRanks(level);
49

50
51
52
53
    // === Create interior boundary data structure. ===
    
    for (int geoPos = 0; geoPos < mesh->getDim(); geoPos++) {
      GeoIndex geoIndex = INDEX_OF_DIM(geoPos, mesh->getDim());
54

55
      elObjDb.resetIterator();
56
      while (elObjDb.iterate(geoIndex)) {
57
	flat_map<int, ElementObjectData>& objData = elObjDb.getIterateData();
58
59

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

63
64
65
	// 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;
66
	if (levelData.getNumberOfLevels() == 1) {
67
68
	  boundaryWithinMpiGroup = true;
	} else {
69
70
71
72
	  TEST_EXIT_DBG(level + 1 < levelData.getNumberOfLevels())
	    ("Level = %d    but number of levels = %d\n",
	     level, levelData.getNumberOfLevels());

73
	  for (flat_map<int, ElementObjectData>::iterator it = objData.begin();
74
	       it != objData.end(); ++it) {
75
	    if (!levelData.rankInLevel(it->first, level + 1)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
76
	      boundaryWithinMpiGroup = true;
77
78
79
80
	      break;
	    }
	  }
	}
81
	
82
83
84
85
	if (!boundaryWithinMpiGroup)
	  continue;

	int owner = elObjDb.getIterateOwner(level);
Thomas Witkowski's avatar
Thomas Witkowski committed
86
	ElementObjectData& rankBoundEl = objData[globalMpiRank];
87
88
89
90
91
92
93
94
	
	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
95
96


97
98
99
100
101
102
103
104
	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));
	  }
	}
105
	
Thomas Witkowski's avatar
Thomas Witkowski committed
106

Thomas Witkowski's avatar
Thomas Witkowski committed
107
	if (owner == globalMpiRank) {
108
	  for (flat_map<int, ElementObjectData>::iterator it2 = objData.begin();
109
	       it2 != objData.end(); ++it2) {
Thomas Witkowski's avatar
Thomas Witkowski committed
110
	    if (it2->first == globalMpiRank)
111
	      continue;
112

113
	    if (!levelData.rankInLevel(it2->first, level))
114
115
	      continue;

116
117
118
119
120
121
122
	    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;
123

124
125
	    int rankInLevel = levelData.mapRank(it2->first, 0, level);
	    AtomicBoundary& b = getNewOwn(rankInLevel);
126
127
128
129
130
131
132
133
134
135
136
137
	    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");
138
	  
139
	  ElementObjectData& ownerBoundEl = objData[owner];
140
  
141
142
143
144
145
146
147
	  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;
148
149
150

	  int rankInLevel = levelData.mapRank(owner, 0, level);	  
	  AtomicBoundary& b = getNewOther(rankInLevel);
151
152
153
154
155
156
157
158
	  b = bound;	    
	  if (geoIndex == EDGE)
	    b.rankObj.reverseMode =
	      elObjDb.getEdgeReverseMode(rankBoundEl, ownerBoundEl);
	  if (geoIndex == FACE)
	    b.rankObj.reverseMode = 
	      elObjDb.getFaceReverseMode(rankBoundEl, ownerBoundEl);
	}
159
      }
160
161
162
    }


163
    // === Create periodic boundary data structure. ===
164

165
166
167
168
    int removePeriodicBoundary = 0;
    Parameters::get("parallel->remove periodic boundary", removePeriodicBoundary);


169
170
    for (PerBoundMap<DegreeOfFreedom>::iterator it = elObjDb.getPeriodicVertices().begin();
	 it != elObjDb.getPeriodicVertices().end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
171
      if (elObjDb.isInRank(it->first.first, globalMpiRank) == false)
172
173
174
	continue;

      ElementObjectData& perDofEl0 = 
Thomas Witkowski's avatar
Thomas Witkowski committed
175
	elObjDb.getElementsInRank(it->first.first)[globalMpiRank];
176

177
      for (flat_map<int, ElementObjectData>::iterator elIt = elObjDb.getElementsInRank(it->first.second).begin();
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
	   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;

196

197
198
199
	if (removePeriodicBoundary) {
	  bound.type = INTERIOR;
	  
200
	  flat_map<int, ElementObjectData> objData = 
201
202
203
204
205
206
	    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));
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237

	  // Consider the following configuration of four elements partitioned
	  // to four ranks:
	  //
	  //  |--------|--------|(*)
	  //  |       /|       /|
	  //  |  3   / |  0   / |
	  //  |     /  |     /  |
	  //  |    /   |    /   |
	  //  |   /    |   /    |
	  //  |  /  2  |  /  1  |
	  //  | /      | /      |
	  //  |/-------|/-------|
	  //
	  // We are in now interested in vertex (*). For the first, rank 1 owns
	  // the interior boundary with rank 0 on this vertex. When we remove 
	  // periodic boundaries, which are applied on the left and right outer
	  // boundary, rank 3 becomes the owner of this vertex. Rank 0 and rank 1
	  // will have an interior boundary with rank 3 on this vertex, but rank
	  // 0 and rank 1 have no more interior boundaries on vertex (*). Thus we
	  // have to search and remove for this scenario.
	  
	  if (owner != globalMpiRank) {
	    removeOwn(bound.rankObj);
	    removeOther(bound.rankObj);
	  }


	  // And than we can add the boundary to either own or other part of the
	  // interior database.

238
	  if (owner == globalMpiRank) {
239
240
	    int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
	    AtomicBoundary& b = getNewOwn(rankInLevel);	  
241
	    b = bound;	  
242
	  } else {
243
	    ElementObjectData& ownerBoundEl = objData[owner];
244
245
246
247
248
	    bound.neighObj.el = elObjDb.getElementPtr(ownerBoundEl.elIndex);
	    bound.neighObj.elIndex = ownerBoundEl.elIndex;
	    bound.neighObj.elType = -1;
	    bound.neighObj.ithObj = ownerBoundEl.ithObject;

249
250
	    int rankInLevel = levelData.mapRank(owner, 0, level);
	    AtomicBoundary& b = getNewOther(rankInLevel);
251
	    b = bound;
252
	  }
253

254
	} else {
255
256
	  bound.type = it->second;
	  
257
258
	  int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
	  AtomicBoundary& b = getNewPeriodic(rankInLevel);
259
	  b = bound;	  
260
	}
261
262
263
264
265
266
      }
    }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
270
271
      ElementObjectData& perEdgeEl0 = 
	elObjDb.getElementsInRank(it->first.first)[globalMpiRank];
272

273
      for (flat_map<int, ElementObjectData>::iterator elIt = elObjDb.getElementsInRank(it->first.second).begin();
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
 	   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;
291

292
293
294
	if (removePeriodicBoundary) {
	  bound.type = INTERIOR;

295
	  flat_map<int, ElementObjectData> objData = 
296
297
298
299
300
301
302
303
	    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];

304
305
306
307
308
309
	  // See comments in the same part of code for VERTEX
	  if (owner != globalMpiRank) {
	    removeOwn(bound.rankObj);
	    removeOther(bound.rankObj);
	  }

310
	  if (owner == globalMpiRank) {
311
312
	    int rankInLevel = levelData.mapRank(owner, 0, level);
	    AtomicBoundary& b = getNewOwn(rankInLevel); 
313
	    b = bound;	  
314
	  } else {
315
316
	    int rankInLevel = levelData.mapRank(owner, 0, level);
	    AtomicBoundary& b = getNewOther(rankInLevel);
317
318
319
320
321
322
323
324
325
326
327
328
	    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;
329
330
331

	  int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
	  AtomicBoundary& b = getNewPeriodic(rankInLevel);
332
333
334
335
336
337
338
339
340
	  b = bound;
	  
	  if (globalMpiRank > otherElementRank)
	    b.neighObj.reverseMode = 
	      elObjDb.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
	  else
	    b.rankObj.reverseMode = 
	      elObjDb.getEdgeReverseMode(perEdgeEl0, perEdgeEl1);
	}
341
342
343
344
345
346
      }
    }


    for (PerBoundMap<DofFace>::iterator it = elObjDb.getPeriodicFaces().begin();
	 it != elObjDb.getPeriodicFaces().end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
347
      if (elObjDb.isInRank(it->first.first, globalMpiRank) == false)
348
349
350
351
352
353
354
	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
355
356
      ElementObjectData& perFaceEl0 = 
	elObjDb.getElementsInRank(it->first.first)[globalMpiRank];
357

358
      for (flat_map<int, ElementObjectData>::iterator elIt = elObjDb.getElementsInRank(it->first.second).begin();
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
 	   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;
	
377
378
379
	if (removePeriodicBoundary) {
	  bound.type = INTERIOR;

380
	  flat_map<int, ElementObjectData> objData = 
381
382
383
384
385
386
387
388
389
	    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) {
390
391
	    int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
	    AtomicBoundary& b = getNewOwn(rankInLevel);
392
	    b = bound;	  
393
	  } else {
394
395
	    int rankInLevel = levelData.mapRank(owner, 0, level);
	    AtomicBoundary& b = getNewOther(rankInLevel);
396
397
398
399
400
401
402
403
404
405
406
407
408
	    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;
	  
409
410
	  int rankInLevel = levelData.mapRank(otherElementRank, 0, level);
	  AtomicBoundary& b = getNewPeriodic(rankInLevel);
411
412
413
414
415
416
417
418
419
	  b = bound;
	  
	  if (globalMpiRank > otherElementRank)
	    b.neighObj.reverseMode = 
	      elObjDb.getFaceReverseMode(perFaceEl0, perFaceEl1);
	  else
	    b.rankObj.reverseMode = 
	      elObjDb.getFaceReverseMode(perFaceEl0, perFaceEl1);
	}
420
421
422
      }
    }

423

424
425
426
427
428
429
    // === 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);
430
431
    stdMpi.send(own);
    stdMpi.recv(other);
432
433
434
435
436
437
438
439
440
441
442
    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) {

443
444
445
446
447
      int rank = rankIt->first;

      // === 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.                                    ===
448
449
450
451
452
     
      for (unsigned int j = 0; j < rankIt->second.size(); j++) {

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

453
	BoundaryObject &receivedBound = 
454
	  stdMpi.getRecvData()[rank][j].rankObj;
455

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

459
460
	  for (; k < rankIt->second.size(); k++) {
 	    if ((rankIt->second)[k].neighObj == receivedBound)
461
	      break;
462
	  }
463
464
465

	  // The element must always be found, because the list is just in
	  // another order.
466
467
468
469
470
471
	  TEST_EXIT_DBG(k < rankIt->second.size())
	    ("Cannot find element %d/%d/%d received from rank %d!\n",
	     receivedBound.elIndex, 
	     receivedBound.subObj, 
	     receivedBound.ithObj,
	     rank);
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490

	  // 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
491
	if (rankIt->first == globalMpiRank)
492
493
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
494
	if (rankIt->first < globalMpiRank)
495
496
497
498
499
500
501
502
503
504
505
506
	  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
507
 	if (rankIt->first <= globalMpiRank)
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
 	  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");
527

528
529
530
531
532
533
534
535
	    // 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
536
537
538
539
540


    // === Run verification procedure. ===
    
    verifyBoundary();
541
542
543
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
  int InteriorBoundary::getDegreeOwn(BoundaryObject &bObj)
  {
    int counter = 0;

    for (RankToBoundMap::iterator it = own.begin(); it != own.end(); ++it) {
      for (vector<AtomicBoundary>::iterator bIt = it->second.begin(); 
	   bIt != it->second.end(); ++bIt) {
	if (bIt->rankObj == bObj) {
	  counter++;
	  break;
	}
      }
    }

    return counter;
  }


562
  void InteriorBoundary::serialize(ostream &out)
563
  {
564
565
566
567
    serialize(out, own);
    serialize(out, other);
    serialize(out, periodic);
  }
568

569

570
571
572
573
574
  void InteriorBoundary::serialize(ostream &out,
				   RankToBoundMap& boundary)
  {
    FUNCNAME("InteriorBoundary::serialize()");

575
    int mSize = boundary.size();
576
    SerUtil::serialize(out, mSize);
577
578
    for (RankToBoundMap::iterator it = boundary.begin(); 
	 it != boundary.end(); ++it) {
579
580
      int rank = it->first;
      int boundSize = it->second.size();
581
582
      SerUtil::serialize(out, rank);
      SerUtil::serialize(out, boundSize);
583
584
585
      for (int i = 0; i < boundSize; i++) {
	AtomicBoundary &bound = (it->second)[i];

586
	SerUtil::serialize(out, bound.rankObj.elIndex);
587
	SerUtil::serialize(out, bound.rankObj.elType);
588
589
	SerUtil::serialize(out, bound.rankObj.subObj);
	SerUtil::serialize(out, bound.rankObj.ithObj);
590
	SerUtil::serialize(out, bound.rankObj.reverseMode);
591
	serializeExcludeList(out, bound.rankObj.excludedSubstructures);
592

593
	SerUtil::serialize(out, bound.neighObj.elIndex);
594
	SerUtil::serialize(out, bound.neighObj.elType);
595
596
	SerUtil::serialize(out, bound.neighObj.subObj);
	SerUtil::serialize(out, bound.neighObj.ithObj);
597
	SerUtil::serialize(out, bound.neighObj.reverseMode);
598
	serializeExcludeList(out, bound.neighObj.excludedSubstructures);
599
600

	SerUtil::serialize(out, bound.type);
601
602
      }
    }
603
604
  }

605

606
  void InteriorBoundary::deserialize(istream &in, Mesh *mesh)				     
607
  {
608
609
610
611
612
613
614
    map<int, Element*> elIndexMap;
    mesh->getElementIndexMap(elIndexMap);

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

616

617
618
619
620
621
622
  void InteriorBoundary::deserialize(istream &in, 
				     RankToBoundMap& boundary,
				     map<int, Element*> &elIndexMap)
  {
    FUNCNAME("InteriorBoundary::deserialize()");

623
    int mSize = 0;
624
    SerUtil::deserialize(in, mSize);
625
626
627
    for (int i = 0; i < mSize; i++) {
      int rank = 0;
      int boundSize = 0;
628
629
      SerUtil::deserialize(in, rank);
      SerUtil::deserialize(in, boundSize);
630
631
632
633
634

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

635
	SerUtil::deserialize(in, bound.rankObj.elIndex);
636
	SerUtil::deserialize(in, bound.rankObj.elType);
637
638
	SerUtil::deserialize(in, bound.rankObj.subObj);
	SerUtil::deserialize(in, bound.rankObj.ithObj);
639
	SerUtil::deserialize(in, bound.rankObj.reverseMode);
640
	deserializeExcludeList(in, bound.rankObj.excludedSubstructures);
641

642
	SerUtil::deserialize(in, bound.neighObj.elIndex);
643
	SerUtil::deserialize(in, bound.neighObj.elType);
644
645
	SerUtil::deserialize(in, bound.neighObj.subObj);
	SerUtil::deserialize(in, bound.neighObj.ithObj);
646
	SerUtil::deserialize(in, bound.neighObj.reverseMode);
647
	deserializeExcludeList(in, bound.neighObj.excludedSubstructures);
648

649
650
	SerUtil::deserialize(in, bound.type);

651
652
653
654
	TEST_EXIT_DBG(elIndexMap.count(bound.rankObj.elIndex) == 1)
	  ("Cannot find element with index %d for deserialization!\n", 
	   bound.rankObj.elIndex);

655
656
657
	TEST_EXIT_DBG(elIndexMap[bound.rankObj.elIndex]->getIndex() == 
		      bound.rankObj.elIndex)("Should not happen!\n");

658
	bound.rankObj.el = elIndexMap[bound.rankObj.elIndex];
659

660
661
	// For the case of periodic interior boundaries, a rank may have an
	// boundary with itself. In this case, also the pointer to the neighbour
662
	// object must be set correctly.
663
664
665
666
	if (elIndexMap.count(bound.neighObj.elIndex))
	  bound.neighObj.el = elIndexMap[bound.neighObj.elIndex];
	else
	  bound.neighObj.el = NULL;
667
668
      }
    }
669
  }
670
671


672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
  void InteriorBoundary::verifyBoundary()
  {
    FUNCNAME("InteriorBoundary::verifyBoundary()");

    // === Check if no other boundery rank object is also included in own ===
    // === boundary part, which would make no sence.                      ===

    std::set<BoundaryObject> allOwnBounds;
    for (map<int, vector<AtomicBoundary> >::iterator it = own.begin(); 
	 it != own.end(); ++it)
      for (vector<AtomicBoundary>::iterator bIt = it->second.begin();
	   bIt != it->second.end(); ++bIt)
	allOwnBounds.insert(bIt->rankObj);

    for (map<int, vector<AtomicBoundary> >::iterator it = other.begin(); 
	 it != other.end(); ++it)
      for (vector<AtomicBoundary>::iterator bIt = it->second.begin();
	   bIt != it->second.end(); ++bIt)
	if (allOwnBounds.count(bIt->rankObj)) {
	  ERROR_EXIT("Boundary %d/%d/%d is included in both, own and other interior boundary part!\n",
		     bIt->rankObj.elIndex,
		     bIt->rankObj.subObj,
		     bIt->rankObj.ithObj);
	}
  }


699
700
  AtomicBoundary& InteriorBoundary::getNewOwn(int rank)
  {
701
702
    FUNCNAME("InteriorBoundary::getNewOwn()");

703
704
705
706
707
708
709
710
    int size = own[rank].size();
    own[rank].resize(size + 1);
    return own[rank][size];
  }


  AtomicBoundary& InteriorBoundary::getNewOther(int rank)
  {
711
712
    FUNCNAME("InteriorBoundary::getNewOther()");

713
714
715
716
717
718
719
720
    int size = other[rank].size();
    other[rank].resize(size + 1);
    return other[rank][size];
  }


  AtomicBoundary& InteriorBoundary::getNewPeriodic(int rank)
  {
721
722
    FUNCNAME("InteriorBoundary::getNewPeriodic()");

723
724
725
726
727
728
    int size = periodic[rank].size();
    periodic[rank].resize(size + 1);
    return periodic[rank][size];
  }


729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
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
787
788
789
790
791
792
793
794
795
796
797
798
799
800
  bool InteriorBoundary::checkOther(AtomicBoundary& bound, int rank)
  {
    FUNCNAME("InteriorBoundary::checkOther()");

    int globalMpiRank = MPI::COMM_WORLD.Get_rank();
    TEST_EXIT(rank > globalMpiRank)
      ("Wrong MPI ranks: %d %d\n", rank, globalMpiRank);

    bool found = false;

    if (other.count(rank)) {
      vector<AtomicBoundary>& otherBounds = other[rank];
      for (int i = 0; i < static_cast<int>(otherBounds.size()); i++) {
	if (otherBounds[i].rankObj == bound.rankObj) {
	  if (otherBounds[i].neighObj != bound.neighObj) {
	    ERROR_EXIT("If this really can happen, I have to thing about this!\n");
	  } else {
	    found = true;
	    break;
	  }
	}
      }
    }

    return found;
  }


  bool InteriorBoundary::removeOwn(BoundaryObject& bound)
  {
    FUNCNAME("InteriorBoundary::removeOwn()");

    bool removed = false;

    for (map<int, vector<AtomicBoundary> >::iterator it = own.begin();
	 it != own.end(); ++it) {
      for (vector<AtomicBoundary>::iterator bIt = it->second.begin();
	   bIt != it->second.end(); ++bIt) {
	if (bIt->rankObj == bound) {
	  it->second.erase(bIt);
	  removed = true;
	  break;
	}
      }
    }

    return removed;
  }


  bool InteriorBoundary::removeOther(BoundaryObject& bound)
  {
    FUNCNAME("InteriorBoundary::removeOther()");

    bool removed = false;

    for (map<int, vector<AtomicBoundary> >::iterator it = other.begin();
	 it != other.end(); ++it) {
      for (vector<AtomicBoundary>::iterator bIt = it->second.begin();
	   bIt != it->second.end(); ++bIt) {
	if (bIt->rankObj == bound) {
	  it->second.erase(bIt);
	  removed = true;
	  break;
	}
      }
    }

    return removed;    
  }


801
  void InteriorBoundary::serializeExcludeList(ostream &out, 
802
					      ExcludeList &list)
803
804
805
806
807
808
809
810
811
812
  {
    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);
    }
  }


813
  void InteriorBoundary::deserializeExcludeList(istream &in, 
814
						ExcludeList &list)
815
816
817
818
819
820
821
822
823
824
825
826
  {
    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);
827
      list.push_back(make_pair(a, b));
828
829
830
    }
  }

831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857

  void MultiLevelInteriorBoundary::create(MeshLevelData &levelData,
					  ElementObjectDatabase &elObjDb)
  {
    FUNCNAME(" MultLevelInteriorBoundary::create()");

    levelIntBound.clear();

    int nLevel = levelData.getNumberOfLevels();
    for (int level = 0; level < nLevel; level++)
      levelIntBound[level].create(levelData, level, elObjDb);
  }


  void MultiLevelInteriorBoundary::serialize(ostream &out)
  {
    FUNCNAME("MultiLevelInteriorBoundary::serialize()");
    ERROR_EXIT("Not yet implemented!\n");
  }


  void MultiLevelInteriorBoundary::deserialize(istream &in, Mesh *mesh)
  {
    FUNCNAME("MultiLevelInteriorBoundary::deserialize()");
    ERROR_EXIT("Not yet implemented!\n");
  }

858
}