ParallelDomainBase.cc 75.3 KB
Newer Older
1001
1002
	    int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()];

1003
1004
1005
1006
1007
1008
1009
1010
1011

	    // === 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.                                         ===

1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
	    if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(i))) {	      
	      // 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;	      
	    }
1023
1024
1025
1026
1027
1028
 	  }
	}
      }

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

1030
1031
1032
1033

    // === 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
1034

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

1037
    // First we communicate myInt/otherIntBoundary, and than the periodic boundaries.
1038
1039
1040
    MPI::Request request[max(myIntBoundary.boundary.size() + 
			     otherIntBoundary.boundary.size(),
			     periodicBoundary.boundary.size())];
Thomas Witkowski's avatar
Thomas Witkowski committed
1041
1042
    int requestCounter = 0;

1043

1044
1045
1046
    // === The owner of the boundaries send their atomic boundary order to the ===
    // === neighbouring ranks.                                                 ===

Thomas Witkowski's avatar
Thomas Witkowski committed
1047
1048
    for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin();
	 rankIt != myIntBoundary.boundary.end(); ++rankIt) {    
Thomas Witkowski's avatar
Thomas Witkowski committed
1049
      int nSendInt = rankIt->second.size();
1050
1051
1052
1053
1054
      int* buffer = new int[nSendInt * 2];
      for (int i = 0; i < nSendInt; i++) {
	buffer[i * 2] = (rankIt->second)[i].rankObject.elIndex;
	buffer[i * 2 + 1] = (rankIt->second)[i].rankObject.ithObjAtBoundary;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
1055
1056
      sendBuffers.push_back(buffer);
      
Thomas Witkowski's avatar
Thomas Witkowski committed
1057
      request[requestCounter++] =
1058
	mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1059
1060
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1061
1062
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1063
      int nRecvInt = rankIt->second.size();
1064
      int *buffer = new int[nRecvInt * 2];
Thomas Witkowski's avatar
Thomas Witkowski committed
1065
      recvBuffers.push_back(buffer);
Thomas Witkowski's avatar
Thomas Witkowski committed
1066

Thomas Witkowski's avatar
Thomas Witkowski committed
1067
      request[requestCounter++] = 
1068
	mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1069
1070
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1071
1072
1073

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

1074
1075
1076
1077
    
    // === 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
1078
1079

    int i = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1080
1081
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1082
1083
1084
1085
1086
1087
1088

      // === 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.
1089
1090
	if ((rankIt->second)[j].neighbourObject.elIndex != recvBuffers[i][j * 2] ||
	    (rankIt->second)[j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1091
1092
	  int k = j + 1;
	  for (; k < static_cast<int>(rankIt->second.size()); k++)
1093
1094
	    if ((rankIt->second)[k].neighbourObject.elIndex == recvBuffers[i][j * 2] && 
		(rankIt->second)[k].neighbourObject.ithObjAtBoundary == recvBuffers[i][j * 2 + 1])
Thomas Witkowski's avatar
Thomas Witkowski committed
1095
1096
1097
	      break;

	  // The element must always be found, because the list is just in another order.
Thomas Witkowski's avatar
Thomas Witkowski committed
1098
	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
1099
	    ("Should never happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
1100
1101
1102
1103
1104
1105
1106
1107

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

1108
1109
      delete [] recvBuffers[i];
      i++;
Thomas Witkowski's avatar
Thomas Witkowski committed
1110
1111
1112
1113
    }

    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
      delete [] sendBuffers[i];
1114
1115
1116
1117
1118

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

1119
    // === Do the same for the periodic boundaries. ===
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185

    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++) {
	  buffer[i * 2] = (rankIt->second)[i].rankObject.elIndex;
	  buffer[i * 2 + 1] = (rankIt->second)[i].rankObject.ithObjAtBoundary;
	}
	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++) 
	  if (periodicBoundary.boundary[rankIt->first][j].neighbourObject.elIndex != recvBuffers[i][j * 2] ||
	      periodicBoundary.boundary[rankIt->first][j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) {
	    int k = j + 1;

	    for (; k < static_cast<int>(rankIt->second.size()); k++)
	      if (periodicBoundary.boundary[rankIt->first][k].neighbourObject.elIndex == recvBuffers[i][j * 2] && 
		  periodicBoundary.boundary[rankIt->first][k].neighbourObject.ithObjAtBoundary == recvBuffers[i][j * 2 + 1])
		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]; 
1186
1187
1188
  }


1189
  void ParallelDomainBase::removeMacroElements()
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
  {
    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
1202
    mesh->removeMacroElements(macrosToRemove, feSpace);
1203
1204
1205
  }


1206
  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDofs,
1207
						      int& nRankDofs, 
Thomas Witkowski's avatar
Thomas Witkowski committed
1208
						      int& nOverallDOFs)
1209
  {
1210
    FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
1211
1212

    // === Get rank information about DOFs. ===
1213
1214
1215

    // 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
1216
    DofContainer rankAllDofs;
1217
    DofToRank boundaryDofs;
1218

1219
    createDofMemberInfo(partitionDOFs, rankDofs, rankAllDofs, boundaryDofs, vertexDof);
1220

1221
    nRankDofs = rankDofs.size();
1222
1223
    nOverallDOFs = partitionDOFs.size();

1224

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

1227
    rstart = 0;
1228
1229
    mpiComm.Scan(&nRankDofs, &rstart, 1, MPI_INT, MPI_SUM);
    rstart -= nRankDofs;
1230

1231
1232
1233
1234
    
    // === 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
1235

Thomas Witkowski's avatar
Thomas Witkowski committed
1236

1237
    // Do not change the indices now, but create a new indexing and store it here.
1238
1239
    DofIndexMap rankDofsNewLocalIndex;
    isRankDof.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1240
1241
1242
    int i = 0;
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
1243
1244
1245
      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
1246
      isRankDof[i] = true;
1247
      i++;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1248
1249
    }

1250
    // === Create for all rank owned dofs a new global indexing. ===
1251

1252
1253
    // Stores for dofs in rank a new global index.
    DofIndexMap rankDofsNewGlobalIndex;
1254
1255
1256
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1257
    i = 0;
1258
1259
    for (DofContainer::iterator dofIt = rankDofs.begin();
	 dofIt != rankDofs.end(); ++dofIt) {
1260
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
1261
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1262
1263
1264
      i++;
    }

1265
1266
1267
1268
1269
    // === 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.
1270
    std::map<int, DofIndexMap > sendNewDofs;
1271

1272
1273
1274
1275
    // 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;
1276

1277
    for (DofToRank::iterator it = boundaryDofs.begin(); it != boundaryDofs.end(); ++it) {
1278
1279
1280
1281
1282
1283
1284
1285

      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
1286
	  if (*itRanks != mpiRank) {
1287
	    TEST_EXIT_DBG(rankDofsNewGlobalIndex.count(it->first) == 1)
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1288
1289
	      ("DOF Key not found!\n");

1290
	    sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first];
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1291
	  }
1292
1293
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1294
1295
	// 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.
1296
1297
1298
1299
	if (recvNewDofs.find(it->second) == recvNewDofs.end()) 
	  recvNewDofs[it->second] = 1;
	else
	  recvNewDofs[it->second]++;
1300
1301
1302
      }
    }

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1303

1304
    // === Send and receive the dof indices at boundary. ===
1305

Thomas Witkowski's avatar
Thomas Witkowski committed
1306
    std::vector<int*> sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size());
1307
1308
1309

    sendDofs.clear();
    recvDofs.clear();
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325

#if 0
    StlMpi stlMpi(mpiComm);
    stlMpi.prepareCommunication(false);
    std::map<

    for (std::map<int, DofIndexMap >::iterator sendIt = sendNewDofs.begin();
	 sendIt != sendNewDofs.end(); ++sendIt, i++)
      stlMpi.send(it->second, it->first);

    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
	 recvIt != recvNewDofs.end(); ++recvIt, i++)
      stlMpi.recv(it->first, it->second * 2); 

    stlMpi.startCommunication();
#endif
1326
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1327
1328
1329
    MPI::Request request[sendNewDofs.size() + recvNewDofs.size()];
    int requestCounter = 0;

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1330
    i = 0;
1331
    for (std::map<int, DofIndexMap>::iterator sendIt = sendNewDofs.begin();
1332
	 sendIt != sendNewDofs.end(); ++sendIt, i++) {
1333
1334
      int nSendDofs = sendIt->second.size() * 2;
      sendBuffers[i] = new int[nSendDofs];
1335

1336
      int c = 0;
1337
1338
      for (DofIndexMap::iterator dofIt = sendIt->second.begin();
	   dofIt != sendIt->second.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1339
	sendBuffers[i][c++] = *(dofIt->first);
1340
	sendBuffers[i][c++] = dofIt->second;
1341

1342
	sendDofs[sendIt->first].push_back(dofIt->first);
1343
1344
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
1345
1346
      request[requestCounter++] =
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
1347
1348
1349
    }

    i = 0;
1350
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
1351
	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
1352
1353
      int nRecvDOFs = recvIt->second * 2;
      recvBuffers[i] = new int[nRecvDOFs];
1354

Thomas Witkowski's avatar
Thomas Witkowski committed
1355
1356
      request[requestCounter++] =
	mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0);
1357
1358
1359
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
1360
    MPI::Request::Waitall(requestCounter, request);
1361
1362

    
1363
    // === Delete send buffers. ===
1364

Thomas Witkowski's avatar
Thomas Witkowski committed
1365
1366
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
1367
1368


1369
    // === Change dof indices at boundary from other ranks. ===
1370

Thomas Witkowski's avatar
Thomas Witkowski committed
1371
1372
    // 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
1373
    // 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
1374
1375
1376
1377
1378
    // 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;
1379
    for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
1380
	 ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
1381
1382
      dofChanged[dofIt->first] = false;

1383
    i = 0;
1384
    for (std::map<int, int>::iterator recvIt = recvNewDofs.begin();
1385
1386
1387
	 recvIt != recvNewDofs.end();
	 ++recvIt, i++) {

1388
      for (int j = 0; j < recvIt->second; j++) {
1389

1390
1391
	DegreeOfFreedom oldDof = recvBuffers[i][j * 2];
	DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1];
1392

1393
	bool found = false;
1394

1395
	// Iterate over all boundary dofs to find the dof which index we have to change.
Thomas Witkowski's avatar
Thomas Witkowski committed
1396

1397
1398
	for (DofToRank::iterator dofIt = boundaryDofs.begin(); 
	     dofIt != boundaryDofs.end(); ++dofIt) {
1399

Thomas Witkowski's avatar
Thomas Witkowski committed
1400
1401
1402
	  if (*(dofIt->first) == oldDof && !dofChanged[dofIt->first]) {
	    dofChanged[dofIt->first] = true;

1403
	    recvDofs[recvIt->first].push_back(dofIt->first);
1404
	    rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof;
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1405
	    isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
1406

1407
	    found = true;
1408
1409
1410
	    break;
	  }
	}
1411

Thomas Witkowski's avatar
Thomas Witkowski committed
1412
	TEST_EXIT_DBG(found)("Should not happen!\n");
1413
1414
1415
1416
      }

      delete [] recvBuffers[i];
    }
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1417

1418

1419
    // === Create now the local to global index and local to dof index mappings.  ===
1420

1421
1422
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
1423
1424
1425
  }


1426
  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDofs, int& nOverallDOFs)
1427
  {
1428
    FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
Thomas Witkowski's avatar
Thomas Witkowski committed
1429

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

1432
    // === Get all DOFs in ranks partition. ===
1433

Thomas Witkowski's avatar
Thomas Witkowski committed
1434
    ElementDofIterator elDofIt(feSpace);
1435
    DofSet rankDofSet;
1436
    
1437
1438
1439
    // 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).
1440
    DofToBool oldVertexDof = vertexDof;
1441
1442
    vertexDof.clear();

1443
1444
1445
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1446
      Element *element = elInfo->getElement();     
Thomas Witkowski's avatar
Thomas Witkowski committed
1447
1448
      elDofIt.reset(element);
      do {
1449
1450
1451
1452
1453
1454
	rankDofSet.insert(elDofIt.getDofPtr());

	if (elDofIt.getCurrentPos() == 0) 
	  vertexDof[elDofIt.getDofPtr()] = true;
	else
	  vertexDof[elDofIt.getDofPtr()] = false;
Thomas Witkowski's avatar
Thomas Witkowski committed
1455
      } while(elDofIt.next());
1456
1457
1458
1459

      elInfo = stack.traverseNext(elInfo);
    }

1460
    DofContainer rankAllDofs;
1461
    for (DofSet::iterator dofIt = rankDofSet.begin(); dofIt != rankDofSet.end(); ++dofIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
1462
      rankAllDofs.push_back(*dofIt);    
1463
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
1464
    DofContainer rankDofs = rankAllDofs;
1465
1466


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

1470
1471
1472
    RankToDofContainer oldSendDofs = sendDofs;
    RankToDofContainer oldRecvDofs = recvDofs;

1473
1474
    sendDofs.clear();
    recvDofs.clear();
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498

    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);    
	}
      }
    }

1499
    
1500
1501
1502
1503
    for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin();
	 it != myIntBoundary.boundary.end(); ++it) {    
      
      DofContainer &dofsToSend = sendDofs[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
1504

Thomas Witkowski's avatar
Thomas Witkowski committed
1505
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1506
	   boundIt != it->second.end(); ++boundIt) {
1507
	DofContainer dofs;	
1508
1509
1510
1511
	addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		      dofs);	
  	addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 
		    dofs);
1512
1513
1514
1515
	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");

1516
	  dofsToSend.push_back(dofs[i]);
1517
	}
1518
      }
1519
    }
1520
   
1521

Thomas Witkowski's avatar
Thomas Witkowski committed
1522
1523
1524
    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end(); ++it) {

1525
      DofContainer &dofsToRecv = recvDofs[it->first];
1526

1527
1528
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
1529

1530
	DofContainer dofs;
1531
1532
1533
1534
	addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 
		    dofs);
	addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
		      dofs);
1535
1536

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

1540
1541
1542
	  DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), dofs[i]);
	  if (eraseIt != rankDofs.end())
	    rankDofs.erase(eraseIt);	  
1543

1544
	  dofsToRecv.push_back(dofs[i]);
1545
	}
1546

Thomas Witkowski's avatar
Thomas Witkowski committed
1547
1548
1549
      }
    }

1550
    nRankDofs = rankDofs.size();
1551
1552
1553

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

1554
1555
    mpiComm.Scan(&nRankDofs, &rstart, 1, MPI_INT, MPI_SUM);
    rstart -= nRankDofs;
1556
1557
1558
1559


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

1560
    mpiComm.Allreduce(&nRankDofs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
1561
1562


1563
1564
1565
    // Do not change the indices now, but create a new indexing a store it here.
    DofIndexMap rankDofsNewLocalIndex;
    isRankDof.clear();
1566
    int i = 0;
1567
1568
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
1569

1570
1571
1572
      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
1573
      isRankDof[i] = true;
1574
      i++;
1575
1576
    }

1577
    // Stores for all rank owned dofs a new global index.
1578
    DofIndexMap rankDofsNewGlobalIndex;
1579
1580
1581
1582
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

    i = 0;
1583
1584
    for (DofContainer::iterator dofIt = rankDofs.begin();
	 dofIt != rankDofs.end(); ++dofIt) {
1585
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
1586
1587
1588
1589
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
      i++;
    }

1590
1591
    // === Send new DOF indices. ===

1592
1593
    std::vector<int*> sendBuffers(sendDofs.size());
    std::vector<int*> recvBuffers(recvDofs.size());
1594

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

1598
    i = 0;
1599
1600
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt, i++) {
1601
1602
1603
      int nSendDofs = sendIt->second.size();
      sendBuffers[i] = new int[nSendDofs];
      int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1604
      for (DofContainer::iterator dofIt = sendIt->second.begin();
1605
	   dofIt != sendIt->second.end(); ++dofIt)
1606
	sendBuffers[i][c++] = rankDofsNewGlobalIndex[*dofIt];
1607

Thomas Witkowski's avatar
Thomas Witkowski committed
1608
      request[requestCounter++] = 
1609
 	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
1610
1611
1612
    }

    i = 0;
1613
1614
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt, i++) {
1615
1616
      int nRecvDofs = recvIt->second.size();
      recvBuffers[i] = new int[nRecvDofs];
1617

Thomas Witkowski's avatar
Thomas Witkowski committed
1618
      request[requestCounter++] = 
1619
 	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
1620
1621
    }

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

1624
1625
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
1626
1627

    i = 0;
1628
1629
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt) {      
1630
      int j = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1631
      for (DofContainer::iterator dofIt = recvIt->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1632
	   dofIt != recvIt->second.end(); ++dofIt) {
1633

1634
	rankDofsNewGlobalIndex[*dofIt] = recvBuffers[i][j];     
1635
	isRankDof[rankDofsNewLocalIndex[*dofIt]] = false;
1636
1637
1638
1639
	j++;
      }

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

1642
    // === Create now the local to global index and local to dof index mappings.  ===
1643

1644
1645
1646
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
  }
1647

1648
1649
1650
1651
  void ParallelDomainBase::createLocalMappings(DofIndexMap &rankDofsNewLocalIndex,
					       DofIndexMap &rankOwnedDofsNewLocalIndex,
					       DofIndexMap &rankDofsNewGlobalIndex)
  {
1652
    mapLocalGlobalDOFs.clear();
1653
    mapLocalToDofIndex.clear();
1654
1655
1656
1657

    for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin();
	 dofIt != rankDofsNewLocalIndex.end(); ++dofIt) {
      DegreeOfFreedom localDof = dofIt->second;
1658
      DegreeOfFreedom globalDof = rankDofsNewGlobalIndex[dofIt->first];
1659
1660
1661
1662
1663
1664
1665
1666

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

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

1669
1670
  void ParallelDomainBase::addVertexDofs(Element *el, int ithEdge, DofContainer& dofs,
					 bool parentVertices)
1671
  {
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
    FUNCNAME("ParallelDomainBase::addVertexDofs()");

    TEST_EXIT_DBG(mesh->getDim() == 2)("Not yet implemented for 3D!\n");

    if (parentVertices) {
      switch (ithEdge) {
      case 0:
	dofs.push_back(el->getDOF(1));
	break;
      case 1:
	dofs.push_back(el->getDOF(0));
	break;
      case 2:
	dofs.push_back(el->getDOF(0));
	break;
    default:
      ERROR_EXIT("Should never happen!\n");
      }
    }
1691

Thomas Witkowski's avatar
Thomas Witkowski committed
1692
1693
1694
    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
1695
	addVertexDofs(el->getSecondChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1696
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
1697
	addVertexDofs(el->getSecondChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1698
1699
1700
1701
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
1702
	addVertexDofs(el->getFirstChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1703
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
1704
	addVertexDofs(el->getFirstChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1705
1706
1707
1708
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
1709
	addVertexDofs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1710
	dofs.push_back(el->getFirstChild()->getDOF(2));
1711
	addVertexDofs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1712
1713
1714
1715
1716
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732

    if (parentVertices) {
      switch (ithEdge) {
      case 0:
	dofs.push_back(el->getDOF(2));
	break;
      case 1:
	dofs.push_back(el->getDOF(2));
	break;
      case 2:
	dofs.push_back(el->getDOF(1));
	break;
    default:
      ERROR_EXIT("Should never happen!\n");
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1733
1734
1735
  }


1736
  void ParallelDomainBase::addEdgeDofs(Element *el, int ithEdge, DofContainer& dofs)
Thomas Witkowski's avatar
Thomas Witkowski committed
1737
  {
1738
1739
1740
    FUNCNAME("ParallelDomainBase::addEdgeDofs()");

    TEST_EXIT_DBG(mesh->getDim() == 2)("Not yet implemented for 3D!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
1741
1742
1743
1744
1745

    bool addThisEdge = false;

    switch (ithEdge) {
    case 0:
1746
      if (el->getSecondChild())
1747
	addEdgeDofs(el->getSecondChild(), 2, dofs);
1748
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
1749
1750
1751
1752
	addThisEdge = true;

      break;
    case 1:
1753
      if (el->getFirstChild())
1754
	addEdgeDofs(el->getFirstChild(), 2, dofs);
1755
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
1756
1757
1758
1759
1760
	addThisEdge = true;

      break;
    case 2:
      if (el->getFirstChild()) {
1761
1762
	addEdgeDofs(el->getFirstChild(), 0, dofs);
	addEdgeDofs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
      } else {
	addThisEdge = true;
      } 

      break;
    default:
      ERROR_EXIT("Should never happen!\n");
    }

    if (addThisEdge) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1773
      ElementDofIterator elDofIter(feSpace, true);
Thomas Witkowski's avatar
Thomas Witkowski committed
1774
1775
      elDofIter.reset(el);
      do {
1776
1777
1778
	if (elDofIter.getCurrentPos() == 1 && 
	    elDofIter.getCurrentElementPos() == ithEdge)
	  dofs.push_back(elDofIter.getDofPtr());	
Thomas Witkowski's avatar
Thomas Witkowski committed
1779
1780
1781
1782
1783
      } while(elDofIter.next());      
    }
  }


1784
  void ParallelDomainBase::createDofMemberInfo(DofToPartitions& partitionDofs,
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1785
1786
					       DofContainer& rankOwnedDofs,
					       DofContainer& rankAllDofs,
1787
1788
					       DofToRank& boundaryDofs,
					       DofToBool& vertexDof)
1789
  {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1790
1791
1792
1793
    partitionDofs.clear();
    rankOwnedDofs.clear();
    rankAllDofs.clear();
    boundaryDofs.clear();
1794
    vertexDof.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1795

1796
    // === Determine to each dof the set of partitions the dof belongs to. ===
1797

Thomas Witkowski's avatar
Thomas Witkowski committed
1798
    ElementDofIterator elDofIt(feSpace);
1799

1800
1801
1802
1803
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *element = elInfo->getElement();
1804
1805
1806
      elDofIt.reset(element);
      do {
	// Determine to each dof the partition(s) it corresponds to.
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1807
	partitionDofs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
1808
1809
1810
1811
1812

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

1815
1816
1817
      elInfo = stack.traverseNext(elInfo);
    }

1818

1819
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
1820

1821
    // iterate over all DOFs
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1822
1823
    for (DofToPartitions::iterator it = partitionDofs.begin();
	 it != partitionDofs.end(); ++it) {
1824

1825
1826
      bool isInRank = false;

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

1831
	if (*itpart1 == mpiRank) {	  
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1832
1833
	  rankAllDofs.push_back(it->first);

1834
	  if (it->second.size() == 1) {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1835
	    rankOwnedDofs.push_back(it->first);
1836
1837
1838
1839
1840
1841
1842
	  } 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
1843
		 itpart2 != it->second.end(); ++itpart2) {
1844
1845
1846
1847
1848
1849
1850
1851
	      if (*itpart2 > mpiRank)
		insert = false;

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

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

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1854
	    boundaryDofs[it->first] = highestRank;
1855
1856
	  }

1857
	  isInRank = true;
1858
1859
	  break;
	}
1860

1861
      }
1862

1863
      if (!isInRank)
1864
	vertexDof.erase(it->first);
1865
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1866

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1867
1868
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
    sort(rankOwnedDofs.begin(), rankOwnedDofs.end(), cmpDofsByValue);
1869
1870
1871
  }


1872
1873
1874
1875
  void ParallelDomainBase::createPeriodicMap()
  {
    FUNCNAME("ParallelDomainBase::createPeriodicMap()");

1876
    // Clear all periodic dof mappings calculated before. We do it from scratch.
1877
1878
    periodicDof.clear();

1879
1880
1881
1882
    MPI::Request request[periodicBoundary.boundary.size() * 2];
    int requestCounter = 0;
    std::vector<int*> sendBuffers, recvBuffers;

1883
1884
1885
    // === Each rank traverse its periodic boundaries and sends the dof indices ===
    // === to the rank "on the other side" of the periodic boundary.            ===

1886
1887
1888
1889
1890
1891
1892
1893
    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;

1894
1895
      // Create dof indices on the boundary. 

1896
1897
1898
1899
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 
		      dofs, true);
1900
	addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 
1901
 		    dofs);
1902
      }
1903

1904
1905
1906
1907
1908
      // 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])];      
1909
1910
1911
1912
1913
1914
      
      request[requestCounter++] = 
	mpiComm.Isend(sendbuf, dofs.size(), MPI_INT, it->first, 0);
      
      sendBuffers.push_back(sendbuf);

1915
1916
      // Receive from this rank the same number of dofs.

1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
      int *recvbuf = new int[dofs.size()];
      request[requestCounter++] = 
	mpiComm.Irecv(recvbuf, dofs.size(), MPI_INT, it->first, 0);

      recvBuffers.push_back(recvbuf);
    }


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


1928
1929
1930
1931
1932
    // === 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.                                        ===

1933
1934
1935
1936
1937
    int i = 0;
    for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
	 it != periodicBoundary.boundary.end(); ++it) {
      DofContainer dofs;
      
1938
      // Create the dofs on the boundary in inverse order.
1939
1940
1941
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	DofContainer tmpdofs;
1942
1943
 	addEdgeDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 
 		    tmpdofs);
1944
1945
	addVertexDofs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 
		      tmpdofs, true);
1946
1947
	for (int j = static_cast<int>(tmpdofs.size()) - 1; j >= 0; j--)
	  dofs.push_back(tmpdofs[j]);	
1948
1949
      }

1950
1951
1952
      // Added the received dofs to the mapping.
      for (int j = 0; j < static_cast<int>(dofs.size()); j++)
	periodicDof[mapLocalGlobalDOFs[*(dofs[j])]].insert(recvBuffers[i][j]);      
1953
1954
1955
1956