ParallelDofMapping.cc 14.4 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
  void ParallelDofMapping::init(MeshLevelData &ldata,
231
				vector<const FiniteElemSpace*> &fe,
232
				vector<const FiniteElemSpace*> &uniqueFe,
233
				bool b)
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
234
  {
235
236
    FUNCNAME("ParallelDofMapping::init()");

237
    levelData = &ldata;
238
    isNonLocal = b;
239

240
    data->init(fe, uniqueFe, isNonLocal, ldata);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
241
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
242
243


244
245
246
247
  void ParallelDofMapping::clear()
  {
    FUNCNAME("ParallelDofMapping::clear()");

248
249
    TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n");

250
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
Thomas Witkowski's avatar
Thomas Witkowski committed
251
      it->clear();
252

253
254
255
256
257
    nRankDofs = -1;
    nLocalDofs = -1;
    nOverallDofs = -1;
    rStartDofs = -1;
    dofToMatIndex.clear();
258
259
  }

260
261
262
263
264
265
  
  void ParallelDofMapping::setMpiComm(MPI::Intracomm &m, int l)
  {
    mpiComm = m;
    meshLevel = l;

266
267
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->setMpiComm(m, l);
268
269
  }

270

271
  void ParallelDofMapping::setDofComm(DofComm &dc)
272
273
274
  {
    FUNCNAME("ParallelDofMapping::setDofComm()");

275
    dofComm = &dc;
276
277

    // Add the DOF communicator also to all FE space DOF mappings.
278
279
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->setDofComm(dc);
280
  }
281

282

283
  int ParallelDofMapping::computeRankDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
284
  {
285
    FUNCNAME("ParallelDofMapping::computeRankDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
286
287
    
    int result = 0;
288
289
290
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nRankDofs;
   
Thomas Witkowski's avatar
Thomas Witkowski committed
291
292
293
294
    return result;
  }


295
  int ParallelDofMapping::computeLocalDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
296
  {
297
    FUNCNAME("ParallelDofMapping::computeLocalDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
    
    int result = 0;
300
301
302
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nLocalDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
305
306
    return result;
  }


307
  int ParallelDofMapping::computeOverallDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
308
  {
309
    FUNCNAME("ParallelDofMapping::computeOverallDofs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
310
311

    int result = 0;
312
313
314
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->nOverallDofs;
       
Thomas Witkowski's avatar
Thomas Witkowski committed
315
316
317
318
    return result;
  }


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

    int result = 0;
324
325
326
    for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next())
      result += it->rStartDofs;

Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
329
330
331
332
333
334
    return result;
  }


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

335
    // First, update all FE space DOF mappings.
336
337
    for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next())
      it->update();
Thomas Witkowski's avatar
Thomas Witkowski committed
338
    
339
340
341
342
343
344
345
    // Compute all numbers from this mappings.
    nRankDofs = computeRankDofs();
    nLocalDofs = computeLocalDofs();
    nOverallDofs = computeOverallDofs();
    rStartDofs = computeStartDofs();
    
    // And finally, compute the matrix indices.
346
    computeMatIndex(needMatIndexFromGlobal);
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
  }
  
349
350
351
352
353

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

354
355
356
357
358
359
360
361
362
363
    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;

364
365
366
367

    update();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
368
  
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. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
384
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
385
386
387

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

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

412
      StdMpi<vector<DegreeOfFreedom> > stdMpi(mpiComm);
413
      for (DofComm::Iterator it(dofComm->getSendDofs(), 0, feSpaces[i]); 
Thomas Witkowski's avatar
Thomas Witkowski committed
414
415
416
417
418
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof())
	  if (dofMap.count(it.getDofIndex()))
419
	    if (globalIndex)
420
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, dofMap[it.getDofIndex()].global));
421
	    else
422
	      sendGlobalDofs.push_back(dofToMatIndex.get(i, 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, feSpaces[i]); 
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
443
444
445
446
447
448
449
450
451
452
453
454
455
      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
456
457
458
	  }
	}
      }
459
460
461
462
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


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

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


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

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

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

547
}