ParallelDofMapping.cc 14.1 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
    : levelData(ld),
      dofComm(NULL),
      feSpace(NULL),
45
      globalMapping(false)
46
47
48
49
50
  {
    clear();
  }


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

55
    nonRankDofs.clear();
56

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


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

68
69
70
71
72
73
74
75
76
77
78
79
80
81
    // === 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. ===
    
82
83
84
    if (globalMapping) {
      computeGlobalMapping();      
      computeNonLocalIndices();
85
    }
86
87
88
  }


89
  void ComponentDofMap::computeGlobalMapping()
90
  {
91
    FUNCNAME("ComponentDofMap::computeGlobalMapping()");
92

93
94
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
      it->second.global = it->second.local + rStartDofs;
95
96
  }

97

98
  void ComponentDofMap::computeNonLocalIndices()
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
99
  {
100
    FUNCNAME("ComponentDofMap::computeNonLocalIndices()");
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
101

102
103
    TEST_EXIT_DBG(dofComm)("No DOF communicator defined!\n");

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

106
107
108
    // === Send all global indices of DOFs that are owned by the rank to all ===
    // === other ranks that also include this DOF.                           ===

109
    StdMpi<vector<int> > stdMpi(mpiComm);
110

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

    stdMpi.updateSendDataSize();

126
127
128

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

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

139
140
141
142
143
144
145
      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
146
147
    }

148
149
150

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

Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
151
152
    stdMpi.startCommunication();

153
154
155

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

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


170
171
172
173
  void FeSpaceData::init(vector<const FiniteElemSpace*> &f0,
			 vector<const FiniteElemSpace*> &f1,
			 bool globalMapping,
			 MeshLevelData &levelData)
174
  {
175
176
177
178
179
180
181
182
183
184
185
    componentSpaces = f0;
    feSpaces = f1;
    for (vector<const FiniteElemSpace*>::iterator it = feSpaces.begin();
     	 it != feSpaces.end(); ++it) {
      if (componentData.count(*it))
	componentData.find(*it)->second.clear();
      else
	componentData.insert(make_pair(*it, ComponentDofMap(&levelData)));   

      componentData[*it].setFeSpace(*it);    
      componentData[*it].setGlobalMapping(globalMapping);
186
187
188
189
    }
  }


190
191
192
193
  void ComponentData::init(vector<const FiniteElemSpace*> &f0,
			   vector<const FiniteElemSpace*> &f1,
			   bool globalMapping,
			   MeshLevelData &levelData)
194
  {
195
196
197
198
199
200
201
202
203
204
205
    componentSpaces = f0;
    feSpaces = f1;

    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
      if (componentData.count(component)) 
	componentData.find(component)->second.clear();
      else
	componentData.insert(make_pair(component, ComponentDofMap(&levelData)));
      
      componentData[component].setFeSpace(componentSpaces[component]);
      componentData[component].setGlobalMapping(globalMapping);
206
207
208
209
    }    
  }


210
211
212
  ParallelDofMapping::ParallelDofMapping(DofMappingMode mode) 
    : levelData(NULL),
      dofComm(NULL),
213
      globalMapping(true),
214
215
216
217
218
219
220
221
      needMatIndexFromGlobal(false),
      nRankDofs(1),
      nLocalDofs(1),
      nOverallDofs(1),
      rStartDofs(1)
  {
    switch (mode) {
    case COMPONENT_WISE:
222
      data = new ComponentData();
223
224
      break;
    case FESPACE_WISE:
225
      data = new FeSpaceData();
226
227
228
229
230
231
232
233
234
235
      break;
    }
    
    nRankDofs = -1;
    nLocalDofs = -1;
    nOverallDofs = -1;
    rStartDofs = -1;
  } 


236
  void ParallelDofMapping::init(MeshLevelData &ldata,
237
				vector<const FiniteElemSpace*> &fe,
238
				vector<const FiniteElemSpace*> &uniqueFe,
239
				bool b)
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
240
  {
241
242
    FUNCNAME("ParallelDofMapping::init()");

243
    levelData = &ldata;
244
245
    globalMapping = b;
    componentSpaces = fe;
246
    feSpaces = uniqueFe;
247

248
    data->init(fe, uniqueFe, globalMapping, ldata);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
249
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
250
251


252
253
254
255
  void ParallelDofMapping::clear()
  {
    FUNCNAME("ParallelDofMapping::clear()");

256
257
    TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n");

258
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
Thomas Witkowski's avatar
Thomas Witkowski committed
259
      it->clear();
260

261
262
263
264
265
    nRankDofs = -1;
    nLocalDofs = -1;
    nOverallDofs = -1;
    rStartDofs = -1;
    dofToMatIndex.clear();
266
267
  }

268
269
270
271
272
273
  
  void ParallelDofMapping::setMpiComm(MPI::Intracomm &m, int l)
  {
    mpiComm = m;
    meshLevel = l;

274
275
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->setMpiComm(m, l);
276
277
  }

278

279
  void ParallelDofMapping::setDofComm(DofComm &dc)
280
281
282
  {
    FUNCNAME("ParallelDofMapping::setDofComm()");

283
    dofComm = &dc;
284
285

    // Add the DOF communicator also to all FE space DOF mappings.
286
287
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->setDofComm(dc);
288
  }
289

290

291
  int ParallelDofMapping::computeRankDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
292
  {
293
    FUNCNAME("ParallelDofMapping::computeRankDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
294
295
    
    int result = 0;
296
297
298
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nRankDofs;
   
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
301
302
    return result;
  }


303
  int ParallelDofMapping::computeLocalDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
304
  {
305
    FUNCNAME("ParallelDofMapping::computeLocalDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
306
307
    
    int result = 0;
308
309
310
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nLocalDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
311
312
313
314
    return result;
  }


315
  int ParallelDofMapping::computeOverallDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
316
  {
317
    FUNCNAME("ParallelDofMapping::computeOverallDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
318
319

    int result = 0;
320
321
322
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nOverallDofs;
       
Thomas Witkowski's avatar
Thomas Witkowski committed
323
324
325
326
    return result;
  }


327
  int ParallelDofMapping::computeStartDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
328
  {
329
    FUNCNAME("ParallelDofMapping::computeStartDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331

    int result = 0;
332
333
334
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->rStartDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
337
338
339
340
341
342
    return result;
  }


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

343
    // First, update all FE space DOF mappings.
344
345
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->update();
Thomas Witkowski's avatar
Thomas Witkowski committed
346
    
347
348
349
350
351
    // Compute all numbers from this mappings.
    nRankDofs = computeRankDofs();
    nLocalDofs = computeLocalDofs();
    nOverallDofs = computeOverallDofs();
    rStartDofs = computeStartDofs();
352

353
    // And finally, compute the matrix indices.
354
    computeMatIndex(needMatIndexFromGlobal);
Thomas Witkowski's avatar
Thomas Witkowski committed
355
356
  }
  
357

358
  void ParallelDofMapping::computeMatIndex(bool globalIndex)
Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
  {
    FUNCNAME("ParallelDofMapping::computeMatIndex()");
361

362
    dofToMatIndex.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
363
    
364
365
366
    // 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.
367
    int offset = rStartDofs;
368
369
370

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

371
372
    for (unsigned int component = 0; component < componentSpaces.size(); 
	 component++) {
373
374
375

      // Traverse all DOFs of the FE space and create for all rank owned DOFs
      // a matrix index.
376
      DofMap& dofMap = (*data)[component].getMap();
377
      for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
378
	if ((*data)[component].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
379
	  int globalMatIndex = it->second.local + offset;
380
	  if (globalIndex)
381
	    dofToMatIndex.add(component, it->second.global, globalMatIndex);
382
	  else
383
	    dofToMatIndex.add(component, it->first, globalMatIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
384
385
	}
      }
386

387
388
      // Increase the offset for the next FE space by the number of DOFs owned 
      // by the rank in the current FE space.
389
      offset += (*data)[component].nRankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
390
	
391
      // If there are no non local DOFs, continue with the next FE space.
392
      if (!globalMapping)
Thomas Witkowski's avatar
Thomas Witkowski committed
393
394
	continue;
      
395
      TEST_EXIT_DBG(dofComm != NULL)("No communicator given!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
396
      
397
398
399
      // === Communicate the matrix indices for all DOFs that are on some ===
      // === interior boundaries.                                         ===

400
      StdMpi<vector<DegreeOfFreedom> > stdMpi(mpiComm);
401
      for (DofComm::Iterator it(dofComm->getSendDofs(), 0, componentSpaces[component]); 
Thomas Witkowski's avatar
Thomas Witkowski committed
402
403
404
405
406
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex()))
407
	    if (globalIndex)
408
409
	      sendGlobalDofs.push_back(dofToMatIndex.get(component, 
							 dofMap[it.getDofIndex()].global));
410
	    else
411
	      sendGlobalDofs.push_back(dofToMatIndex.get(component, it.getDofIndex()));
Thomas Witkowski's avatar
Thomas Witkowski committed
412
	
413
414
415
416
417
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	stdMpi.send(rank, sendGlobalDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
418
419
      }
      
420
      for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, componentSpaces[component]); 
421
422
423
424
425
426
427
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);

	stdMpi.recv(rank);
      }
428

Thomas Witkowski's avatar
Thomas Witkowski committed
429
430
      stdMpi.startCommunication();
      
431
      for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, componentSpaces[component]); 
432
433
434
435
436
437
438
439
440
441
	   !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)
442
	      dofToMatIndex.add(component, dofMap[it.getDofIndex()].global, d);
443
	    else
444
	      dofToMatIndex.add(component, it.getDofIndex(), d);
Thomas Witkowski's avatar
Thomas Witkowski committed
445
446
447
	  }
	}
      }
448
449
450
451
452
453


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

      stdMpi.clear();

454
      for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, componentSpaces[component]); 
455
456
457
458
459
460
461
	   !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);
462
463
	      sendGlobalDofs.push_back(dofToMatIndex.get(component, 
							 dofMap[it.getDofIndex()].global));
464
	    } else {
465
	      sendGlobalDofs.push_back(dofToMatIndex.get(component, it.getDofIndex()));
466
467
468
469
470
471
472
473
474
475
476
477
478
	    }
	  }
	
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	stdMpi.send(rank, sendGlobalDofs);
	stdMpi.recv(rank);
      }

      stdMpi.startCommunication();

479
      for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, componentSpaces[component]); 
480
481
482
483
484
485
486
487
488
489
490
	   !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) {
491
492
493
	      TEST_EXIT_DBG(counter + 2 <= stdMpi.getRecvData(it.getRank()).size())
		("Should not happen!\n");

494
	      dofToMatIndex.add(component, 
495
496
497
				stdMpi.getRecvData(it.getRank())[counter],
				stdMpi.getRecvData(it.getRank())[counter + 1]);
	      counter += 2;
498
	    } else {
499
500
	      dofToMatIndex.add(component, 
				it.getDofIndex(),
501
502
503
504
				stdMpi.getRecvData(it.getRank())[counter++]);
	    }
	  }
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
505
506
    }
  }
507
508
509
510
511
512
513
514
515


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

    int firstRankDof = -1;
516
517
    ComponentDofMap &compMap = (*data)[firstComponent];
    DofMap &dofMap = compMap.getMap();
518
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
519
      if (compMap.isRankDof(it->first)) {
520
521
522
523
524
525
526
527
528
529
530
531
532
	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++)
533
      nRankRows += (*data)[i].nRankDofs;   
534
535
536
537

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

538
}