ParallelDofMapping.cc 13.6 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"
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
14
#include "parallel/StdMpi.h"
15
16
17
18
19

namespace AMDiS {

  using namespace std;

20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

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

    for (map<int, map<DegreeOfFreedom, int> >::iterator it0 = data.begin();
	 it0 != data.end(); ++it0)
      for (map<DegreeOfFreedom, int>::iterator it1 = it0->second.begin();
	   it1 != it0->second.end(); ++it1)
	if (it1->second == rowIndex) {
	  component = it0->first;
	  dofIndex = it1->first;
	  return;
	}
  
    component = -1;
    dofIndex = -1;
  }


40
  void FeSpaceDofMap::clear()
41
42
  {
    dofMap.clear();
43

44
    nonRankDofs.clear();
45

46
47
48
49
    nRankDofs = 0;
    nLocalDofs = 0;
    nOverallDofs = 0;
    rStartDofs = 0;
50
51
52
  }


53
  void FeSpaceDofMap::update()
54
  {
55
    FUNCNAME("FeSpaceDofMap::update()");
56

57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
    // === 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();
73
      
74
      if (isNonLocal)
75
	computeNonLocalIndices();
76
    }
77
78
79
  }


80
  void FeSpaceDofMap::computeGlobalMapping()
81
  {
82
83
    FUNCNAME("FeSpaceDofMap::computeGlobalMapping()");

84
85
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
      it->second.global = it->second.local + rStartDofs;
86
87
  }

88

89
  void FeSpaceDofMap::computeNonLocalIndices()
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
90
  {
91
    FUNCNAME("FeSpaceDofMap::computeNonLocalIndices()");
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
92

93
94
    TEST_EXIT_DBG(dofComm)("No DOF communicator defined!\n");

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

97
98
99
    // === Send all global indices of DOFs that are owned by the rank to all ===
    // === other ranks that also include this DOF.                           ===

100
    StdMpi<vector<int> > stdMpi(mpiComm);
101

102
    for (DofComm::Iterator it(dofComm->getSendDofs(), 0, feSpace); 
103
104
105
106
	 !it.end(); it.nextRank()) {
      int rank = it.getRank();
      if (meshLevel > 0)
	rank = levelData->mapRank(rank, 0, meshLevel);
107
      
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
108
      for (; !it.endDofIter(); it.nextDof())
109
110
	if (dofMap.count(it.getDofIndex()) && 
	    !nonRankDofs.count(it.getDofIndex()))
111
	  stdMpi.getSendData(rank).
112
	    push_back(dofMap[it.getDofIndex()].global);
113
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
114
115
116

    stdMpi.updateSendDataSize();

117
118
119

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

120
    for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpace); 
121
	 !it.end(); it.nextRank()) {
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
122
123
      bool recvFromRank = false;
      for (; !it.endDofIter(); it.nextDof()) {
124
	if (nonRankDofs.count(it.getDofIndex())) {
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
125
126
127
128
129
	  recvFromRank = true;
	  break;
	}
      }

130
131
132
133
134
135
136
      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
137
138
    }

139
140
141

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

Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
142
143
    stdMpi.startCommunication();

144
145
146

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

Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
153
154
      int i = 0;
      for (; !it.endDofIter(); it.nextDof())
155
156
	if (nonRankDofs.count(it.getDofIndex()))
	  dofMap[it.getDofIndex()].global = stdMpi.getRecvData(rank)[i++];
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
157
158
159
160
    }
  }


161
  void ParallelDofMapping::init(MeshLevelData &ldata,
162
				vector<const FiniteElemSpace*> &fe,
163
				vector<const FiniteElemSpace*> &uniqueFe,
164
				bool b)
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
165
  {
166
167
    FUNCNAME("ParallelDofMapping::init()");

168
    levelData = &ldata;
169
    feSpaces = fe;
170
    feSpacesUnique = uniqueFe;
171
    isNonLocal = b;
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
172
    
173
174
175

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

176
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
177
178
	 it != feSpacesUnique.end(); ++it) {
      addFeSpace(*it);
179
180
      data[*it].setNeedGlobalMapping(isNonLocal);
      data[*it].setNonLocal(isNonLocal);
181
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
182
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
183
184


185
186
187
188
  void ParallelDofMapping::clear()
  {
    FUNCNAME("ParallelDofMapping::clear()");

189
190
    TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n");

191
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
192
193
194
	 it != feSpacesUnique.end(); ++it)
      data[*it].clear();

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

202
203
204
205
206
207
208
209
210
211
212
  
  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);
  }

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
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
222
	 it != feSpacesUnique.end(); ++it)
223
      data[*it].setDofComm(dc);
224
  }
225

226

Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
229
230
231
232
233
  void ParallelDofMapping::addFeSpace(const FiniteElemSpace* feSpace)
  {
    FUNCNAME("ParallelDofMapping::addFeSpace()");
    
    if (data.count(feSpace))
      data.find(feSpace)->second.clear();
    else
234
      data.insert(make_pair(feSpace, FeSpaceDofMap(levelData)));
Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
237
238
239
    
    data.find(feSpace)->second.setFeSpace(feSpace);
  }    


240
  int ParallelDofMapping::computeRankDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
241
  {
242
    FUNCNAME("ParallelDofMapping::computeRankDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
243
244
    
    int result = 0;
245
246
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
247
      result += data[feSpaces[i]].nRankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
248
249
250
251
252
253
    }
    
    return result;
  }


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


268
  int ParallelDofMapping::computeOverallDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
269
  {
270
    FUNCNAME("ParallelDofMapping::computeOverallDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
271
272
273
274

    int result = 0;
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n");
275
      result += data.find(feSpaces[i])->second.nOverallDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
278
279
280
281
    }
    
    return result;
  }


282
  int ParallelDofMapping::computeStartDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
283
  {
284
    FUNCNAME("ParallelDofMapping::computeStartDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
287
288

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


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

300
    // First, update all FE space DOF mappings.
301
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
304
	 it != feSpacesUnique.end(); ++it)
      data[*it].update();
    
305
306
307
308
309
310
311
    // Compute all numbers from this mappings.
    nRankDofs = computeRankDofs();
    nLocalDofs = computeLocalDofs();
    nOverallDofs = computeOverallDofs();
    rStartDofs = computeStartDofs();
    
    // And finally, compute the matrix indices.
312
    computeMatIndex(needMatIndexFromGlobal);
Thomas Witkowski's avatar
Thomas Witkowski committed
313
314
  }
  
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329

  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
330
  
331
  void ParallelDofMapping::computeMatIndex(bool globalIndex)
Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
  {
    FUNCNAME("ParallelDofMapping::computeMatIndex()");
334

335
    dofToMatIndex.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    
337
338
339
    // 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.
340
    int offset = rStartDofs;
341
342
343

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

Thomas Witkowski's avatar
Thomas Witkowski committed
344
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
345
346
347

      // Traverse all DOFs of the FE space and create for all rank owned DOFs
      // a matrix index.
348
      DofMap& dofMap = data[feSpaces[i]].getMap();
349
      for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
350
	if (data[feSpaces[i]].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
351
	  int globalMatIndex = it->second.local + offset;
352
	  if (globalIndex)
353
	    dofToMatIndex.add(i, it->second.global, globalMatIndex);
354
	  else
355
	    dofToMatIndex.add(i, it->first, globalMatIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
356
357
	}
      }
358

359
360
      // Increase the offset for the next FE space by the number of DOFs owned 
      // by the rank in the current FE space.
361
      offset += data[feSpaces[i]].nRankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
362
	
363
      // If there are no non local DOFs, continue with the next FE space.
364
      if (!isNonLocal)
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
	continue;
      
367
      TEST_EXIT_DBG(dofComm != NULL)("No communicator given!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
368
      
369
370
371
      // === Communicate the matrix indices for all DOFs that are on some ===
      // === interior boundaries.                                         ===

372
      StdMpi<vector<DegreeOfFreedom> > stdMpi(mpiComm);
373
      for (DofComm::Iterator it(dofComm->getSendDofs(), 0, feSpaces[i]); 
Thomas Witkowski's avatar
Thomas Witkowski committed
374
375
376
377
378
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex()))
379
	    if (globalIndex)
380
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
381
	    else
382
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex()));
Thomas Witkowski's avatar
Thomas Witkowski committed
383
	
384
385
386
387
388
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	stdMpi.send(rank, sendGlobalDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
      }
      
391
      for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]); 
392
393
394
395
396
397
398
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);

	stdMpi.recv(rank);
      }
399

Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
      stdMpi.startCommunication();
      
402
403
404
405
406
407
408
409
410
411
412
413
414
415
      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
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
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460


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

464
	      dofToMatIndex.add(i, 
465
466
467
				stdMpi.getRecvData(it.getRank())[counter],
				stdMpi.getRecvData(it.getRank())[counter + 1]);
	      counter += 2;
468
469
470
471
472
473
	    } else {
	      dofToMatIndex.add(i, it.getDofIndex(),
				stdMpi.getRecvData(it.getRank())[counter++]);
	    }
	  }
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
474
475
    }
  }
476
477
478
479
480
481
482
483
484
485
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


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

513
}