ParallelDomainBase.cc 50.3 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1001
1002
	 it != myIntBoundary.boundary.end(); ++it) {

Thomas Witkowski's avatar
Thomas Witkowski committed
1003
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1004
1005
	   boundIt != it->second.end(); ++boundIt) {

1006
1007
	DofContainer dofs;
	DofContainer &dofsToSend = sendDofs[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
1008
1009
1010

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
1011
1012
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1013
1014
	  break;
	case 1:
1015
1016
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1017
1018
	  break;
	case 2:
1019
1020
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
Thomas Witkowski's avatar
Thomas Witkowski committed
1021
1022
1023
1024
1025
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

1026
1027
1028
	for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt)
	  if (find(dofsToSend.begin(), dofsToSend.end(), *dofIt) == dofsToSend.end())
	    dofsToSend.push_back(*dofIt);	
1029

1030
1031
1032
1033
1034
	dofs.clear();
	addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
			 dofs);	
  	addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
  		       dofs);
1035
	for (int i = 0; i < static_cast<int>(dofs.size()); i++)
1036
	  dofsToSend.push_back(dofs[i]);
1037
      }
1038
1039
1040
    }

    
1041

Thomas Witkowski's avatar
Thomas Witkowski committed
1042
1043
1044
    for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin();
	 it != otherIntBoundary.boundary.end(); ++it) {

1045
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1046
1047
	   boundIt != it->second.end(); ++boundIt) {

1048
1049
	DofContainer dofs;
	DofContainer &dofsToRecv = recvDofs[it->first];
1050
1051
1052

	switch (boundIt->rankObject.ithObjAtBoundary) {
	case 0:
1053
1054
1055
1056
1057
1058
1059
	  if (boundIt->neighbourObject.ithObjAtBoundary == 0) {
	    dofs.push_back(boundIt->rankObject.el->getDOF(2));
	    dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  } else {
	    dofs.push_back(boundIt->rankObject.el->getDOF(1));
	    dofs.push_back(boundIt->rankObject.el->getDOF(2));
	  }
1060
1061
	  break;
	case 1:
1062
1063
1064
1065
1066
1067
1068
	  if (boundIt->neighbourObject.ithObjAtBoundary == 0) {
	    dofs.push_back(boundIt->rankObject.el->getDOF(0));
	    dofs.push_back(boundIt->rankObject.el->getDOF(2));
	  } else {
	    dofs.push_back(boundIt->rankObject.el->getDOF(2));
	    dofs.push_back(boundIt->rankObject.el->getDOF(0));
	  }
1069
1070
	  break;
	case 2:
1071
1072
	  dofs.push_back(boundIt->rankObject.el->getDOF(1));
	  dofs.push_back(boundIt->rankObject.el->getDOF(0));
1073
1074
1075
1076
1077
	  break;
	default:
	  ERROR_EXIT("Should never happen!\n");
	}

1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
	for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt) {
	  DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), *dofIt);
	  if (eraseIt != rankDOFs.end())
	    rankDOFs.erase(eraseIt);
	  if (find(dofsToRecv.begin(), dofsToRecv.end(), *dofIt) == dofsToRecv.end())
	    dofsToRecv.push_back(*dofIt);
	}

	dofs.clear();
	addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 
		       dofs);
	addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
			 dofs);

	for (int i = static_cast<int>(dofs.size()) - 1; i >= 0; i--) {
	  TEST_EXIT_DBG(find(rankDOFs.begin(), rankDOFs.end(), dofs[i]) != rankDOFs.end())
1094
1095
	    ("Should never happen!\n");

1096
	  DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), dofs[i]);
1097
1098
1099
	  if (eraseIt != rankDOFs.end())
	    rankDOFs.erase(eraseIt);

1100
	  dofsToRecv.push_back(dofs[i]);
1101
	}
1102

Thomas Witkowski's avatar
Thomas Witkowski committed
1103
1104
1105
      }
    }

1106
1107
1108
1109
    nRankDOFs = rankDOFs.size();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1110
    mpiComm.Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM);
1111
1112
1113
1114
1115
    rstart -= nRankDOFs;


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

Thomas Witkowski's avatar
Thomas Witkowski committed
1116
    mpiComm.Allreduce(&nRankDOFs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
1117
1118


1119
1120
1121
    // Do not change the indices now, but create a new indexing a store it here.
    DofIndexMap rankDofsNewLocalIndex;
    isRankDof.clear();
1122
    int i = 0;
1123
1124
    for (DofContainer::iterator dofIt = rankAllDofs.begin();
	 dofIt != rankAllDofs.end(); ++dofIt) {
1125

1126
1127
1128
      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
1129
      isRankDof[i] = true;
1130
      i++;
1131
1132
    }

1133
    // Stores for all rank owned dofs a new global index.
1134
    DofIndexMap rankDofsNewGlobalIndex;
1135
1136
1137
1138
1139
1140
    // Stores for all rank owned dofs a continues local index.
    DofIndexMap rankOwnedDofsNewLocalIndex;

    i = 0;
    for (DofContainer::iterator dofIt = rankDOFs.begin();
	 dofIt != rankDOFs.end(); ++dofIt) {
1141
      rankDofsNewGlobalIndex[*dofIt] = i + rstart;
1142
1143
1144
1145
1146
      rankOwnedDofsNewLocalIndex[*dofIt] = i;
      i++;
    }


1147
1148
    // === Send new DOF indices. ===

1149
1150
    std::vector<int*> sendBuffers(sendDofs.size());
    std::vector<int*> recvBuffers(recvDofs.size());
1151

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

1155
    i = 0;
1156
1157
    for (RankToDofContainer::iterator sendIt = sendDofs.begin();
	 sendIt != sendDofs.end(); ++sendIt, i++) {
1158
1159
1160
      int nSendDofs = sendIt->second.size();
      sendBuffers[i] = new int[nSendDofs];
      int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1161
      for (DofContainer::iterator dofIt = sendIt->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1162
	   dofIt != sendIt->second.end(); ++dofIt)
1163
	sendBuffers[i][c++] = rankDofsNewGlobalIndex[*dofIt];
1164

Thomas Witkowski's avatar
Thomas Witkowski committed
1165
1166
      request[requestCounter++] = 
	mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
1167
1168
1169
    }

    i = 0;
1170
1171
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt, i++) {
1172
1173
1174
      int nRecvDofs = recvIt->second.size();
      recvBuffers[i] = new int[nRecvDofs];
	
Thomas Witkowski's avatar
Thomas Witkowski committed
1175
1176
      request[requestCounter++] = 
	mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
1177
1178
    }

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

1181
1182
    for (int j = 0; j < static_cast<int>(sendBuffers.size()); j++)
      delete [] sendBuffers[j];
1183
1184

    i = 0;
1185
1186
    for (RankToDofContainer::iterator recvIt = recvDofs.begin();
	 recvIt != recvDofs.end(); ++recvIt) {      
1187
      int j = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1188
      for (DofContainer::iterator dofIt = recvIt->second.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1189
	   dofIt != recvIt->second.end(); ++dofIt) {
1190

1191
1192
	rankDofsNewGlobalIndex[*dofIt] = recvBuffers[i][j];
	isRankDof[rankDofsNewLocalIndex[*dofIt]] = false;
1193
1194
1195
1196
	j++;
      }

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

1199

1200
    // === Create now the local to global index and local to dof index mappings.  ===
1201

1202
1203
1204
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
  }
1205

1206
1207
1208
1209
  void ParallelDomainBase::createLocalMappings(DofIndexMap &rankDofsNewLocalIndex,
					       DofIndexMap &rankOwnedDofsNewLocalIndex,
					       DofIndexMap &rankDofsNewGlobalIndex)
  {
1210
    mapLocalGlobalDOFs.clear();
1211
    mapLocalToDofIndex.clear();
1212
1213
1214
1215

    for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin();
	 dofIt != rankDofsNewLocalIndex.end(); ++dofIt) {
      DegreeOfFreedom localDof = dofIt->second;
1216
      DegreeOfFreedom globalDof = rankDofsNewGlobalIndex[dofIt->first];
1217
1218
1219
1220
1221
1222
1223
1224

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1227
1228
  void ParallelDomainBase::addAllVertexDOFs(Element *el, int ithEdge, 
					    DofContainer& dofs)
1229
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1230
    FUNCNAME("ParallelDomainBase::addAllVertexDOFs()");
1231

Thomas Witkowski's avatar
Thomas Witkowski committed
1232
1233
1234
    switch (ithEdge) {
    case 0:
      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1235
	addAllVertexDOFs(el->getSecondChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1236
	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1237
	addAllVertexDOFs(el->getSecondChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1238
1239
1240
1241
      }
      break;
    case 1:
      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1242
	addAllVertexDOFs(el->getFirstChild()->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1243
	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1244
	addAllVertexDOFs(el->getFirstChild()->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1245
1246
1247
1248
      }
      break;
    case 2:      
      if (el->getFirstChild()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1249
	addAllVertexDOFs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1250
	dofs.push_back(el->getFirstChild()->getDOF(2));
Thomas Witkowski's avatar
Thomas Witkowski committed
1251
	addAllVertexDOFs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1252
1253
1254
1255
1256
1257
1258
1259
      }
      break;      
    default:
      ERROR_EXIT("Should never happen!\n");
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
1260
1261
1262
1263
1264
1265
1266
1267
1268
  void ParallelDomainBase::addAllEdgeDOFs(Element *el, int ithEdge,
					  DofContainer& dofs)
  {
    FUNCNAME("ParallelDomainBase::addAllEdgeDOFs()");

    bool addThisEdge = false;

    switch (ithEdge) {
    case 0:
1269
      if (el->getSecondChild())
Thomas Witkowski's avatar
Thomas Witkowski committed
1270
	addAllEdgeDOFs(el->getSecondChild(), 2, dofs);
1271
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
1272
1273
1274
1275
	addThisEdge = true;

      break;
    case 1:
1276
      if (el->getFirstChild())
Thomas Witkowski's avatar
Thomas Witkowski committed
1277
	addAllEdgeDOFs(el->getFirstChild(), 2, dofs);
1278
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
1279
1280
1281
1282
1283
1284
	addThisEdge = true;

      break;
    case 2:
      if (el->getFirstChild()) {
	addAllEdgeDOFs(el->getFirstChild(), 0, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1285
	addAllEdgeDOFs(el->getSecondChild(), 1, dofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
      } else {
	addThisEdge = true;
      } 

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

    if (addThisEdge) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1296
      ElementDofIterator elDofIter(feSpace, true);
Thomas Witkowski's avatar
Thomas Witkowski committed
1297
1298
      elDofIter.reset(el);
      do {
1299
1300
1301
	if (elDofIter.getCurrentPos() == 1 && 
	    elDofIter.getCurrentElementPos() == ithEdge)
	  dofs.push_back(elDofIter.getDofPtr());	
Thomas Witkowski's avatar
Thomas Witkowski committed
1302
1303
1304
1305
1306
      } while(elDofIter.next());      
    }
  }


Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1307
1308
1309
1310
  void ParallelDomainBase::createDOFMemberInfo(DofToPartitions& partitionDofs,
					       DofContainer& rankOwnedDofs,
					       DofContainer& rankAllDofs,
					       DofToRank& boundaryDofs)
1311
  {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1312
1313
1314
1315
1316
    partitionDofs.clear();
    rankOwnedDofs.clear();
    rankAllDofs.clear();
    boundaryDofs.clear();

1317
    // === Determine to each dof the set of partitions the dof belongs to. ===
1318

Thomas Witkowski's avatar
Thomas Witkowski committed
1319
    ElementDofIterator elDofIt(feSpace);
1320

1321
1322
1323
1324
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *element = elInfo->getElement();
1325
1326
1327
      elDofIt.reset(element);
      do {
	// Determine to each dof the partition(s) it corresponds to.
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1328
	partitionDofs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
1329
1330
      } while(elDofIt.next());

1331
1332
1333
      elInfo = stack.traverseNext(elInfo);
    }

1334

1335
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
1336

1337
    // iterate over all DOFs
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1338
1339
    for (DofToPartitions::iterator it = partitionDofs.begin();
	 it != partitionDofs.end(); ++it) {
1340
1341

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

1345
	if (*itpart1 == mpiRank) {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1346
1347
	  rankAllDofs.push_back(it->first);

1348
	  if (it->second.size() == 1) {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1349
	    rankOwnedDofs.push_back(it->first);
1350
1351
1352
1353
1354
1355
1356
	  } 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
1357
		 itpart2 != it->second.end(); ++itpart2) {
1358
1359
1360
1361
1362
1363
1364
1365
	      if (*itpart2 > mpiRank)
		insert = false;

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

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

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1368
	    boundaryDofs[it->first] = highestRank;
1369
1370
1371
1372
	  }

	  break;
	}
1373

1374
1375
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1376

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
1377
1378
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
    sort(rankOwnedDofs.begin(), rankOwnedDofs.end(), cmpDofsByValue);
1379
1380
1381
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
  void ParallelDomainBase::DbgCreateElementMap(ElementIdxToDofs &elMap)
  {
    FUNCNAME("ParallelDomainBase::DbgCreateElementMap()");

    TEST_EXIT(mesh->getDim() == 2)("Works only for 2d!\n");

    elMap.clear();
       
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *el = elInfo->getElement();
      orderDOFs(el->getDOF(0), el->getDOF(1), el->getDOF(2), elMap[el->getIndex()]);
      elInfo = stack.traverseNext(elInfo);
    }
  }

  void ParallelDomainBase::DbgTestElementMap(ElementIdxToDofs &elMap)
  {
    FUNCNAME("ParallelDomainbase::DbgTestElementMap()");

    TEST_EXIT(mesh->getDim() == 2)("Works only for 2d!\n");

    DofContainer vec(3);
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *el = elInfo->getElement();
      orderDOFs(el->getDOF(0), el->getDOF(1), el->getDOF(2), vec);
      for (int i = 0; i < 3; i++) {
	if (elMap[el->getIndex()][i] != vec[i]) {
	  std::cout << "[DBG " << mpiRank 
		    << "]: Wrong new dof numeration in element = " 
		    << el->getIndex() << std::endl;
	  std::cout << "[DBG " << mpiRank << "]: Old numeration was: ";
	  for (int j = 0; j < 3; j++)
	    std::cout << elMap[el->getIndex()][j] << " = " 
		      << *(elMap[el->getIndex()][j]) << "  ";
	  std::cout << std::endl;
	  std::cout << "[DBG " << mpiRank << "]: New numeration is:  ";
	  for (int j = 0; j < 3; j++)
	    std::cout << vec[j] << " = "  << *(vec[j]) << "  ";
	  std::cout << std::endl;
	  ERROR_EXIT("WRONG NEW DOF NUMERATION!\n");
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
  void ParallelDomainBase::DbgTestInteriorBoundary()
  {
    FUNCNAME("ParallelDomainBase::DbgTestInteriorBoundary()");

    std::vector<int*> sendBuffers, recvBuffers;

    MPI::Request request[myIntBoundary.boundary.size() + 
			 otherIntBoundary.boundary.size()];
    int requestCounter = 0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1445
1446
      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
1447
      for (int i = 0; i < nSendInt; i++)
1448
	buffer[i] = (rankIt->second)[i].rankObject.elIndex;
1449
      
Thomas Witkowski's avatar
Thomas Witkowski committed
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
      sendBuffers.push_back(buffer);
      
      request[requestCounter++] =
	mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0);
    }

    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {
      int nRecvInt = rankIt->second.size();
      int *buffer = new int[nRecvInt];
      recvBuffers.push_back(buffer);

      request[requestCounter++] = 
	mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0);
    }


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

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

    int bufCounter = 0;
    for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin();
	 rankIt != otherIntBoundary.boundary.end(); ++rankIt) {

      TEST_EXIT(rankIt->second.size() == otherIntBoundary.boundary[rankIt->first].size())
	("Boundaries does not fit together!\n");      

      for (int i = 0; i < static_cast<int>(rankIt->second.size()); i++) {
	int elIndex1 = recvBuffers[bufCounter][i];
1481
	int elIndex2 = otherIntBoundary.boundary[rankIt->first][i].neighbourObject.elIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
1482
1483
1484
1485
1486
1487
1488
1489
1490

	TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at interior boundary!\n");
      }

      delete [] recvBuffers[bufCounter++];
    }
  }
  

1491
  void ParallelDomainBase::DbgTestCommonDofs(bool printCoords)
Thomas Witkowski's avatar
Thomas Witkowski committed
1492
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1493
    FUNCNAME("ParallelDomainBase::DbgTestCommonDofs()");
1494

Thomas Witkowski's avatar
Thomas Witkowski committed
1495
1496
    // Maps to each neighbour rank an array of WorldVectors. This array contain the 
    // coordinates of all dofs this rank shares on the interior boundary with the 
1497
1498
    // neighbour rank. A rank sends the coordinates to another rank, if it owns the
    // boundarys dof.
Thomas Witkowski's avatar
Thomas Witkowski committed
1499
    RankToCoords sendCoords;
1500
1501
1502
    // A rank receives all boundary dofs that are at its interior boundaries but are
    // not owned by the rank. This map stores for each rank the the coordinates of dofs
    // this rank expectes to receive from.
Thomas Witkowski's avatar
Thomas Witkowski committed
1503
    RankToCoords recvCoords;
1504
1505
1506

    WorldVector<double> coords;
   
Thomas Witkowski's avatar
Thomas Witkowski committed
1507
    for (RankToDofContainer::iterator it = sendDofs.begin();
1508
1509
1510
1511
1512
	 it != sendDofs.end(); ++it) {
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
	TEST_EXIT(b)("Cannot find DOF in mesh!\n");
1513

1514
1515
1516
1517
	sendCoords[it->first].push_back(coords);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1518
    for (RankToDofContainer::iterator it = recvDofs.begin();
1519
1520
1521
1522
1523
1524
1525
1526
1527
	 it != recvDofs.end(); ++it) {
      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords);
	TEST_EXIT(b)("Cannot find DOF in mesh!\n");
	recvCoords[it->first].push_back(coords);
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1528
1529
1530
1531
1532
1533
    std::vector<int> sendSize(mpiSize, 0);
    std::vector<int> recvSize(mpiSize, 0);
    std::vector<int> recvSizeBuffer(mpiSize, 0);
    MPI::Request request[(mpiSize - 1) * 2];
    int requestCounter = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
1534
    for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
1535
1536
      sendSize[it->first] = it->second.size();

Thomas Witkowski's avatar
Thomas Witkowski committed
1537
    for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
      recvSize[it->first] = it->second.size();

    for (int i = 0; i < mpiSize; i++) {
      if (i == mpiRank)
	continue;

      request[requestCounter++] = mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0);
    }   

    for (int i = 0; i < mpiSize; i++) {
      if (i == mpiRank)
	continue;

      request[requestCounter++] = mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0);
    }

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


    for (int i = 0; i < mpiSize; i++) {
      if (i == mpiRank)
	continue;

      if (recvSize[i] != recvSizeBuffer[i]) {
	std::cout << "Error: MPI rank " << mpiRank << " expectes to receive " 
		  << recvSize[i] << " DOFs from rank " << i << ". But this rank sends "
		  << recvSizeBuffer[i] << " DOFs!" << std::endl;

	ERROR_EXIT("Should not happen!\n");
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1570
    std::map<int, double*> sendCoordsBuffer, recvCoordsBuffer;
Thomas Witkowski's avatar
Thomas Witkowski committed
1571
1572
    requestCounter = 0;

1573
    int dimOfWorld = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
1574

1575
1576
    for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) {
      sendCoordsBuffer[it->first] = new double[it->second.size() * dimOfWorld];
Thomas Witkowski's avatar
Thomas Witkowski committed
1577

1578
1579
1580
1581
      for (int i = 0; i < static_cast<int>(it->second.size()); i++) 
	for (int j = 0; j < dimOfWorld; j++)
	  sendCoordsBuffer[it->first][i * dimOfWorld + j] = (it->second)[i][j];
      
Thomas Witkowski's avatar
Thomas Witkowski committed
1582
      request[requestCounter++] = mpiComm.Isend(sendCoordsBuffer[it->first],
1583
						it->second.size() * dimOfWorld,
Thomas Witkowski's avatar
Thomas Witkowski committed
1584
						MPI_DOUBLE, it->first, 0);
1585
1586
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1587
    for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) {
1588
      recvCoordsBuffer[it->first] = new double[it->second.size() * dimOfWorld];
Thomas Witkowski's avatar
Thomas Witkowski committed
1589
1590

      request[requestCounter++] = mpiComm.Irecv(recvCoordsBuffer[it->first], 
1591
						it->second.size() * dimOfWorld,
Thomas Witkowski's avatar
Thomas Witkowski committed
1592
1593
1594
1595
1596
						MPI_DOUBLE, it->first, 0);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1597
    for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
1598
1599
      delete [] sendCoordsBuffer[it->first];

1600
    double eps = 1e-13;
1601

Thomas Witkowski's avatar
Thomas Witkowski committed
1602
    for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1603
      for (int i = 0; i < static_cast<int>(it->second.size()); i++) {
1604
1605
1606
	for (int j = 0; j < dimOfWorld; j++) {
	  bool c = fabs((it->second)[i][j] - 
			recvCoordsBuffer[it->first][i * dimOfWorld + j]) < eps;
1607
	  
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
	  if (printCoords && !c) {
	    std::cout << "[DBG] i = " << i << std::endl;
	    
	    std::cout.precision(5);
	    std::cout << "[DBG]  "
		      << "Rank " << mpiRank << " from rank " << it->first
		      << " expect coords (";
	    for (int k = 0; k < dimOfWorld; k++)
	      std::cout << (it->second)[i][k] << " , ";
	    std::cout << ")  received coords (";
	    for (int k = 0; k < dimOfWorld; k++)
	      std::cout << recvCoordsBuffer[it->first][i * dimOfWorld + k] << " , ";
	    std::cout << ")" << std::endl;
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
1622

1623
1624
	  TEST_EXIT(c)("Wrong DOFs in rank %d!\n", mpiRank);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
1625
1626
1627
      } 

      delete [] recvCoordsBuffer[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
1628
1629
1630
    }
  }

1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
  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;
  }


1651
}
For faster browsing, not all history is shown. View entire blame