ParallelDofMapping.cc 15.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
    : levelData(ld),
Thomas Witkowski's avatar
Thomas Witkowski committed
43
      meshLevel(0),
44
45
      dofComm(NULL),
      feSpace(NULL),
46
      globalMapping(false)
47
48
49
50
51
  {
    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
    // === 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. ===
    
83
84
85
    if (globalMapping) {
      computeGlobalMapping();      
      computeNonLocalIndices();
86
    }
87
88
89
  }


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

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

98

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

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

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

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

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

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

    stdMpi.updateSendDataSize();

127
128
129

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

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

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

149
150
151

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

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

154
155
156

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

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


171
172
173
174
  void FeSpaceData::init(vector<const FiniteElemSpace*> &f0,
			 vector<const FiniteElemSpace*> &f1,
			 bool globalMapping,
			 MeshLevelData &levelData)
175
  {
176
177
178
179
180
181
182
183
184
185
186
    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);
187
188
189
190
    }
  }


191
192
193
194
  void ComponentData::init(vector<const FiniteElemSpace*> &f0,
			   vector<const FiniteElemSpace*> &f1,
			   bool globalMapping,
			   MeshLevelData &levelData)
195
  {
196
197
198
199
200
201
202
203
204
205
206
    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);
207
208
209
210
    }    
  }


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


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

246
    levelData = &ldata;
247
248
    globalMapping = b;
    componentSpaces = fe;
249
    feSpaces = uniqueFe;
250

251
    data->init(fe, uniqueFe, globalMapping, ldata);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
252
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
253
254


255
256
257
258
  void ParallelDofMapping::clear()
  {
    FUNCNAME("ParallelDofMapping::clear()");

259
260
    TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n");

261
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
Thomas Witkowski's avatar
Thomas Witkowski committed
262
      it->clear();
263

264
265
266
267
268
    nRankDofs = -1;
    nLocalDofs = -1;
    nOverallDofs = -1;
    rStartDofs = -1;
    dofToMatIndex.clear();
269
270
  }

271
272
273
274
275
276
  
  void ParallelDofMapping::setMpiComm(MPI::Intracomm &m, int l)
  {
    mpiComm = m;
    meshLevel = l;

277
278
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->setMpiComm(m, l);
279
280
  }

281

282
  void ParallelDofMapping::setDofComm(DofComm &dc)
283
284
285
  {
    FUNCNAME("ParallelDofMapping::setDofComm()");

286
    dofComm = &dc;
287
288

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

293

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


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

Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
316
317
    return result;
  }


318
  int ParallelDofMapping::computeOverallDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
319
  {
320
    FUNCNAME("ParallelDofMapping::computeOverallDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
321
322

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


330
  int ParallelDofMapping::computeStartDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
331
  {
332
    FUNCNAME("ParallelDofMapping::computeStartDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334

    int result = 0;
335
336
337
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->rStartDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
340
341
342
343
344
345
    return result;
  }


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

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

356
    // And finally, compute the matrix indices.
357
    computeMatIndex(needMatIndexFromGlobal);
Thomas Witkowski's avatar
Thomas Witkowski committed
358
359
  }
  
360

361
362
363
364
365
366
367
368
  void ParallelDofMapping::updateMatIndex()
  {
    FUNCNAME("ParallelDofMapping::updateMatIndex()");

    computeMatIndex(needMatIndexFromGlobal);
  }


369
  void ParallelDofMapping::computeMatIndex(bool globalIndex)
Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
  {
    FUNCNAME("ParallelDofMapping::computeMatIndex()");
372

373
    dofToMatIndex.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
374
    
375
376
377
    // 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.
378
    int offset = rStartDofs;
379
380
381

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

382
383
    for (unsigned int component = 0; component < componentSpaces.size(); 
	 component++) {
384
385
386

      // Traverse all DOFs of the FE space and create for all rank owned DOFs
      // a matrix index.
387
      DofMap& dofMap = (*data)[component].getMap();
388
      for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) {
389
	if ((*data)[component].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
390
	  int globalMatIndex = it->second.local + offset;
391
	  if (globalIndex)
392
	    dofToMatIndex.add(component, it->second.global, globalMatIndex);
393
	  else
394
	    dofToMatIndex.add(component, it->first, globalMatIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
	}
      }
397

398
399
      // Increase the offset for the next FE space by the number of DOFs owned 
      // by the rank in the current FE space.
400
      offset += (*data)[component].nRankDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
401
	
402
      // If there are no non local DOFs, continue with the next FE space.
403
      if (!globalMapping)
Thomas Witkowski's avatar
Thomas Witkowski committed
404
405
	continue;
      
406
      TEST_EXIT_DBG(dofComm != NULL)("No communicator given!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
407
      
408
409
410
      // === Communicate the matrix indices for all DOFs that are on some ===
      // === interior boundaries.                                         ===

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

	stdMpi.recv(rank);
      }
439

Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
      stdMpi.startCommunication();
      
442
      for (DofComm::Iterator it(dofComm->getRecvDofs(), 0, componentSpaces[component]); 
443
444
445
446
447
448
449
450
451
452
	   !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)
453
	      dofToMatIndex.add(component, dofMap[it.getDofIndex()].global, d);
454
	    else
455
	      dofToMatIndex.add(component, it.getDofIndex(), d);
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
458
	  }
	}
      }
459
460
461
462
463
464


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

      stdMpi.clear();

465
      for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, componentSpaces[component]); 
466
467
468
469
470
471
472
	   !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);
473
474
	      sendGlobalDofs.push_back(dofToMatIndex.get(component, 
							 dofMap[it.getDofIndex()].global));
475
	    } else {
476
	      sendGlobalDofs.push_back(dofToMatIndex.get(component, it.getDofIndex()));
477
478
479
480
481
482
483
484
485
486
487
488
489
	    }
	  }
	
	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);
	
	stdMpi.send(rank, sendGlobalDofs);
	stdMpi.recv(rank);
      }

      stdMpi.startCommunication();

490
      for (DofComm::Iterator it(dofComm->getPeriodicDofs(), 0, componentSpaces[component]); 
491
492
493
494
495
496
497
498
	   !it.end(); it.nextRank()) {

	int rank = it.getRank();
	if (meshLevel > 0)
	  rank = levelData->mapRank(rank, 0, meshLevel);

	int counter = 0;

499
	for (; !it.endDofIter(); it.nextDof()) {
500
501
	  if (dofMap.count(it.getDofIndex())) {
	    if (globalIndex) {
502
503
504
	      TEST_EXIT_DBG(counter + 2 <= stdMpi.getRecvData(it.getRank()).size())
		("Should not happen!\n");

505
	      dofToMatIndex.add(component, 
506
507
508
				stdMpi.getRecvData(it.getRank())[counter],
				stdMpi.getRecvData(it.getRank())[counter + 1]);
	      counter += 2;
509
	    } else {
510
511
	      dofToMatIndex.add(component, 
				it.getDofIndex(),
512
513
514
				stdMpi.getRecvData(it.getRank())[counter++]);
	    }
	  }
515
	}
516
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
517
518
    }
  }
519
520
521
522
523
524
525
526
527


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

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

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

Thomas Witkowski's avatar
bo eh    
Thomas Witkowski committed
550
551
552
553
554

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

555
    MSG("=== Parallel DOF mapping debug information ===\n");
Thomas Witkowski's avatar
bo eh    
Thomas Witkowski committed
556
    if (mode == COMPONENT_WISE) {
557
      MSG("  mapping is defined by component numbers!\n");
Thomas Witkowski's avatar
bo eh    
Thomas Witkowski committed
558
    } else {
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
      MSG("  mapping is defined by FE spaces!\n");
    }
    MSG("  matrix index is based on global DOF indices: %d\n", 
	needMatIndexFromGlobal);

    MSG("  nRankDofs = %d   nLocalDofs = %d   nOverallDofs = %d  rStartDofs = %d\n",
	nRankDofs, nLocalDofs, nOverallDofs, rStartDofs);

    int nComponents = componentSpaces.size();
    int nFeSpaces = feSpaces.size();
    MSG("  number of components: %d   number of different FE spaces: %d\n",
	nComponents, nFeSpaces);

    for (int i = 0; i < nComponents; i++) {
      MSG("  component %d:\n", i);
      MSG("    dof-to-mat-index has %d mappings\n", dofToMatIndex.getSize(i));
      if (dofToMatIndex.getSize(i) > 0) {
	MSG("    dof-to-mat-index starts with (%d -> %d) and ends with (%d -> %d)\n",
	    dofToMatIndex.getData(i).begin()->first,
	    dofToMatIndex.getData(i).begin()->second,
	    (dofToMatIndex.getData(i).end() - 1)->first,
	    (dofToMatIndex.getData(i).end() - 1)->second);
      }
Thomas Witkowski's avatar
bo eh    
Thomas Witkowski committed
582
583
    }
  }
584
}