ParallelDofMapping.cc 14.5 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
174
175
176
177
  void ComponentDataEqFeSpace::init(vector<const FiniteElemSpace*> &f0,
				    vector<const FiniteElemSpace*> &f1,
				    bool isNonLocal,
				    MeshLevelData &levelData)
  {
178
179
    feSpaces = f0;
    feSpacesUnique = f1;
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
    for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin();
     	 it != feSpacesUnique.end(); ++it) {
      addFeSpace(*it, levelData);
      componentData[*it].setNeedGlobalMapping(isNonLocal);
      componentData[*it].setNonLocal(isNonLocal);
    }
  }


  void ComponentDataEqFeSpace::addFeSpace(const FiniteElemSpace* feSpace,
					  MeshLevelData &levelData)
  {
    if (componentData.count(feSpace))
      componentData.find(feSpace)->second.clear();
    else
      componentData.insert(make_pair(feSpace, ComponentDofMap(&levelData)));
    
    componentData.find(feSpace)->second.setFeSpace(feSpace);    
  }


  void ComponentDataDiffFeSpace::init(vector<const FiniteElemSpace*> &f0,
				      vector<const FiniteElemSpace*> &f1,
				      bool isNonLocal,
				      MeshLevelData &levelData)
  {
206
207
    feSpaces = f0;
    feSpacesUnique = f1;
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229

    for (unsigned int component = 0; component < feSpaces.size(); component++) {
      addComponent(component, feSpaces[component], levelData);
      componentData[component].setNeedGlobalMapping(isNonLocal);
      componentData[component].setNonLocal(isNonLocal);
    }    
  }


  void ComponentDataDiffFeSpace::addComponent(unsigned int component,
					      const FiniteElemSpace* feSpace,
					      MeshLevelData &levelData)
  {
    if (componentData.count(component)) 
      componentData.find(component)->second.clear();
    else
      componentData.insert(make_pair(component, ComponentDofMap(&levelData)));

    componentData.find(component)->second.setFeSpace(feSpace);
  }


230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
  ParallelDofMapping::ParallelDofMapping(DofMappingMode mode) 
    : levelData(NULL),
      dofComm(NULL),
      isNonLocal(true),
      needMatIndexFromGlobal(false),
      nRankDofs(1),
      nLocalDofs(1),
      nOverallDofs(1),
      rStartDofs(1)
  {
    switch (mode) {
    case COMPONENT_WISE:
      data = new ComponentDataDiffFeSpace();
      break;
    case FESPACE_WISE:
      data = new ComponentDataEqFeSpace();
      break;
    }
    
    nRankDofs = -1;
    nLocalDofs = -1;
    nOverallDofs = -1;
    rStartDofs = -1;
  } 


256
  void ParallelDofMapping::init(MeshLevelData &ldata,
257
				vector<const FiniteElemSpace*> &fe,
258
				vector<const FiniteElemSpace*> &uniqueFe,
259
				bool b)
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
260
  {
261
262
    FUNCNAME("ParallelDofMapping::init()");

263
    levelData = &ldata;
264
    isNonLocal = b;
265

266
    data->init(fe, uniqueFe, isNonLocal, ldata);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
267
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
268
269


270
271
272
273
  void ParallelDofMapping::clear()
  {
    FUNCNAME("ParallelDofMapping::clear()");

274
275
    TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n");

276
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
Thomas Witkowski's avatar
Thomas Witkowski committed
277
      it->clear();
278

279
280
281
282
283
    nRankDofs = -1;
    nLocalDofs = -1;
    nOverallDofs = -1;
    rStartDofs = -1;
    dofToMatIndex.clear();
284
285
  }

286
287
288
289
290
291
  
  void ParallelDofMapping::setMpiComm(MPI::Intracomm &m, int l)
  {
    mpiComm = m;
    meshLevel = l;

292
293
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->setMpiComm(m, l);
294
295
  }

296

297
  void ParallelDofMapping::setDofComm(DofComm &dc)
298
299
300
  {
    FUNCNAME("ParallelDofMapping::setDofComm()");

301
    dofComm = &dc;
302
303

    // Add the DOF communicator also to all FE space DOF mappings.
304
305
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->setDofComm(dc);
306
  }
307

308

309
  int ParallelDofMapping::computeRankDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
310
  {
311
    FUNCNAME("ParallelDofMapping::computeRankDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
312
313
    
    int result = 0;
314
315
316
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nRankDofs;
   
Thomas Witkowski's avatar
Thomas Witkowski committed
317
318
319
320
    return result;
  }


321
  int ParallelDofMapping::computeLocalDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
322
  {
323
    FUNCNAME("ParallelDofMapping::computeLocalDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
324
325
    
    int result = 0;
326
327
328
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nLocalDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
331
332
    return result;
  }


333
  int ParallelDofMapping::computeOverallDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
334
  {
335
    FUNCNAME("ParallelDofMapping::computeOverallDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
336
337

    int result = 0;
338
339
340
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nOverallDofs;
       
Thomas Witkowski's avatar
Thomas Witkowski committed
341
342
343
344
    return result;
  }


345
  int ParallelDofMapping::computeStartDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
346
  {
347
    FUNCNAME("ParallelDofMapping::computeStartDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
348
349

    int result = 0;
350
351
352
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->rStartDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
353
354
355
356
357
358
359
360
    return result;
  }


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

361
    // First, update all FE space DOF mappings.
362
363
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->update();
Thomas Witkowski's avatar
Thomas Witkowski committed
364
    
365
366
367
368
369
    // Compute all numbers from this mappings.
    nRankDofs = computeRankDofs();
    nLocalDofs = computeLocalDofs();
    nOverallDofs = computeOverallDofs();
    rStartDofs = computeStartDofs();
370

371
    // And finally, compute the matrix indices.
372
    computeMatIndex(needMatIndexFromGlobal);
Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
  }
  
375

376
  void ParallelDofMapping::computeMatIndex(bool globalIndex)
Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
  {
    FUNCNAME("ParallelDofMapping::computeMatIndex()");
379

380
    dofToMatIndex.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
381
    
382
383
384
    // 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.
385
    int offset = rStartDofs;
386
387
388

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
391
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
392
393
394

      // Traverse all DOFs of the FE space and create for all rank owned DOFs
      // a matrix index.
395
      DofMap& dofMap = (*data)[i].getMap();
396
      for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
397
	if ((*data)[i].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
398
	  int globalMatIndex = it->second.local + offset;
399
	  if (globalIndex)
400
	    dofToMatIndex.add(i, it->second.global, globalMatIndex);
401
	  else
402
	    dofToMatIndex.add(i, it->first, globalMatIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
403
404
	}
      }
405

406
407
      // Increase the offset for the next FE space by the number of DOFs owned 
      // by the rank in the current FE space.
408
      offset += (*data)[i].nRankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
409
	
410
      // If there are no non local DOFs, continue with the next FE space.
411
      if (!isNonLocal)
Thomas Witkowski's avatar
Thomas Witkowski committed
412
413
	continue;
      
414
      TEST_EXIT_DBG(dofComm != NULL)("No communicator given!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
415
      
416
417
418
      // === Communicate the matrix indices for all DOFs that are on some ===
      // === interior boundaries.                                         ===

419
      StdMpi<vector<DegreeOfFreedom> > stdMpi(mpiComm);
420
      for (DofComm::Iterator it(dofComm->getSendDofs(), 0, feSpaces[i]); 
Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
423
424
425
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex()))
426
	    if (globalIndex)
427
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
428
	    else
429
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex()));
Thomas Witkowski's avatar
Thomas Witkowski committed
430
	
431
432
433
434
435
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	stdMpi.send(rank, sendGlobalDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
      }
      
438
      for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]); 
439
440
441
442
443
444
445
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);

	stdMpi.recv(rank);
      }
446

Thomas Witkowski's avatar
Thomas Witkowski committed
447
448
      stdMpi.startCommunication();
      
449
450
451
452
453
454
455
456
457
458
459
460
461
462
      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
463
464
465
	  }
	}
      }
466
467
468
469
470
471
472
473
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


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

511
	      dofToMatIndex.add(i, 
512
513
514
				stdMpi.getRecvData(it.getRank())[counter],
				stdMpi.getRecvData(it.getRank())[counter + 1]);
	      counter += 2;
515
516
517
518
519
520
	    } else {
	      dofToMatIndex.add(i, it.getDofIndex(),
				stdMpi.getRecvData(it.getRank())[counter++]);
	    }
	  }
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
521
522
    }
  }
523
524
525
526
527
528
529
530
531


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

    int firstRankDof = -1;
532
533
    ComponentDofMap &compMap = (*data)[firstComponent];
    DofMap &dofMap = compMap.getMap();
534
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
535
      if (compMap.isRankDof(it->first)) {
536
537
538
539
540
541
542
543
544
545
546
547
548
	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++)
549
      nRankRows += (*data)[i].nRankDofs;   
550
551
552
553

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

554
}