ParallelDomainBase.cc 90.1 KB
Newer Older
2001

2002
    // === Create now the local to global index and local to dof index mappings.  ===
2003

2004
2005
    createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
			rankDofsNewGlobalIndex);
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019

    nRankRows = nRankDofs * nComponents;
    nOverallRows = nOverallDofs * nComponents;

    // === Update dof admins due to new number of dofs. ===
  
    updateDofAdmins();

    lastMeshChangeIndex = mesh->getChangeIndex();

#if (DEBUG != 0)
    dbgTestElementMap(elMap);
    dbgTestCommonDofs(true);
#endif
2020
  }
2021

2022
2023
2024
2025
  void ParallelDomainBase::createLocalMappings(DofIndexMap &rankDofsNewLocalIndex,
					       DofIndexMap &rankOwnedDofsNewLocalIndex,
					       DofIndexMap &rankDofsNewGlobalIndex)
  {
2026
    mapLocalGlobalDofs.clear();
2027
    mapLocalToDofIndex.clear();
2028
2029
2030
2031

    for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin();
	 dofIt != rankDofsNewLocalIndex.end(); ++dofIt) {
      DegreeOfFreedom localDof = dofIt->second;
2032
      DegreeOfFreedom globalDof = rankDofsNewGlobalIndex[dofIt->first];
2033
2034

      *const_cast<DegreeOfFreedom*>(dofIt->first) = localDof;
2035
      mapLocalGlobalDofs[localDof] = globalDof;
2036
2037
2038
2039
2040
    }

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

2043
  
2044
  void ParallelDomainBase::createDofMemberInfo(DofToPartitions& partitionDofs,
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
2045
2046
					       DofContainer& rankOwnedDofs,
					       DofContainer& rankAllDofs,
2047
2048
					       DofToRank& boundaryDofs,
					       DofToBool& vertexDof)
2049
  {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
2050
2051
2052
2053
    partitionDofs.clear();
    rankOwnedDofs.clear();
    rankAllDofs.clear();
    boundaryDofs.clear();
2054
    vertexDof.clear();
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
2055

2056
    // === Determine to each dof the set of partitions the dof belongs to. ===
2057

Thomas Witkowski's avatar
Thomas Witkowski committed
2058
    ElementDofIterator elDofIt(feSpace);
2059

2060
2061
2062
2063
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *element = elInfo->getElement();
2064
2065
2066
      elDofIt.reset(element);
      do {
	// Determine to each dof the partition(s) it corresponds to.
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
2067
	partitionDofs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]);
2068
2069
2070
2071
2072

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

2075
2076
2077
      elInfo = stack.traverseNext(elInfo);
    }

2078

2079
    // === Determine the set of ranks dofs and the dofs ownership at the boundary. ===
2080

2081
    // iterate over all DOFs
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
2082
2083
    for (DofToPartitions::iterator it = partitionDofs.begin();
	 it != partitionDofs.end(); ++it) {
2084

2085
2086
      bool isInRank = false;

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

2091
	if (*itpart1 == mpiRank) {	  
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
2092
2093
	  rankAllDofs.push_back(it->first);

2094
	  if (it->second.size() == 1) {
Thomas Witkowski's avatar
n    
Thomas Witkowski committed
2095
	    rankOwnedDofs.push_back(it->first);
2096
2097
2098
2099
2100
2101
2102
	  } 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
2103
		 itpart2 != it->second.end(); ++itpart2) {
2104
2105
2106
2107
2108
2109
2110
2111
	      if (*itpart2 > mpiRank)
		insert = false;

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

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

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
2114
	    boundaryDofs[it->first] = highestRank;
2115
2116
	  }

2117
	  isInRank = true;
2118
2119
	  break;
	}
2120

2121
      }
2122

2123
      if (!isInRank)
2124
	vertexDof.erase(it->first);
2125
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
2126

Thomas Witkowski's avatar
n    
Thomas Witkowski committed
2127
2128
    sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue);
    sort(rankOwnedDofs.begin(), rankOwnedDofs.end(), cmpDofsByValue);
2129
2130
  }

2131
2132
2133
2134
  void ParallelDomainBase::createPeriodicMap()
  {
    FUNCNAME("ParallelDomainBase::createPeriodicMap()");

2135
2136
2137
    if (periodicBoundary.boundary.size() == 0)
      return;

2138
    // Clear all periodic dof mappings calculated before. We do it from scratch.
2139
2140
    periodicDof.clear();

2141
    StdMpi<std::vector<int> > stdMpi(mpiComm, false);
2142

2143
2144
2145
    // === Each rank traverse its periodic boundaries and sends the dof indices ===
    // === to the rank "on the other side" of the periodic boundary.            ===

2146
2147
2148
2149
2150
    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");

2151
      // Create dof indices on the boundary. 
2152
      DofContainer dofs;
2153
2154
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
2155
2156
	boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, dofs);
	boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, dofs);
2157
      }
2158

2159
      // Send the global indices to the rank on the other side.
2160
      stdMpi.getSendData(it->first).reserve(dofs.size());
2161
      for (int i = 0; i < static_cast<int>(dofs.size()); i++)
2162
	stdMpi.getSendData(it->first).push_back(mapLocalGlobalDofs[*(dofs[i])]);
2163

2164
      // Receive from this rank the same number of dofs.
2165
      stdMpi.recv(it->first, dofs.size());     
2166
2167
    }

2168
2169
    stdMpi.updateSendDataSize();
    stdMpi.startCommunication<int>(MPI_INT);
2170

2171
2172
2173
2174
2175
    // === 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.                                        ===

2176
2177

    std::map<DegreeOfFreedom, std::set<int> > dofFromRank;
2178
2179
2180
2181
    for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin();
	 it != periodicBoundary.boundary.end(); ++it) {
      DofContainer dofs;
      
2182
      // Create the dofs on the boundary in inverse order.
2183
2184
2185
      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
	   boundIt != it->second.end(); ++boundIt) {
	DofContainer tmpdofs;
2186
2187
2188
	boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, tmpdofs);
	boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, tmpdofs);

2189
2190
	for (int j = static_cast<int>(tmpdofs.size()) - 1; j >= 0; j--)
	  dofs.push_back(tmpdofs[j]);	
2191
2192
      }

2193
      // Added the received dofs to the mapping.
2194
      for (int j = 0; j < static_cast<int>(dofs.size()); j++) {
2195
	int globalDofIndex = mapLocalGlobalDofs[*(dofs[j])];
2196
	periodicDof[globalDofIndex].insert(stdMpi.getRecvData(it->first)[j]);
2197
2198
	dofFromRank[globalDofIndex].insert(it->first);
      }
2199
    }
2200

2201
2202
2203
    if (dofFromRank.size() > 0) 
      TEST_EXIT_DBG(mesh->getDim() == 2)
	("Periodic boundary corner problem must be generalized to 3d!\n");
2204

2205
2206
2207
2208
    MPI::Request request[min(static_cast<int>(periodicBoundary.boundary.size() * 2), 4)];
    int requestCounter = 0;
    std::vector<int*> sendBuffers, recvBuffers;

2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
    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);

2240
    int i = 0;
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
    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++;
      }
    }
2252
2253
2254
  }


2255
2256
2257
2258
  Flag ParallelDomainBase::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
  {
    FUNCNAME("ParallelDomainBase::oneIteration()");

2259
    Flag flag = buildAndAdapt(adaptInfo, toDo);
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273

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


2274
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
  Flag ParallelDomainBase::buildAndAdapt(AdaptInfo *adaptInfo, Flag toDo)
  {
        FUNCNAME("StandardProblemIteration::buildAndAdapt()");

    Flag flag = 0, markFlag = 0;
    ProblemStatBase *problem = iterationIF->getProblem();

    if (toDo.isSet(MARK))
      markFlag = problem->markElements(adaptInfo);
    else
      markFlag = 3;

    if (toDo.isSet(BUILD))
      problem->buildBeforeRefine(adaptInfo, markFlag);

    // refine
    if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
      flag = problem->refineMesh(adaptInfo);

    if (toDo.isSet(BUILD))
      problem->buildBeforeCoarsen(adaptInfo, markFlag);

    // coarsen
    if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))
      flag |= problem->coarsenMesh(adaptInfo);

    checkMeshChange();

    if (toDo.isSet(BUILD))
      problem->buildAfterCoarsen(adaptInfo, markFlag, true, true);

    if (toDo.isSet(BUILD_RHS))
      problem->buildAfterCoarsen(adaptInfo, markFlag, false, true);

    return flag;
  }


2312
2313
2314
  void ParallelDomainBase::serialize(std::ostream &out)
  {
    SerUtil::serialize(out, elemWeights);
2315
    SerUtil::serialize(out, initialPartitionMesh);
2316
2317
    SerUtil::serialize(out, partitionVec);
    SerUtil::serialize(out, oldPartitionVec);
2318
    SerUtil::serialize(out, nRankDofs);
2319
2320
2321
2322
    
    myIntBoundary.serialize(out);
    otherIntBoundary.serialize(out);

2323
2324
    serialize(out, sendDofs);
    serialize(out, recvDofs);
2325

2326
    SerUtil::serialize(out, mapLocalGlobalDofs);
2327
2328
2329
    SerUtil::serialize(out, mapLocalToDofIndex);
    SerUtil::serialize(out, isRankDof);

2330
    serialize(out, vertexDof);
2331
    serialize(out, periodicDof);
2332

2333
2334
2335
    SerUtil::serialize(out, rstart);
    SerUtil::serialize(out, nRankRows);
    SerUtil::serialize(out, nOverallRows);
2336
2337
2338
2339
2340
  }


  void ParallelDomainBase::deserialize(std::istream &in)
  {
2341
    SerUtil::deserialize(in, elemWeights);
2342
    SerUtil::deserialize(in, initialPartitionMesh);
2343
2344
    SerUtil::deserialize(in, partitionVec);
    SerUtil::deserialize(in, oldPartitionVec);
2345
    SerUtil::deserialize(in, nRankDofs);    
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
2362
2363
2364
2365
2366
2367
2368
2369
2370

    // Create two maps: one from from element indices to the corresponding element 
    // pointers, and one map from Dof indices to the corresponding dof pointers.
    std::map<int, Element*> elIndexMap;
    std::map<int, const DegreeOfFreedom*> dofMap;
    ElementDofIterator elDofIter(feSpace);
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *el = elInfo->getElement();
      elIndexMap[el->getIndex()] = el;
      elInfo = stack.traverseNext(elInfo);

      elDofIter.reset(el);
      do {
	dofMap[elDofIter.getDof()] = elDofIter.getDofPtr();
      } while(elDofIter.next());      
    }

    myIntBoundary.deserialize(in, elIndexMap);
    otherIntBoundary.deserialize(in, elIndexMap);

    deserialize(in, sendDofs, dofMap);
    deserialize(in, recvDofs, dofMap);

2371
    SerUtil::deserialize(in, mapLocalGlobalDofs);
2372
2373
2374
    SerUtil::deserialize(in, mapLocalToDofIndex);
    SerUtil::deserialize(in, isRankDof);

2375
    deserialize(in, vertexDof, dofMap);
2376

2377
2378
2379
2380
2381
    SerUtil::deserialize(in, rstart);
    SerUtil::deserialize(in, nRankRows);
    SerUtil::deserialize(in, nOverallRows);

    deserialized = true;
2382
2383
  }

2384
2385
2386
2387
  void ParallelDomainBase::serialize(std::ostream &out, PeriodicDofMap &data)
  {
    int mapSize = data.size();
    SerUtil::serialize(out, mapSize);
2388

2389
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410
    for (PeriodicDofMap::iterator it = data.begin(); it != data.end(); ++it) {
      DegreeOfFreedom dof = it->first;
      std::set<DegreeOfFreedom> dofSet = it->second;
      SerUtil::serialize(out, dof);
      SerUtil::serialize(out, dofSet);
    }
  }

  void ParallelDomainBase::deserialize(std::istream &in, PeriodicDofMap &data)
  {
    data.clear();

    int mapSize = 0;
    SerUtil::deserialize(in, mapSize);
    for (int i = 0; i < mapSize; i++) {
      DegreeOfFreedom dof = 0;
      std::set<DegreeOfFreedom> dofSet;

      SerUtil::deserialize(in, dof);
      SerUtil::deserialize(in, dofSet);

      data[dof] = dofSet;
2411
    }    
2412
2413
  }
  
2414
  void ParallelDomainBase::dbgCreateElementMap(ElementIdxToDofs &elMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
2415
  {
2416
    FUNCNAME("ParallelDomainBase::dbgCreateElementMap()");
Thomas Witkowski's avatar
Thomas Witkowski committed
2417

2418
    int dim = mesh->getDim();
Thomas Witkowski's avatar
Thomas Witkowski committed
2419
2420
2421
2422
2423
2424
    elMap.clear();
       
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *el = elInfo->getElement();
2425
2426
2427
2428
2429
2430
2431
2432
2433
2434
      switch (dim) {
      case 2:
	orderDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), elMap[el->getIndex()]);
	break;
      case 3:
	orderDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), el->getDOF(3), elMap[el->getIndex()]);
	break;
      default:
	ERROR_EXIT("What is this?\n");
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
2435
2436
2437
2438
      elInfo = stack.traverseNext(elInfo);
    }
  }

2439
  void ParallelDomainBase::dbgTestElementMap(ElementIdxToDofs &elMap)
Thomas Witkowski's avatar
Thomas Witkowski committed
2440
  {
2441
    FUNCNAME("ParallelDomainbase::dbgTestElementMap()");
Thomas Witkowski's avatar
Thomas Witkowski committed
2442

2443
2444
2445
    int dim = mesh->getDim();
    int nVertex = Global::getGeo(VERTEX, dim);
    DofContainer vec(nVertex);
Thomas Witkowski's avatar
Thomas Witkowski committed
2446
2447
2448
2449
2450

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      Element *el = elInfo->getElement();
2451
2452
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463

      switch (dim) {
      case 2:
	orderDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), vec);
	break;
      case 3:
	orderDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), el->getDOF(3), vec);
	break;
      default:
	ERROR_EXIT("What is this?\n");
      }

      for (int i = 0; i < nVertex; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2464
2465
2466
2467
2468
	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: ";
2469
	  for (int j = 0; j < nVertex; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
2470
2471
2472
2473
	    std::cout << elMap[el->getIndex()][j] << " = " 
		      << *(elMap[el->getIndex()][j]) << "  ";
	  std::cout << std::endl;
	  std::cout << "[DBG " << mpiRank << "]: New numeration is:  ";
2474
	  for (int j = 0; j < nVertex; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
2475
2476
2477
2478
2479
2480
2481
2482
2483
	    std::cout << vec[j] << " = "  << *(vec[j]) << "  ";
	  std::cout << std::endl;
	  ERROR_EXIT("WRONG NEW DOF NUMERATION!\n");
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }
  }

2484
  void ParallelDomainBase::dbgTestInteriorBoundary()
Thomas Witkowski's avatar
Thomas Witkowski committed
2485
  {
2486
    FUNCNAME("ParallelDomainBase::dbgTestInteriorBoundary()");
Thomas Witkowski's avatar
Thomas Witkowski committed
2487
2488
2489
2490
2491
2492
2493
2494
2495

    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) {
2496

Thomas Witkowski's avatar
Thomas Witkowski committed
2497
2498
      int nSendInt = rankIt->second.size();
      int* buffer = new int[nSendInt];
2499
      for (int i = 0; i < nSendInt; i++)
2500
	buffer[i] = (rankIt->second)[i].rankObj.elIndex;
2501
      
Thomas Witkowski's avatar
Thomas Witkowski committed
2502
2503
2504
2505
2506
2507
2508
2509
2510
2511
2512
2513
2514
2515
2516
2517
2518
2519
2520
2521
2522
2523
2524
2525
2526
2527
2528
2529
2530
2531
      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];
2532
	int elIndex2 = otherIntBoundary.boundary[rankIt->first][i].neighObj.elIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
2533
2534
2535
2536
2537
2538
2539
2540
2541

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

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

2542
  void ParallelDomainBase::dbgTestCommonDofs(bool printCoords)
Thomas Witkowski's avatar
Thomas Witkowski committed
2543
  {
2544
    FUNCNAME("ParallelDomainBase::dbgTestCommonDofs()");
2545

2546
2547
2548
2549
2550
2551
2552
2553
2554
    clock_t first = clock();

    int testCommonDofs = 1;
    GET_PARAMETER(0, "dbg->test common dofs", "%d", &testCommonDofs);
    if (testCommonDofs == 0) {
      MSG("Skip test common dofs!\n");
      return;
    }

2555
    // Maps to each neighbour rank an array of WorldVectors. This array contains the 
Thomas Witkowski's avatar
Thomas Witkowski committed
2556
    // coordinates of all dofs this rank shares on the interior boundary with the 
2557
2558
    // neighbour rank. A rank sends the coordinates to another rank, if it owns the
    // boundarys dof.
Thomas Witkowski's avatar
Thomas Witkowski committed
2559
    RankToCoords sendCoords;
2560

2561
    // A rank receives all boundary dofs that are at its interior boundaries but are
2562
    // not owned by the rank. This map stores for each rank the coordinates of dofs
2563
    // this rank expectes to receive from.
Thomas Witkowski's avatar
Thomas Witkowski committed
2564
    RankToCoords recvCoords;
2565

2566
2567
2568
    DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds");
    mesh->getDofIndexCoords(feSpace, coords);
  
Thomas Witkowski's avatar
Thomas Witkowski committed
2569
    for (RankToDofContainer::iterator it = sendDofs.begin();
2570
2571
2572
2573
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin(); 
	   dofIt != it->second.end(); ++dofIt)
	sendCoords[it->first].push_back(coords[**dofIt]);
2574

Thomas Witkowski's avatar
Thomas Witkowski committed
2575
    for (RankToDofContainer::iterator it = recvDofs.begin();
2576
2577
2578
2579
	 it != recvDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	recvCoords[it->first].push_back(coords[**dofIt]);
2580

Thomas Witkowski's avatar
Thomas Witkowski committed
2581
2582
2583
2584
2585
2586
    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
2587
    for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
2588
2589
      sendSize[it->first] = it->second.size();

Thomas Witkowski's avatar
Thomas Witkowski committed
2590
    for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
2591
2592
2593
2594
2595
2596
2597
2598
2599
2600
2601
2602
2603
2604
2605
2606
2607
2608
2609
2610
2611
2612
2613
2614
2615
2616
2617
2618
2619
2620
2621
2622
      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");
      }
    }

2623
2624
2625
    // === Now we know that the number of send and received DOFs fits together. ===
    // === So we can check if also the coordinates of the communicated DOFs are ===
    // === the same on both corresponding ranks.                                ===
Thomas Witkowski's avatar
Thomas Witkowski committed
2626

2627
    typedef std::vector<WorldVector<double> > CoordsVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
2628
    StdMpi<CoordsVec> stdMpi(mpiComm, true);
2629
2630
2631
    stdMpi.send(sendCoords);
    stdMpi.recv(recvCoords);   
    stdMpi.startCommunication<double>(MPI_DOUBLE);
Thomas Witkowski's avatar
Thomas Witkowski committed
2632

2633
    double eps = 1e-13;
2634
    int dimOfWorld = Global::getGeo(WORLD);
2635

2636
2637
2638
2639
    // === Compare the received with the expected coordinates. ===

    for (RankToCoords::iterator it = stdMpi.getRecvData().begin(); 
	 it != stdMpi.getRecvData().end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2640
      for (int i = 0; i < static_cast<int>(it->second.size()); i++) {
2641
	for (int j = 0; j < dimOfWorld; j++) {
2642
	  bool c = fabs((it->second)[i][j] - recvCoords[it->first][i][j]) < eps;
2643
	  
2644
2645
	  // === Print error message if the coordinates are not the same. ===

2646
2647
	  if (printCoords && !c) {
	    std::cout.precision(5);
2648
2649
	    std::cout << "[DBG] i = " << i << std::endl;    
	    std::cout << "[DBG] Rank " << mpiRank << " from rank " << it->first
2650
		      << " expect coords (";
2651
	    for (int k = 0; k < dimOfWorld; k++) {
2652
	      std::cout << recvCoords[it->first][i][k];
2653
2654
2655
	      if (k + 1 < dimOfWorld)
		std::cout << " / ";
	    }
2656
	    std::cout << ")  received coords (";
2657
	    for (int k = 0; k < dimOfWorld; k++) {
2658
	      std::cout << (it->second)[i][k];
2659
2660
2661
	      if (k + 1 < dimOfWorld)
		std::cout << " / ";
	    }
2662
	    std::cout << ")" << std::endl;
2663
2664

	    debug::printInfoByDof(feSpace, *(recvDofs[it->first][i]));
2665
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
2666

2667
2668
	  TEST_EXIT(c)("Wrong DOFs in rank %d!\n", mpiRank);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
2669
      } 
Thomas Witkowski's avatar
Thomas Witkowski committed
2670
    }
2671
2672

    INFO(info, 8)("Test common dofs needed %.5f seconds\n", TIME_USED(first, clock()));
Thomas Witkowski's avatar
Thomas Witkowski committed
2673
2674
  }

2675

2676
  void ParallelDomainBase::printMapLocalGlobal(int rank)
2677
2678
2679
2680
  {    
    if (rank == -1 || mpiRank == rank) {
      std::cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << std::endl;
      
2681
2682
      for (DofMapping::iterator it = mapLocalGlobalDofs.begin();
	   it != mapLocalGlobalDofs.end(); it++) {
2683
2684
2685
2686
2687
2688
2689
2690
2691
2692
2693
2694
2695
2696
2697
2698
2699
2700
2701
2702
2703
2704
2705
2706
2707
2708
2709
2710
2711
2712
	DegreeOfFreedom localdof = -1;
	if (mapLocalToDofIndex.count(it->first) > 0)
	  localdof = mapLocalToDofIndex[it->first];
	
	std::cout << "DOF " << it->first << " " 
		  << it->second << " " 
		  << localdof << std::endl;
	WorldVector<double> coords;
	mesh->getDofIndexCoords(it->first, feSpace, coords);
	coords.print();
	for (RankToDofContainer::iterator rankit = sendDofs.begin();
	     rankit != sendDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second.begin();
	       dofit != rankit->second.end(); ++dofit)
	    if (**dofit == it->first)
	      std::cout << "SEND DOF TO " << rankit->first << std::endl;	  
	}
	for (RankToDofContainer::iterator rankit = recvDofs.begin();
	     rankit != recvDofs.end(); ++rankit) {
	  for (DofContainer::iterator dofit = rankit->second.begin();
	       dofit != rankit->second.end(); ++dofit)
	    if (**dofit == it->first)
	      std::cout << "RECV DOF FROM " << rankit->first << std::endl;	  
	}
	std::cout << "------" << std::endl;
      }
    }
  }


2713
  void ParallelDomainBase::printMapPeriodic(int rank)
2714
  {
2715
2716
    FUNCNAME("ParallelDomainBase::printMapPeriodic()");

2717
2718
2719
2720
2721
2722
2723
2724
2725
2726
2727
    if (rank == -1 || mpiRank == rank) {
      std::cout << "====== DOF MAP PERIODIC ====== " << std::endl;

      for (PeriodicDofMap::iterator it = periodicDof.begin();
	   it != periodicDof.end(); ++it) {
	std::cout << "DOF MAP " << it->first << ": ";
	for (std::set<DegreeOfFreedom>::iterator dofit = it->second.begin();
	     dofit != it->second.end(); ++dofit)
	  std::cout << *dofit << " ";
	std::cout << std::endl;

2728
	DegreeOfFreedom localdof = -1;
2729
2730
	for (DofMapping::iterator dofIt = mapLocalGlobalDofs.begin();
	     dofIt != mapLocalGlobalDofs.end(); ++dofIt)
2731
2732
2733
	  if (dofIt->second == it->first)
	    localdof = dofIt->first;

2734
2735
	TEST_EXIT(localdof != -1)("There is something wrong!\n");

2736
2737
2738
2739
2740
	WorldVector<double> coords;
	mesh->getDofIndexCoords(localdof, feSpace, coords);
	coords.print();
      }
    }
2741
2742
2743
2744
2745
2746
2747
2748
2749
2750
2751
2752
2753
2754
2755
2756
2757
  }

  
  void ParallelDomainBase::printRankDofs(int rank, DofContainer& rankDofs,
					 DofContainer& rankAllDofs)
  {
    if (rank == -1 || mpiRank == rank) {
      std::cout << "====== RANK DOF INFORMATION ====== " << std::endl;

      std::cout << "  RANK OWNED DOFS: " << std::endl;
      for (DofContainer::iterator dofit = rankDofs.begin();
	   dofit != rankDofs.end(); ++dofit) {
	std::cout << "    " << **dofit << std::endl;
	WorldVector<double> coords;
	mesh->getDofIndexCoords(*dofit, feSpace, coords);
	coords.print();
      }
2758

2759
2760
2761
2762
2763
2764
2765
2766
2767
      std::cout << "  RANK ALL DOFS: " << std::endl;
      for (DofContainer::iterator dofit = rankAllDofs.begin();
	   dofit != rankAllDofs.end(); ++dofit) {
	std::cout << "    " << **dofit << std::endl;
	WorldVector<double> coords;
	mesh->getDofIndexCoords(*dofit, feSpace, coords);
	coords.print();
      }      
    }
2768
2769
2770
  }


2771
2772
2773
2774
2775
2776
2777
2778
2779
2780
2781
2782
2783
2784
2785
2786
2787
2788
  void ParallelDomainBase::writePartitioningMesh(std::string filename)
  {
    FUNCNAME("ParallelDomainBase::writePartitioningMesh()");

    std::map<int, double> vec;    
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, 
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    
    while (elInfo) {		  
      int index = elInfo->getElement()->getIndex();
      vec[index] = partitionVec[index];
      elInfo = stack.traverseNext(elInfo);
    }

    ElementFileWriter::writeFile(vec, feSpace, filename);
  }

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