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
}