ParallelDofMapping.cc 12.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
  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
194
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      data->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
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
326
327
328

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

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

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

	stdMpi.recv(rank);
      }
380

Thomas Witkowski's avatar
Thomas Witkowski committed
381
382
      stdMpi.startCommunication();
      
383
384
385
386
387
388
389
390
391
392
393
394
395
396
      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
397
398
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


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

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


  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;
472
473
    ComponentDofMap &compMap = (*data)[firstComponent];
    DofMap &dofMap = compMap.getMap();
474
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
475
      if (compMap.isRankDof(it->first)) {
476
477
478
479
480
481
482
483
484
485
486
487
488
	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++)
489
      nRankRows += (*data)[i].nRankDofs;   
490
491
492
493

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

494
}