ParallelDomainBase.cc 75 KB
Newer Older
1001
1002
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>(element->getElementData(PARTITION_ED));   
1003

1004
      // Check, if the element is within rank's partition.
1005
      if (partitionData->getPartitionStatus() == IN) {
1006
	for (int i = 0; i < nNeighbours; i++) {
1007
1008
1009
1010
	  if (!elInfo->getNeighbour(i))
	    continue;

	  PartitionElementData *neighbourPartitionData =
1011
1012
	    dynamic_cast<PartitionElementData*>(elInfo->getNeighbour(i)->
						getElementData(PARTITION_ED));
1013

1014
 	  if (neighbourPartitionData->getPartitionStatus() == OUT) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1015

1016
1017
1018
1019
1020
	    // We have found an element that is in rank's partition, but has a 
	    // neighbour outside of the rank's partition.

	    // === Create information about the boundary between the two elements. ===

1021
	    AtomicBoundary bound;	    	    
1022
1023
	    bound.rankObj.el = element;
	    bound.rankObj.elIndex = element->getIndex();
1024
	    bound.rankObj.subObj = subObj;
1025
	    bound.rankObj.ithObj = i;
1026
1027
	    // Do not set a pointer to the element, because if will be deleted from
	    // mesh after partitioning and the pointer would become invalid.
1028
1029
	    bound.neighObj.el = NULL;
	    bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex();
1030
	    bound.neighObj.subObj = subObj;
1031
	    bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i);
1032
	    
1033
1034
1035
1036
1037
	    if (dim == 2) {
	      // i == 2  =>  getSideOfNeighbour(i) == 2
	      TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2)
		("Should not happen!\n");
	    }
1038

1039
	    // Get rank number of the neighbouring element.
1040
1041
	    int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];

1042
1043
1044
1045
1046
1047
1048
1049
1050

	    // === Add the boundary information object to the corresponding overall ===
	    // === boundary. There are three possibilities: if the boundary is a    ===
	    // === periodic boundary, just add it to \ref periodicBounadry. Here it ===
	    // === does not matter which rank is responsible for this boundary.     ===
	    // === Otherwise, the boundary is added either to \ref myIntBoundary or ===
	    // === to \ref otherIntBoundary. It dependes on which rank is respon-   ===
	    // === sible for this boundary.                                         ===

1051
	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(subObj, i))) {	      
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
	      // We have found an element that is at an interior, periodic boundary.
	      AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank);
	      b = bound;
	    } else {
	      // We have found an element that is at an interior, non-periodic boundary.
	      AtomicBoundary& b = ((mpiRank > otherElementRank) ?
				   myIntBoundary.getNewAtomic(otherElementRank) :
				   otherIntBoundary.getNewAtomic(otherElementRank));
	      b = bound;	      
	    }
1062
1063
1064
1065
1066
1067
 	  }
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1068

1069
1070
1071
1072

    // === Once we have this information, we must care about the order of the atomic ===
    // === bounds in the three boundary handling object. Eventually all the bound-   ===
    // === aries have to be in the same order on both ranks that share the bounday.  ===
Thomas Witkowski's avatar
Thomas Witkowski committed
1073

Thomas Witkowski's avatar
Thomas Witkowski committed
1074
    std::vector<int*> sendBuffers, recvBuffers;
Thomas Witkowski's avatar
Thomas Witkowski committed
1075

1076
    // First we communicate myInt/otherIntBoundary, and than the periodic boundaries.
1077
1078
1079
    MPI::Request request[max(myIntBoundary.boundary.size() + 
			     otherIntBoundary.boundary.size(),
			     periodicBoundary.boundary.size())];
Thomas Witkowski's avatar
Thomas Witkowski committed
1080
1081
    int requestCounter = 0;

1082

1083
1084
1085
    // === The owner of the boundaries send their atomic boundary order to the ===
    // === neighbouring ranks.                                                 ===

Thomas Witkowski's avatar
Thomas Witkowski committed
1086
1087
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
1088
      int nSendInt = rankIt->second.size();
1089
1090
      int* buffer = new int[nSendInt * 2];
      for (int i = 0; i < nSendInt; i++) {
1091
1092
	buffer[i * 2] = (rankIt->second)[i].rankObj.elIndex;
	buffer[i * 2 + 1] = (rankIt->second)[i].rankObj.ithObj;
1093
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
1094
1095
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
1096
      request[requestCounter++] =
1097
	mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1098
1099
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1100
1101
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1102
      int nRecvInt = rankIt->second.size();
1103
      int *buffer = new int[nRecvInt * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
1104
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
1105

Thomas Witkowski's avatar
Thomas Witkowski committed
1106
      request[requestCounter++] = 
1107
	mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1108
1109
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1110
1111
1112

    MPI::Request::Waitall(requestCounter, request);

1113
1114
1115
1116
    
    // === 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.             ===
Thomas Witkowski's avatar
Thomas Witkowski committed
1117
1118

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1119
1120
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1121
1122
1123
1124
1125
1126
1127

      // === We have received from rank "rankIt->first" the ordered list of element ===
      // === indices. We now have to sort the corresponding list in this rank to    ===
      // === get the same order.                                                    ===
      
      for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
	// If the expected object is not at place, search for it.
1128
1129
	if ((rankIt->second)[j].neighObj.elIndex != recvBuffers[i][j * 2] ||
	    (rankIt->second)[j].neighObj.ithObj != recvBuffers[i][j * 2 + 1]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1130
1131
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
1132
1133
	    if ((rankIt->second)[k].neighObj.elIndex == recvBuffers[i][j * 2] && 
		(rankIt->second)[k].neighObj.ithObj == recvBuffers[i][j * 2 + 1])
Thomas Witkowski's avatar
Thomas Witkowski committed
1134
1135
1136
	      break;

	  // The element must always be found, because the list is just in another order.
Thomas Witkowski's avatar
Thomas Witkowski committed
1137
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
1138
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
1139
1140
1141
1142
1143
1144
1145
1146

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

1147
1148
      delete [] recvBuffers[i];
      i++;
Thomas Witkowski's avatar
Thomas Witkowski committed
1149
1150
1151
1152
    }

    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];
1153
1154
1155
1156
1157

    recvBuffers.clear();
    sendBuffers.clear();
 

1158
    // === Do the same for the periodic boundaries. ===
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171

    requestCounter = 0;

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

      TEST_EXIT_DBG(rankIt->first != mpiRank)
	("It is no possible to have an interior boundary within a rank partition!\n");

      if (rankIt->first < mpiRank) {
	int nSendInt = rankIt->second.size();
	int* buffer = new int[nSendInt * 2];
	for (int i = 0; i < nSendInt; i++) {
1172
1173
	  buffer[i * 2] = (rankIt->second)[i].rankObj.elIndex;
	  buffer[i * 2 + 1] = (rankIt->second)[i].rankObj.ithObj;
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
	}
	sendBuffers.push_back(buffer);

	request[requestCounter++] =
	  mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
      } else {
	int nRecvInt = rankIt->second.size();
	int *buffer = new int[nRecvInt * 2];
	recvBuffers.push_back(buffer);
	
	request[requestCounter++] = 
	  mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);	
      }
    }


    MPI::Request::Waitall(requestCounter, request);


    i = 0;
    for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin();
	 rankIt != periodicBoundary.boundary.end(); ++rankIt) {
      if (rankIt->first > mpiRank) {

	for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) 
1199
1200
	  if (periodicBoundary.boundary[rankIt->first][j].neighObj.elIndex != recvBuffers[i][j * 2] ||
	      periodicBoundary.boundary[rankIt->first][j].neighObj.ithObj != recvBuffers[i][j * 2 + 1]) {
1201
1202
1203
	    int k = j + 1;

	    for (; k < static_cast<int>(rankIt->second.size()); k++)
1204
1205
	      if (periodicBoundary.boundary[rankIt->first][k].neighObj.elIndex == recvBuffers[i][j * 2] && 
		  periodicBoundary.boundary[rankIt->first][k].neighObj.ithObj == recvBuffers[i][j * 2 + 1])
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
		break;

	    // The element must always be found, because the list is just in another order.
	    TEST_EXIT_DBG(k < static_cast<int>(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;	
	  }

	delete [] recvBuffers[i];
	i++;
      }
    }

    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i]; 
1225
1226
1227
  }


1228
  void ParallelDomainBase::removeMacroElements()
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
  {
    std::vector<MacroElement*> macrosToRemove;
    for (std::deque<MacroElement*>::iterator it = mesh->firstMacroElement();
	 it != mesh->endOfMacroElements();
	 ++it) {
      PartitionElementData *partitionData = 
	dynamic_cast<PartitionElementData*>
	((*it)->getElement()->getElementData(PARTITION_ED));
      if (partitionData->getPartitionStatus() != IN)
	macrosToRemove.push_back(*it);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1241
    mesh->removeMacroElements(macrosToRemove, feSpace);
1242
1243
1244
  }


1245
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDofs,
1246
						      int& nRankDofs, 
Thomas Witkowski's avatar
Thomas Witkowski committed
1247
						      int& nOverallDOFs)
1248
  {
1249
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
1250
1251

    // === Get rank information about DOFs. ===
1252
1253
1254

    // Stores to each DOF pointer the set of ranks the DOF is part of.
    std::map<const DegreeOfFreedom*, std::set<int> > partitionDOFs;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1255
    DofContainer rankAllDofs;
1256
    DofToRank boundaryDofs;
1257

1258
    createDofMemberInfo(partitionDOFs, rankDofs, rankAllDofs, boundaryDofs, vertexDof);
1259

1260
    nRankDofs = rankDofs.size();
1261
1262
    nOverallDOFs = partitionDOFs.size();

1263

1264
    // === Get starting position for global rank dof ordering. ====
1265

1266
    rstart = 0;
1267
1268
    mpiComm.Scan(&nRankDofs, &rstart, 1, MPI_INT, MPI_SUM);
    rstart -= nRankDofs;
1269

1270
1271
1272
1273
    
    // === Create for all dofs in rank new indices. The new index must start at ===
    // === index 0, must be continues and have the same order as the indices    ===
    // === had before.                                                          ===
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1274

Thomas Witkowski's avatar
Thomas Witkowski committed
1275

1276
    // Do not change the indices now, but create a new indexing and store it here.
1277
1278
    DofIndexMap rankDofsNewLocalIndex;
    isRankDof.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1279
1280
1281
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
1282
1283
1284
      rankDofsNewLocalIndex[*dofIt] = i;
      // First, we set all dofs in ranks partition to be owend by the rank. Later,
      // the dofs in ranks partition that are owned by other rank are set to false.
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1285
      isRankDof[i] = true;
1286
      i++;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1287
1288
    }

1289
    // === Create for all rank owned dofs a new global indexing. ===
1290

1291
1292
    // Stores for dofs in rank a new global index.
    DofIndexMap rankDofsNewGlobalIndex;
1293
1294
1295
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1296
    i = 0;
1297
1298
    for (DofContainer::iterator dofIt = rankDofs.begin();
	 dofIt != rankDofs.end(); ++dofIt) {
1299
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
1300
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1301
1302
1303
      i++;
    }

1304
1305
1306
1307
1308
    // === Create information which dof indices must be send and which must  ===
    // === be received.                                                      ===

    // Stores to each rank a map from DOF pointers (which are all owned by the rank
    // and lie on an interior boundary) to new global DOF indices.
1309
    std::map<int, DofIndexMap> sendNewDofs;
1310

1311
1312
1313
1314
    // Maps to each rank the number of new DOF indices this rank will receive from.
    // All these DOFs are at an interior boundary on this rank, but are owned by
    // another rank.
    std::map<int, int> recvNewDofs;
1315

1316
    for (DofToRank::iterator it = boundaryDofs.begin(); it != boundaryDofs.end(); ++it) {
1317
1318
1319
1320
1321
1322
1323
1324

      if (it->second == mpiRank) {
	// If the boundary dof is a rank dof, it must be send to other ranks.

	// Search for all ranks that have this dof too.
	for (std::set<int>::iterator itRanks = partitionDOFs[it->first].begin();
	     itRanks != partitionDOFs[it->first].end();
	     ++itRanks) {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1325
	  if (*itRanks != mpiRank) {
1326
	    TEST_EXIT_DBG(rankDofsNewGlobalIndex.count(it->first) == 1)
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1327
1328
	      ("DOF Key not found!\n");

1329
	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1330
	  }
1331
1332
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1333
1334
	// If the boundary dof is not a rank dof, its new dof index (and later
	// also the dof values) must be received from another rank.
1335
1336
1337
1338
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
1339
1340
1341
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1342

1343
    // === Send and receive the dof indices at boundary. ===
1344

1345
    sendDofs.clear();
1346
    for (std::map<int, DofIndexMap>::iterator sendIt = sendNewDofs.begin();
1347
	 sendIt != sendNewDofs.end(); ++sendIt)
1348
      for (DofIndexMap::iterator dofIt = sendIt->second.begin();
1349
	   dofIt != sendIt->second.end(); ++dofIt)
1350
	sendDofs[sendIt->first].push_back(dofIt->first);
1351
1352


1353
    StdMpi<DofIndexMap, std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > > stdMpi(mpiComm);
1354
1355
    stdMpi.prepareCommunication(false);
    stdMpi.send(sendNewDofs);
1356
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
1357
1358
1359
	 recvIt != recvNewDofs.end(); ++recvIt, i++)
      stdMpi.recv(recvIt->first, recvIt->second * 2); 
    stdMpi.startCommunication();
1360
    std::map<int, std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > > &dofMap = stdMpi.getRecvData();
1361

1362
    // === Change dof indices at boundary from other ranks. ===
1363

Thomas Witkowski's avatar
Thomas Witkowski committed
1364
1365
    // Within this small data structure we track which dof index was already changed.
    // This is used to avoid the following situation: Assume, there are two dof indices
1366
    // a and b in boundaryDofs. Then we have to change index a to b and b to c. When
Thomas Witkowski's avatar
Thomas Witkowski committed
1367
1368
1369
1370
1371
    // the second rule applies, we have to avoid that not the first b, resulted from
    // changing a to b, is set to c, but the second one. Therefore, after the first
    // rule was applied, the dof pointer is set to false in this data structure and 
    // is not allowed to be changed anymore.
    std::map<const DegreeOfFreedom*, bool> dofChanged;
1372
1373
    recvDofs.clear();

1374
    for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
1375
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
1376
1377
      dofChanged[dofIt->first] = false;

1378
    for (std::map<int, std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > >::iterator rankIt = dofMap.begin();
1379
1380
 	 rankIt != dofMap.end(); ++rankIt) {
      
1381
      for (std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> >::iterator recvDofIt = rankIt->second.begin();
1382
	   recvDofIt != rankIt->second.end(); ++recvDofIt) {
1383

1384
1385
1386
	DegreeOfFreedom oldDof = recvDofIt->first;
	DegreeOfFreedom newGlobalDof = recvDofIt->second;
	
1387
	bool found = false;
1388
	
1389
	// Iterate over all boundary dofs to find the dof which index we have to change.
1390
	
1391
	for (DofToRank::iterator dofIt = boundaryDofs.begin(); 
1392
1393
	     dofIt != boundaryDofs.end(); ++dofIt) {   
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
1394
1395
	  if (*(dofIt->first) == oldDof && !dofChanged[dofIt->first]) {
	    dofChanged[dofIt->first] = true;
1396
1397
	    
	    recvDofs[rankIt->first].push_back(dofIt->first);
1398
	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1399
	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
1400
	    
1401
	    found = true;
1402
1403
1404
	    break;
	  }
	}
1405
	
Thomas Witkowski's avatar
Thomas Witkowski committed
1406
	TEST_EXIT_DBG(found)("Should not happen!\n");
1407
      }      
1408
    }
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1409

1410
    // === Create now the local to global index and local to dof index mappings.  ===
1411

1412
1413
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
1414
1415
1416
  }


1417
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDofs, int& nOverallDOFs)
1418
  {
1419
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
1420

1421
    typedef std::set<const DegreeOfFreedom*> DofSet;
Thomas Witkowski's avatar
Thomas Witkowski committed
1422

1423
    // === Get all DOFs in ranks partition. ===
1424

Thomas Witkowski's avatar
Thomas Witkowski committed
1425
    ElementDofIterator elDofIt(feSpace);
1426
    DofSet rankDofSet;
1427
    
1428
1429
1430
    // The vertexDof list must be recreated from the scratch. Otherwise, it is possible
    // that it maps dofs, that were removed (this is also possible, if the mesh was
    // refined, e.g., center dofs of an element are not dofs of the children).
1431
    DofToBool oldVertexDof = vertexDof;
1432
1433
    vertexDof.clear();

1434
1435
1436
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1437
      Element *element = elInfo->getElement();     
Thomas Witkowski's avatar
Thomas Witkowski committed
1438
1439
      elDofIt.reset(element);
      do {
1440
1441
1442
1443
1444
1445
	rankDofSet.insert(elDofIt.getDofPtr());

	if (elDofIt.getCurrentPos() == 0) 
	  vertexDof[elDofIt.getDofPtr()] = true;
	else
	  vertexDof[elDofIt.getDofPtr()] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
1446
      } while(elDofIt.next());
1447
1448
1449
1450

      elInfo = stack.traverseNext(elInfo);
    }

1451
    DofContainer rankAllDofs;
1452
    for (DofSet::iterator dofIt = rankDofSet.begin(); dofIt != rankDofSet.end(); ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
1453
      rankAllDofs.push_back(*dofIt);    
1454
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
1455
    DofContainer rankDofs = rankAllDofs;
1456
1457


Thomas Witkowski's avatar
Thomas Witkowski committed
1458
    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
1459
    // === rankDofs to boundaryDOFs.                                               ===
1460

1461
1462
1463
    RankToDofContainer oldSendDofs = sendDofs;
    RankToDofContainer oldRecvDofs = recvDofs;

1464
1465
    sendDofs.clear();
    recvDofs.clear();
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489

    for (RankToDofContainer::iterator it = oldSendDofs.begin();
	 it != oldSendDofs.end(); ++it) {
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	if (oldVertexDof[*dofIt])
	  sendDofs[it->first].push_back(*dofIt);
      }
    }

    for (RankToDofContainer::iterator it = oldRecvDofs.begin();
	 it != oldRecvDofs.end(); ++it) {
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	if (oldVertexDof[*dofIt]) {
	  recvDofs[it->first].push_back(*dofIt);

	  DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), *dofIt);
	  if (eraseIt != rankDofs.end())
	    rankDofs.erase(eraseIt);    
	}
      }
    }

1490
    
1491
1492
1493
1494
    for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end(); ++it) {    
      
      DofContainer &dofsToSend = sendDofs[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
1495

Thomas Witkowski's avatar
Thomas Witkowski committed
1496
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1497
	   boundIt != it->second.end(); ++boundIt) {
1498
	DofContainer dofs;	
1499
1500
1501
	boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);
	boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);

1502
1503
1504
1505
	for (int i = 0; i < static_cast<int>(dofs.size()); i++) {
	  TEST_EXIT_DBG(find(dofsToSend.begin(), dofsToSend.end(), dofs[i]) == dofsToSend.end())
	    ("Should not happen!\n");

1506
	  dofsToSend.push_back(dofs[i]);
1507
	}
1508
      }
1509
    }
1510
   
1511

Thomas Witkowski's avatar
Thomas Witkowski committed
1512
1513
1514
    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end(); ++it) {

1515
      DofContainer &dofsToRecv = recvDofs[it->first];
1516

1517
1518
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
1519

1520
	DofContainer dofs;
1521
1522
	boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);
	boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);
1523
1524

	for (int i = static_cast<int>(dofs.size()) - 1; i >= 0; i--) {
1525
1526
1527
	  TEST_EXIT_DBG(find(dofsToRecv.begin(), dofsToRecv.end(), dofs[i]) == dofsToRecv.end())
	    ("Should not happen!\n");

1528
1529
1530
	  DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), dofs[i]);
	  if (eraseIt != rankDofs.end())
	    rankDofs.erase(eraseIt);	  
1531

1532
	  dofsToRecv.push_back(dofs[i]);
1533
	}
1534

Thomas Witkowski's avatar
Thomas Witkowski committed
1535
1536
1537
      }
    }

1538
    nRankDofs = rankDofs.size();
1539
1540
1541

    // === Get starting position for global rank dof ordering. ====

1542
1543
    mpiComm.Scan(&nRankDofs, &rstart, 1, MPI_INT, MPI_SUM);
    rstart -= nRankDofs;
1544
1545
1546
1547


    // === Calculate number of overall DOFs of all partitions. ===

1548
    mpiComm.Allreduce(&nRankDofs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
1549
1550


1551
1552
1553
    // Do not change the indices now, but create a new indexing a store it here.
    DofIndexMap rankDofsNewLocalIndex;
    isRankDof.clear();
1554
    int i = 0;
1555
1556
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
1557

1558
1559
1560
      rankDofsNewLocalIndex[*dofIt] = i;
      // First, we set all dofs in ranks partition to be owend by the rank. Later,
      // the dofs in ranks partition that are owned by other rank are set to false.
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1561
      isRankDof[i] = true;
1562
      i++;
1563
1564
    }

1565
    // Stores for all rank owned dofs a new global index.
1566
    DofIndexMap rankDofsNewGlobalIndex;
1567
1568
1569
1570
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

    i = 0;
1571
1572
    for (DofContainer::iterator dofIt = rankDofs.begin();
	 dofIt != rankDofs.end(); ++dofIt) {
1573
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
1574
1575
1576
1577
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
      i++;
    }

1578
1579
    // === Send new DOF indices. ===

1580
1581
    std::vector<int*> sendBuffers(sendDofs.size());
    std::vector<int*> recvBuffers(recvDofs.size());
1582

1583
    MPI::Request request[sendDofs.size() + recvDofs.size()];
Thomas Witkowski's avatar
Thomas Witkowski committed
1584
1585
    int requestCounter = 0;

1586
    i = 0;
1587
1588
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt, i++) {
1589
1590
1591
      int nSendDofs = sendIt->second.size();
      sendBuffers[i] = new int[nSendDofs];
      int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1592
      for (DofContainer::iterator dofIt = sendIt->second.begin();
1593
	   dofIt != sendIt->second.end(); ++dofIt)
1594
	sendBuffers[i][c++] = rankDofsNewGlobalIndex[*dofIt];
1595

Thomas Witkowski's avatar
Thomas Witkowski committed
1596
      request[requestCounter++] = 
1597
 	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
1598
1599
1600
    }

    i = 0;
1601
1602
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt, i++) {
1603
1604
      int nRecvDofs = recvIt->second.size();
      recvBuffers[i] = new int[nRecvDofs];
1605

Thomas Witkowski's avatar
Thomas Witkowski committed
1606
      request[requestCounter++] = 
1607
 	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
1608
1609
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1610
1611
    MPI::Request::Waitall(requestCounter, request);

1612
1613
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
1614
1615

    i = 0;
1616
1617
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt) {      
1618
      int j = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1619
      for (DofContainer::iterator dofIt = recvIt->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1620
	   dofIt != recvIt->second.end(); ++dofIt) {
1621

1622
	rankDofsNewGlobalIndex[*dofIt] = recvBuffers[i][j];     
1623
	isRankDof[rankDofsNewLocalIndex[*dofIt]] = false;
1624
1625
1626
1627
	j++;
      }

      delete [] recvBuffers[i++];
1628
1629
    }

1630
    // === Create now the local to global index and local to dof index mappings.  ===
1631

1632
1633
1634
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
  }
1635

1636
1637
1638
1639
  void ParallelDomainBase::createLocalMappings(DofIndexMap &rankDofsNewLocalIndex,
					       DofIndexMap &rankOwnedDofsNewLocalIndex,
					       DofIndexMap &rankDofsNewGlobalIndex)
  {
1640
    mapLocalGlobalDOFs.clear();
1641
    mapLocalToDofIndex.clear();
1642
1643
1644
1645

    for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin();
	 dofIt != rankDofsNewLocalIndex.end(); ++dofIt) {
      DegreeOfFreedom localDof = dofIt->second;
1646
      DegreeOfFreedom globalDof = rankDofsNewGlobalIndex[dofIt->first];
1647
1648
1649
1650
1651
1652
1653
1654

      *const_cast<DegreeOfFreedom*>(dofIt->first) = localDof;
      mapLocalGlobalDOFs[localDof] = globalDof;
    }

    for (DofIndexMap::iterator dofIt = rankOwnedDofsNewLocalIndex.begin();
	 dofIt != rankOwnedDofsNewLocalIndex.end(); ++dofIt)
      mapLocalToDofIndex[dofIt->second] = *(dofIt->first);
1655
1656
  }

1657
  
1658
  void ParallelDomainBase::createDofMemberInfo(DofToPartitions& partitionDofs,
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1659
1660
					       DofContainer& rankOwnedDofs,
					       DofContainer& rankAllDofs,
1661
1662
					       DofToRank& boundaryDofs,
					       DofToBool& vertexDof)
1663
  {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1664
1665
1666
1667
    partitionDofs.clear();
    rankOwnedDofs.clear();
    rankAllDofs.clear();
    boundaryDofs.clear();
1668
    vertexDof.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1669

1670
    // === Determine to each dof the set of partitions the dof belongs to. ===
1671

Thomas Witkowski's avatar
Thomas Witkowski committed
1672
    ElementDofIterator elDofIt(feSpace);
1673

1674
1675
1676
1677
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *element = elInfo->getElement();
1678
1679
1680
      elDofIt.reset(element);
      do {
	// Determine to each dof the partition(s) it corresponds to.
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1681
	partitionDofs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
1682
1683
1684
1685
1686

	if (elDofIt.getCurrentPos() == 0) 
	  vertexDof[elDofIt.getDofPtr()] = true;
	else
	  vertexDof[elDofIt.getDofPtr()] = false;
1687
      } while (elDofIt.next());
1688

1689
1690
1691
      elInfo = stack.traverseNext(elInfo);
    }

1692

1693
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
1694

1695
    // iterate over all DOFs
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1696
1697
    for (DofToPartitions::iterator it = partitionDofs.begin();
	 it != partitionDofs.end(); ++it) {
1698

1699
1700
      bool isInRank = false;

1701
      // iterate over all partition the current DOF is part of.
1702
      for (std::set<int>::iterator itpart1 = it->second.begin();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1703
	   itpart1 != it->second.end(); ++itpart1) {
1704

1705
	if (*itpart1 == mpiRank) {	  
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1706
1707
	  rankAllDofs.push_back(it->first);

1708
	  if (it->second.size() == 1) {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1709
	    rankOwnedDofs.push_back(it->first);
1710
1711
1712
1713
1714
1715
1716
	  } else {	    
	    // This dof is at the ranks boundary. It is owned by the rank only if
	    // the rank number is the highest of all ranks containing this dof.

	    bool insert = true;
	    int highestRank = mpiRank;
	    for (std::set<int>::iterator itpart2 = it->second.begin();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1717
		 itpart2 != it->second.end(); ++itpart2) {
1718
1719
1720
1721
1722
1723
1724
1725
	      if (*itpart2 > mpiRank)
		insert = false;

	      if (*itpart2 > highestRank)
		highestRank = *itpart2;
	    }

	    if (insert)
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1726
	      rankOwnedDofs.push_back(it->first);
1727

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1728
	    boundaryDofs[it->first] = highestRank;
1729
1730
	  }

1731
	  isInRank = true;
1732
1733
	  break;
	}
1734

1735
      }
1736

1737
      if (!isInRank)
1738
	vertexDof.erase(it->first);
1739
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1740

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1741
1742
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
    sort(rankOwnedDofs.begin(), rankOwnedDofs.end(), cmpDofsByValue);
1743
1744
1745
  }


1746
1747
1748
1749
  void ParallelDomainBase::createPeriodicMap()
  {
    FUNCNAME("ParallelDomainBase::createPeriodicMap()");

1750
    // Clear all periodic dof mappings calculated before. We do it from scratch.
1751
1752
    periodicDof.clear();

1753
    MPI::Request request[min(static_cast<int>(periodicBoundary.boundary.size() * 2), 4)];
1754
1755
1756
    int requestCounter = 0;
    std::vector<int*> sendBuffers, recvBuffers;

1757
1758
1759
    // === Each rank traverse its periodic boundaries and sends the dof indices ===
    // === to the rank "on the other side" of the periodic boundary.            ===

1760
1761
1762
1763
1764
1765
1766
1767
    for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
	 it != periodicBoundary.boundary.end(); ++it) {

      TEST_EXIT_DBG(it->first != mpiRank)
	("Periodic interior boundary within the rank itself is not possible!\n");

      DofContainer dofs;

1768
1769
      // Create dof indices on the boundary. 

1770
1771
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
1772
1773
	boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs, true);
	boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, dofs);
1774
      }
1775

1776
1777
1778
1779
1780
      // Send the global indices to the rank on the other side.

      int *sendbuf = new int[dofs.size()];
      for (int i = 0; i < static_cast<int>(dofs.size()); i++)
	sendbuf[i] = mapLocalGlobalDOFs[*(dofs[i])];      
1781
1782
1783
1784
1785
1786
      
      request[requestCounter++] = 
	mpiComm.Isend(sendbuf, dofs.size(), MPI_INT, it->first, 0);
      
      sendBuffers.push_back(sendbuf);

1787
1788
      // Receive from this rank the same number of dofs.

1789
1790
1791
1792
1793
1794
1795
1796
      int *recvbuf = new int[dofs.size()];
      request[requestCounter++] = 
	mpiComm.Irecv(recvbuf, dofs.size(), MPI_INT, it->first, 0);

      recvBuffers.push_back(recvbuf);
    }


1797
    MPI::Request::Waitall(requestCounter, request);    
1798

1799
1800
1801
1802
1803
    // === The rank has received the dofs from the rank on the other side of ===
    // === the boundary. Now it can use them to create the mapping between   ===
    // === the periodic dofs in this rank and the corresponding periodic     ===
    // === dofs from the other ranks.                                        ===

1804
1805
1806

    std::map<DegreeOfFreedom, std::set<int> > dofFromRank;

1807
1808
1809
1810
1811
    int i = 0;
    for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
	 it != periodicBoundary.boundary.end(); ++it) {
      DofContainer dofs;
      
1812
      // Create the dofs on the boundary in inverse order.
1813
1814
1815
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	DofContainer tmpdofs;
1816
1817
	boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, tmpdofs);
	boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, tmpdofs, true);
1818
1819
	for (int j = static_cast<int>(tmpdofs.size()) - 1; j >= 0; j--)
	  dofs.push_back(tmpdofs[j]);	
1820
1821
      }

1822
      // Added the received dofs to the mapping.
1823
1824
1825
1826
1827
      for (int j = 0; j < static_cast<int>(dofs.size()); j++) {
	int globalDofIndex = mapLocalGlobalDOFs[*(dofs[j])];
	periodicDof[globalDofIndex].insert(recvBuffers[i][j]);      
	dofFromRank[globalDofIndex].insert(it->first);
      }
1828
1829
1830
1831
1832
      
      delete [] sendBuffers[i];
      delete [] recvBuffers[i];
      i++;
    }
1833

1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
    sendBuffers.clear();
    recvBuffers.clear();

    TEST_EXIT_DBG(mesh->getDim() == 2)
      ("Periodic boundary corner problem must be generalized to 3d!\n");

    requestCounter = 0;
    for (std::map<DegreeOfFreedom, std::set<int> >::iterator it = dofFromRank.begin();
	 it != dofFromRank.end(); ++it) {
      if (it->second.size() == 2) {
	TEST_EXIT_DBG(periodicDof[it->first].size() == 2)("Missing periodic dof!\n");
	
	int *sendbuf = new int[2];
	sendbuf[0] = *(periodicDof[it->first].begin());
	sendbuf[1] = *(++(periodicDof[it->first].begin()));
	
	request[requestCounter++] = 
	  mpiComm.Isend(sendbuf, 2, MPI_INT, *(it->second.begin()), 0);
	request[requestCounter++] = 
	  mpiComm.Isend(sendbuf, 2, MPI_INT, *(++(it->second.begin())), 0);
	
	sendBuffers.push_back(sendbuf);

	int *recvbuf1 = new int[2];
	int *recvbuf2 = new int[2];

	request[requestCounter++] = 
	  mpiComm.Irecv(recvbuf1, 2, MPI_INT, *(it->second.begin()), 0);
	request[requestCounter++] = 
	  mpiComm.Irecv(recvbuf2, 2, MPI_INT, *(++(it->second.begin())), 0);

	recvBuffers.push_back(recvbuf1);
	recvBuffers.push_back(recvbuf2);
      }
    }

    MPI::Request::Waitall(requestCounter, request);

    i = 0;
    for (std::map<DegreeOfFreedom, std::set<int> >::iterator it = dofFromRank.begin();
	 it != dofFromRank.end(); ++it) {
      if (it->second.size() == 2) {
	for (int k = 0; k < 2; k++)
	  for (int j = 0; j < 2; j++)
	    if (recvBuffers[i + k][j] != it->first)
	      periodicDof[it->first].insert(recvBuffers[i + k][j]);

	i++;
      }
    }
1884
1885
1886
  }


1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
  Flag ParallelDomainBase::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    FUNCNAME("ParallelDomainBase::oneIteration()");

    Flag flag = dynamic_cast<StandardProblemIteration*>(iterationIF)->
      buildAndAdapt(adaptInfo, toDo);

    if (toDo.isSet(SOLVE))
      solve();

    if (toDo.isSet(SOLVE_RHS))
      ERROR_EXIT("Not yet implemented!\n");

    if (toDo.isSet(ESTIMATE))
      iterationIF->getProblem()->estimate(adaptInfo);

    return flag;
  }


  void ParallelDomainBase::serialize(std::ostream &out)
  {
    SerUtil::serialize(out, elemWeights);
1910
    SerUtil::serialize(out, initialPartitionMesh);
1911
1912
    SerUtil::serialize(out, partitionVec);
    SerUtil::serialize(out, oldPartitionVec);
1913
    SerUtil::serialize(out, nRankDofs);
1914
1915
1916
1917
    
    myIntBoundary.serialize(out);
    otherIntBoundary.serialize(out);

1918
1919
    serialize(out, sendDofs);
    serialize(out, recvDofs);
1920
1921
1922
1923
1924

    SerUtil::serialize(out, mapLocalGlobalDOFs);
    SerUtil::serialize(out, mapLocalToDofIndex);
    SerUtil::serialize(out, isRankDof);

1925
    serialize(out, vertexDof);
1926
    serialize(out, periodicDof);
1927

1928
1929
1930
    SerUtil::serialize(out, rstart);
    SerUtil::serialize(out, nRankRows);
    SerUtil::serialize(out, nOverallRows);
1931
1932
1933
1934
1935
  }


  void ParallelDomainBase::deserialize(std::istream &in)
  {