ParallelDofMapping.cc 12.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include "parallel/ParallelDofMapping.h"
14
#include "parallel/MeshLevelData.h"
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
15
#include "parallel/StdMpi.h"
16 17 18 19 20

namespace AMDiS {

  using namespace std;

21 22 23 24 25

  void DofToMatIndex::getReverse(int rowIndex, int &component, int &dofIndex)
  {
    FUNCNAME("DofToMatIndex::getReverse()");

26
    for (map<int, boost::container::flat_map<DegreeOfFreedom, int> >::iterator it0 = data.begin();
27
	 it0 != data.end(); ++it0)
28
      for (boost::container::flat_map<DegreeOfFreedom, int>::iterator it1 = it0->second.begin();
29 30 31 32 33 34 35 36 37 38 39 40
	   it1 != it0->second.end(); ++it1)
	if (it1->second == rowIndex) {
	  component = it0->first;
	  dofIndex = it1->first;
	  return;
	}
  
    component = -1;
    dofIndex = -1;
  }


41
  ComponentDofMap::ComponentDofMap(MeshLevelData* ld)
42 43 44 45 46 47 48 49 50 51
    : levelData(ld),
      dofComm(NULL),
      feSpace(NULL),
      needGlobalMapping(false),
      isNonLocal(false)
  {
    clear();
  }


52
  void ComponentDofMap::clear()
53 54
  {
    dofMap.clear();
55

56
    nonRankDofs.clear();
57

58 59 60 61
    nRankDofs = 0;
    nLocalDofs = 0;
    nOverallDofs = 0;
    rStartDofs = 0;
62 63 64
  }


65
  void ComponentDofMap::update()
66
  {
67
    FUNCNAME("ComponentDofMap::update()");
68

69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
    // === Compute local indices for all rank owned DOFs. ===
    
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
      if (it->second.local == -1 && nonRankDofs.count(it->first) == 0)
	it->second.local = nRankDofs++;
    
    // === Compute number of local and global DOFs in the mapping. ===
    
    nOverallDofs = 0;
    rStartDofs = 0;
    mpi::getDofNumbering(mpiComm, nRankDofs, rStartDofs, nOverallDofs);
    
    // === If required, compute also the global indices. ===
    
    if (needGlobalMapping) {
      computeGlobalMapping();
85
      
86
      if (isNonLocal)
87
	computeNonLocalIndices();
88
    }
89 90 91
  }


92
  void ComponentDofMap::computeGlobalMapping()
93
  {
94
    FUNCNAME("ComponentDofMap::computeGlobalMapping()");
95

96 97
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
      it->second.global = it->second.local + rStartDofs;
98 99
  }

100

101
  void ComponentDofMap::computeNonLocalIndices()
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
102
  {
103
    FUNCNAME("ComponentDofMap::computeNonLocalIndices()");
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
104

105 106
    TEST_EXIT_DBG(dofComm)("No DOF communicator defined!\n");

Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
107 108
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;

109 110 111
    // === Send all global indices of DOFs that are owned by the rank to all ===
    // === other ranks that also include this DOF.                           ===

112
    StdMpi<vector<int> > stdMpi(mpiComm);
113

114
    for (DofComm::Iterator it(dofComm->getSendDofs(), 0, feSpace); 
115 116 117 118
	 !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (meshLevel > 0)
	rank = levelData->mapRank(rank, 0, meshLevel);
119
      
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
120
      for (; !it.endDofIter(); it.nextDof())
121 122
	if (dofMap.count(it.getDofIndex()) && 
	    !nonRankDofs.count(it.getDofIndex()))
123
	  stdMpi.getSendData(rank).
124
	    push_back(dofMap[it.getDofIndex()].global);
125
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
126 127 128

    stdMpi.updateSendDataSize();

129 130 131

    // === Check from which ranks this rank must receive some data. ===

132
    for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpace); 
133
	 !it.end(); it.nextRank()) {
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
134 135
      bool recvFromRank = false;
      for (; !it.endDofIter(); it.nextDof()) {
136
	if (nonRankDofs.count(it.getDofIndex())) {
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
137 138 139 140 141
	  recvFromRank = true;
	  break;
	}
      }

142 143 144 145 146 147 148
      if (recvFromRank) {
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);

	stdMpi.recv(rank);
      }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
149 150
    }

151 152 153

    // === Start communication to exchange global indices. ===

Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
154 155
    stdMpi.startCommunication();

156 157 158

    // === And set the global indices for all DOFs that are not owned by rank. ===
    
159
    for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpace);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
160
	 !it.end(); it.nextRank()) {
161 162 163 164
      int rank = it.getRank();
      if (meshLevel > 0)
	rank = levelData->mapRank(rank, 0, meshLevel);

Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
165 166
      int i = 0;
      for (; !it.endDofIter(); it.nextDof())
167 168
	if (nonRankDofs.count(it.getDofIndex()))
	  dofMap[it.getDofIndex()].global = stdMpi.getRecvData(rank)[i++];
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
169 170 171 172
    }
  }


173
  void ParallelDofMapping::init(MeshLevelData &ldata,
174
				vector<const FiniteElemSpace*> &fe,
175
				vector<const FiniteElemSpace*> &uniqueFe,
176
				bool b)
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
177
  {
178 179
    FUNCNAME("ParallelDofMapping::init()");

180
    levelData = &ldata;
181
    isNonLocal = b;
182

183
    data->init(fe, uniqueFe, isNonLocal);   
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
184
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
185 186


187 188 189 190
  void ParallelDofMapping::clear()
  {
    FUNCNAME("ParallelDofMapping::clear()");

191 192
    TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n");

193
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
Thomas Witkowski's avatar
Thomas Witkowski committed
194
      it->clear();
195

196 197 198 199 200
    nRankDofs = -1;
    nLocalDofs = -1;
    nOverallDofs = -1;
    rStartDofs = -1;
    dofToMatIndex.clear();
201 202
  }

203 204 205 206 207 208
  
  void ParallelDofMapping::setMpiComm(MPI::Intracomm &m, int l)
  {
    mpiComm = m;
    meshLevel = l;

209 210
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->setMpiComm(m, l);
211 212
  }

213

214
  void ParallelDofMapping::setDofComm(DofComm &dc)
215 216 217
  {
    FUNCNAME("ParallelDofMapping::setDofComm()");

218
    dofComm = &dc;
219 220

    // Add the DOF communicator also to all FE space DOF mappings.
221 222
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->setDofComm(dc);
223
  }
224

225

226
  int ParallelDofMapping::computeRankDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
227
  {
228
    FUNCNAME("ParallelDofMapping::computeRankDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
229 230
    
    int result = 0;
231 232 233
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nRankDofs;
   
Thomas Witkowski's avatar
Thomas Witkowski committed
234 235 236 237
    return result;
  }


238
  int ParallelDofMapping::computeLocalDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
239
  {
240
    FUNCNAME("ParallelDofMapping::computeLocalDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
241 242
    
    int result = 0;
243 244 245
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nLocalDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
246 247 248 249
    return result;
  }


250
  int ParallelDofMapping::computeOverallDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
251
  {
252
    FUNCNAME("ParallelDofMapping::computeOverallDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
253 254

    int result = 0;
255 256 257
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nOverallDofs;
       
Thomas Witkowski's avatar
Thomas Witkowski committed
258 259 260 261
    return result;
  }


262
  int ParallelDofMapping::computeStartDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
263
  {
264
    FUNCNAME("ParallelDofMapping::computeStartDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
265 266

    int result = 0;
267 268 269
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->rStartDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
270 271 272 273 274 275 276 277
    return result;
  }


  void ParallelDofMapping::update()
  {
    FUNCNAME("ParallelDofMapping::update()");

278
    // First, update all FE space DOF mappings.
279 280
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->update();
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    
282 283 284 285 286 287 288
    // Compute all numbers from this mappings.
    nRankDofs = computeRankDofs();
    nLocalDofs = computeLocalDofs();
    nOverallDofs = computeOverallDofs();
    rStartDofs = computeStartDofs();
    
    // And finally, compute the matrix indices.
289
    computeMatIndex(needMatIndexFromGlobal);
Thomas Witkowski's avatar
Thomas Witkowski committed
290 291
  }
  
292 293 294 295 296

  void ParallelDofMapping::update(vector<const FiniteElemSpace*>& fe)
  {
    FUNCNAME("ParallelDofMapping::update()");

297 298 299 300 301 302 303 304 305 306
    ERROR_EXIT("DAS MUSS ICH MIR MAL UEBERLEGEN!\n");

    // for (vector<const FiniteElemSpace*>::iterator it = fe.begin();
    // 	 it != fe.end(); ++it)
    //   if (find(feSpacesUnique.begin(), feSpacesUnique.end(), *it) == 
    // 	  feSpacesUnique.end())
    // 	ERROR_EXIT("Wrong FE space!\n");

    // feSpaces = fe;

307 308 309 310

    update();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
311
  
312
  void ParallelDofMapping::computeMatIndex(bool globalIndex)
Thomas Witkowski's avatar
Thomas Witkowski committed
313 314
  {
    FUNCNAME("ParallelDofMapping::computeMatIndex()");
315

316
    dofToMatIndex.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
317
    
318 319 320
    // The offset is always added to the local matrix index. The offset for the
    // DOFs in the first FE spaces is the smalled global index of a DOF that is
    // owned by the rank.
321
    int offset = rStartDofs;
322 323 324

    // === Create the matrix indices for all component FE spaces. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
325 326
    vector<const FiniteElemSpace*> &feSpaces = data->getFeSpaces();

Thomas Witkowski's avatar
Thomas Witkowski committed
327
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
328 329 330

      // Traverse all DOFs of the FE space and create for all rank owned DOFs
      // a matrix index.
331
      DofMap& dofMap = (*data)[i].getMap();
332
      for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
333
	if ((*data)[i].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
334
	  int globalMatIndex = it->second.local + offset;
335
	  if (globalIndex)
336
	    dofToMatIndex.add(i, it->second.global, globalMatIndex);
337
	  else
338
	    dofToMatIndex.add(i, it->first, globalMatIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
339 340
	}
      }
341

342 343
      // Increase the offset for the next FE space by the number of DOFs owned 
      // by the rank in the current FE space.
344
      offset += (*data)[i].nRankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
345
	
346
      // If there are no non local DOFs, continue with the next FE space.
347
      if (!isNonLocal)
Thomas Witkowski's avatar
Thomas Witkowski committed
348 349
	continue;
      
350
      TEST_EXIT_DBG(dofComm != NULL)("No communicator given!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
351
      
352 353 354
      // === Communicate the matrix indices for all DOFs that are on some ===
      // === interior boundaries.                                         ===

355
      StdMpi<vector<DegreeOfFreedom> > stdMpi(mpiComm);
356
      for (DofComm::Iterator it(dofComm->getSendDofs(), 0, feSpaces[i]); 
Thomas Witkowski's avatar
Thomas Witkowski committed
357 358 359 360 361
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex()))
362
	    if (globalIndex)
363
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
364
	    else
365
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex()));
Thomas Witkowski's avatar
Thomas Witkowski committed
366
	
367 368 369 370 371
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	stdMpi.send(rank, sendGlobalDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
372 373
      }
      
374
      for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]); 
375 376 377 378 379 380 381
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);

	stdMpi.recv(rank);
      }
382

Thomas Witkowski's avatar
Thomas Witkowski committed
383 384
      stdMpi.startCommunication();
      
385 386 387 388 389 390 391 392 393 394 395 396 397 398
      for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]); 
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	int counter = 0;
	for (; !it.endDofIter(); it.nextDof()) {
	  if (dofMap.count(it.getDofIndex())) {
	    DegreeOfFreedom d = stdMpi.getRecvData(rank)[counter++];
	    if (globalIndex)
	      dofToMatIndex.add(i, dofMap[it.getDofIndex()].global, d);
	    else
	      dofToMatIndex.add(i, it.getDofIndex(), d);
Thomas Witkowski's avatar
Thomas Witkowski committed
399 400 401
	  }
	}
      }
402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443


      // === Communicate DOFs on periodic boundaries. ===

      stdMpi.clear();

      for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, feSpaces[i]); 
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex())) {
	    if (globalIndex) {
	      sendGlobalDofs.push_back(dofMap[it.getDofIndex()].global);
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
	    } else {
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex()));
	    }
	  }
	
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	stdMpi.send(rank, sendGlobalDofs);
	stdMpi.recv(rank);
      }

      stdMpi.startCommunication();

      for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, feSpaces[i]); 
	   !it.end(); it.nextRank()) {

	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);

	int counter = 0;

	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex())) {
	    if (globalIndex) {
444 445 446
	      TEST_EXIT_DBG(counter + 2 <= stdMpi.getRecvData(it.getRank()).size())
		("Should not happen!\n");

447
	      dofToMatIndex.add(i, 
448 449 450
				stdMpi.getRecvData(it.getRank())[counter],
				stdMpi.getRecvData(it.getRank())[counter + 1]);
	      counter += 2;
451 452 453 454 455 456
	    } else {
	      dofToMatIndex.add(i, it.getDofIndex(),
				stdMpi.getRecvData(it.getRank())[counter++]);
	    }
	  }
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
457 458
    }
  }
459 460 461 462 463 464 465 466 467


  void ParallelDofMapping::createIndexSet(IS &is, 
					  int firstComponent, 
					  int nComponents)
  {
    FUNCNAME("ParallelDofMapping::createIndexSet()");

    int firstRankDof = -1;
468 469
    ComponentDofMap &compMap = (*data)[firstComponent];
    DofMap &dofMap = compMap.getMap();
470
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
471
      if (compMap.isRankDof(it->first)) {
472 473 474 475 476 477 478 479 480 481 482 483 484
	if (needMatIndexFromGlobal)
	  firstRankDof = it->second.global;
	else
	  firstRankDof = it->first;
	break;
      }
    }
    
    TEST_EXIT_DBG(firstRankDof >= 0)("No rank DOF found!\n");    
    int firstMatIndex = dofToMatIndex.get(firstComponent, firstRankDof);

    int nRankRows = 0;
    for (int i = firstComponent; i < firstComponent + nComponents; i++)
485
      nRankRows += (*data)[i].nRankDofs;   
486 487 488 489

    ISCreateStride(mpiComm, nRankRows, firstMatIndex, 1, &is);
  }

490
}