ParallelDofMapping.cc 14.9 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
178
179
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
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
  void ComponentDataEqFeSpace::init(vector<const FiniteElemSpace*> &f0,
				    vector<const FiniteElemSpace*> &f1,
				    bool isNonLocal,
				    MeshLevelData &levelData)
  {
    feSpaces = f1;
    feSpacesUnique = f0;
    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)
  {
    feSpaces = f1;
    feSpacesUnique = f0;

    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
370
371
    // Compute all numbers from this mappings.
    nRankDofs = computeRankDofs();
    nLocalDofs = computeLocalDofs();
    nOverallDofs = computeOverallDofs();
    rStartDofs = computeStartDofs();
    
    // And finally, compute the matrix indices.
372
    computeMatIndex(needMatIndexFromGlobal);
Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
  }
  
375
376
377
378
379

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

380
381
382
383
384
385
386
387
388
389
    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;

390
391
392
393

    update();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
394
  
395
  void ParallelDofMapping::computeMatIndex(bool globalIndex)
Thomas Witkowski's avatar
Thomas Witkowski committed
396
397
  {
    FUNCNAME("ParallelDofMapping::computeMatIndex()");
398

399
    dofToMatIndex.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
400
    
401
402
403
    // 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.
404
    int offset = rStartDofs;
405
406
407

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
410
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
411
412
413

      // Traverse all DOFs of the FE space and create for all rank owned DOFs
      // a matrix index.
414
      DofMap& dofMap = (*data)[i].getMap();
415
      for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
416
	if ((*data)[i].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
417
	  int globalMatIndex = it->second.local + offset;
418
	  if (globalIndex)
419
	    dofToMatIndex.add(i, it->second.global, globalMatIndex);
420
	  else
421
	    dofToMatIndex.add(i, it->first, globalMatIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
	}
      }
424

425
426
      // Increase the offset for the next FE space by the number of DOFs owned 
      // by the rank in the current FE space.
427
      offset += (*data)[i].nRankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
428
	
429
      // If there are no non local DOFs, continue with the next FE space.
430
      if (!isNonLocal)
Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
	continue;
      
433
      TEST_EXIT_DBG(dofComm != NULL)("No communicator given!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
434
      
435
436
437
      // === Communicate the matrix indices for all DOFs that are on some ===
      // === interior boundaries.                                         ===

438
      StdMpi<vector<DegreeOfFreedom> > stdMpi(mpiComm);
439
      for (DofComm::Iterator it(dofComm->getSendDofs(), 0, feSpaces[i]); 
Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
442
443
444
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex()))
445
	    if (globalIndex)
446
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
447
	    else
448
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, it.getDofIndex()));
Thomas Witkowski's avatar
Thomas Witkowski committed
449
	
450
451
452
453
454
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	stdMpi.send(rank, sendGlobalDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
455
456
      }
      
457
      for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, feSpaces[i]); 
458
459
460
461
462
463
464
	   !it.end(); it.nextRank()) {
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);

	stdMpi.recv(rank);
      }
465

Thomas Witkowski's avatar
Thomas Witkowski committed
466
467
      stdMpi.startCommunication();
      
468
469
470
471
472
473
474
475
476
477
478
479
480
481
      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
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
513
514
515
516
517
518
519
520
521
522
523
524
525
526


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

530
	      dofToMatIndex.add(i, 
531
532
533
				stdMpi.getRecvData(it.getRank())[counter],
				stdMpi.getRecvData(it.getRank())[counter + 1]);
	      counter += 2;
534
535
536
537
538
539
	    } else {
	      dofToMatIndex.add(i, it.getDofIndex(),
				stdMpi.getRecvData(it.getRank())[counter++]);
	    }
	  }
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
540
541
    }
  }
542
543
544
545
546
547
548
549
550


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

    int firstRankDof = -1;
551
552
    ComponentDofMap &compMap = (*data)[firstComponent];
    DofMap &dofMap = compMap.getMap();
553
    for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
554
      if (compMap.isRankDof(it->first)) {
555
556
557
558
559
560
561
562
563
564
565
566
567
	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++)
568
      nRankRows += (*data)[i].nRankDofs;   
569
570
571
572

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

573
}