ParallelDofMapping.cc 13.8 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
42
43
44
45
46
47
48
49
50
51
  FeSpaceDofMap::FeSpaceDofMap(MeshLevelData* ld)
    : levelData(ld),
      dofComm(NULL),
      feSpace(NULL),
      needGlobalMapping(false),
      isNonLocal(false)
  {
    clear();
  }


52
  void FeSpaceDofMap::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 FeSpaceDofMap::update()
66
  {
67
    FUNCNAME("FeSpaceDofMap::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 FeSpaceDofMap::computeGlobalMapping()
93
  {
94
95
    FUNCNAME("FeSpaceDofMap::computeGlobalMapping()");

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

100

101
  void FeSpaceDofMap::computeNonLocalIndices()
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
102
  {
103
    FUNCNAME("FeSpaceDofMap::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
    feSpaces = fe;
182
    feSpacesUnique = uniqueFe;
183
    isNonLocal = b;
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
184
    
185
186
187

    // === Init the mapping for all different FE spaces. ===

188
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
189
190
	 it != feSpacesUnique.end(); ++it) {
      addFeSpace(*it);
191
192
      data[*it].setNeedGlobalMapping(isNonLocal);
      data[*it].setNonLocal(isNonLocal);
193
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
194
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
195
196


197
198
199
200
  void ParallelDofMapping::clear()
  {
    FUNCNAME("ParallelDofMapping::clear()");

201
202
    TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n");

203
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
204
205
206
	 it != feSpacesUnique.end(); ++it)
      data[*it].clear();

207
208
209
210
211
    nRankDofs = -1;
    nLocalDofs = -1;
    nOverallDofs = -1;
    rStartDofs = -1;
    dofToMatIndex.clear();
212
213
  }

214
215
216
217
218
219
220
221
222
223
224
  
  void ParallelDofMapping::setMpiComm(MPI::Intracomm &m, int l)
  {
    mpiComm = m;
    meshLevel = l;

    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
	 it != feSpacesUnique.end(); ++it)
      data[*it].setMpiComm(m, l);
  }

225

226
  void ParallelDofMapping::setDofComm(DofComm &dc)
227
228
229
  {
    FUNCNAME("ParallelDofMapping::setDofComm()");

230
    dofComm = &dc;
231
232

    // Add the DOF communicator also to all FE space DOF mappings.
233
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
234
	 it != feSpacesUnique.end(); ++it)
235
      data[*it].setDofComm(dc);
236
  }
237

238

Thomas Witkowski's avatar
Thomas Witkowski committed
239
240
241
242
243
244
245
  void ParallelDofMapping::addFeSpace(const FiniteElemSpace* feSpace)
  {
    FUNCNAME("ParallelDofMapping::addFeSpace()");
    
    if (data.count(feSpace))
      data.find(feSpace)->second.clear();
    else
246
      data.insert(make_pair(feSpace, FeSpaceDofMap(levelData)));
Thomas Witkowski's avatar
Thomas Witkowski committed
247
248
249
250
251
    
    data.find(feSpace)->second.setFeSpace(feSpace);
  }    


252
  int ParallelDofMapping::computeRankDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
253
  {
254
    FUNCNAME("ParallelDofMapping::computeRankDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
255
256
    
    int result = 0;
257
258
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
259
      result += data[feSpaces[i]].nRankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
260
261
262
263
264
265
    }
    
    return result;
  }


266
  int ParallelDofMapping::computeLocalDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
267
  {
268
    FUNCNAME("ParallelDofMapping::computeLocalDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270
    
    int result = 0;
271
272
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
273
      result += data[feSpaces[i]].nLocalDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
276
277
278
279
    }
    
    return result;
  }


280
  int ParallelDofMapping::computeOverallDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
281
  {
282
    FUNCNAME("ParallelDofMapping::computeOverallDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
283
284
285
286

    int result = 0;
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
287
      result += data.find(feSpaces[i])->second.nOverallDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
290
291
292
293
    }
    
    return result;
  }


294
  int ParallelDofMapping::computeStartDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
295
  {
296
    FUNCNAME("ParallelDofMapping::computeStartDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
297
298
299
300

    int result = 0;
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
301
      result += data.find(feSpaces[i])->second.rStartDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
304
305
306
307
308
309
310
311
    }
    
    return result;
  }


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

312
    // First, update all FE space DOF mappings.
313
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
316
	 it != feSpacesUnique.end(); ++it)
      data[*it].update();
    
317
318
319
320
321
322
323
    // Compute all numbers from this mappings.
    nRankDofs = computeRankDofs();
    nLocalDofs = computeLocalDofs();
    nOverallDofs = computeOverallDofs();
    rStartDofs = computeStartDofs();
    
    // And finally, compute the matrix indices.
324
    computeMatIndex(needMatIndexFromGlobal);
Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
  }
  
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
342
  
343
  void ParallelDofMapping::computeMatIndex(bool globalIndex)
Thomas Witkowski's avatar
Thomas Witkowski committed
344
345
  {
    FUNCNAME("ParallelDofMapping::computeMatIndex()");
346

347
    dofToMatIndex.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
348
    
349
350
351
    // 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.
352
    int offset = rStartDofs;
353
354
355

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

Thomas Witkowski's avatar
Thomas Witkowski committed
356
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
357
358
359

      // Traverse all DOFs of the FE space and create for all rank owned DOFs
      // a matrix index.
360
      DofMap& dofMap = data[feSpaces[i]].getMap();
361
      for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
362
	if (data[feSpaces[i]].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
363
	  int globalMatIndex = it->second.local + offset;
364
	  if (globalIndex)
365
	    dofToMatIndex.add(i, it->second.global, globalMatIndex);
366
	  else
367
	    dofToMatIndex.add(i, it->first, globalMatIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
368
369
	}
      }
370

371
372
      // Increase the offset for the next FE space by the number of DOFs owned 
      // by the rank in the current FE space.
373
      offset += data[feSpaces[i]].nRankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
374
	
375
      // If there are no non local DOFs, continue with the next FE space.
376
      if (!isNonLocal)
Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
	continue;
      
379
      TEST_EXIT_DBG(dofComm != NULL)("No communicator given!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
380
      
381
382
383
      // === Communicate the matrix indices for all DOFs that are on some ===
      // === interior boundaries.                                         ===

384
      StdMpi<vector<DegreeOfFreedom> > stdMpi(mpiComm);
385
      for (DofComm::Iterator it(dofComm->getSendDofs(), 0, feSpaces[i]); 
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
389
390
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex()))
391
	    if (globalIndex)
392
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
393
	    else
394
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex()));
Thomas Witkowski's avatar
Thomas Witkowski committed
395
	
396
397
398
399
400
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	stdMpi.send(rank, sendGlobalDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
401
402
      }
      
403
      for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]); 
404
405
406
407
408
409
410
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);

	stdMpi.recv(rank);
      }
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
      stdMpi.startCommunication();
      
414
415
416
417
418
419
420
421
422
423
424
425
426
427
      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
428
429
430
	  }
	}
      }
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472


      // === 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) {
473
474
475
	      TEST_EXIT_DBG(counter + 2 <= stdMpi.getRecvData(it.getRank()).size())
		("Should not happen!\n");

476
	      dofToMatIndex.add(i, 
477
478
479
				stdMpi.getRecvData(it.getRank())[counter],
				stdMpi.getRecvData(it.getRank())[counter + 1]);
	      counter += 2;
480
481
482
483
484
485
	    } else {
	      dofToMatIndex.add(i, it.getDofIndex(),
				stdMpi.getRecvData(it.getRank())[counter++]);
	    }
	  }
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
486
487
    }
  }
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524


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

    TEST_EXIT_DBG(firstComponent + nComponents <= feSpaces.size())
      ("Should not happen!\n");

    TEST_EXIT_DBG(data.count(feSpaces[firstComponent]))
      ("No data for FE space at address %p!\n", feSpaces[firstComponent]);
    
    int firstRankDof = -1;
    FeSpaceDofMap &feMap = data.find(feSpaces[firstComponent])->second;
    DofMap &dofMap = feMap.getMap();
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
      if (feMap.isRankDof(it->first)) {
	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++)
      nRankRows += data.find(feSpaces[firstComponent])->second.nRankDofs;

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

525
}