PetscSolverFeti.cc 63.7 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 "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/PetscSolverFeti.h"
16
17
#include "parallel/PetscSolverFetiDebug.h"
#include "parallel/PetscSolverFetiMonitor.h"
18
#include "parallel/PetscSolverFetiStructs.h"
19
20
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
21
22
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
23
#include "parallel/PetscSolverGlobalMatrix.h"
24
#include "io/VtkWriter.h"
25
26
27
28
29

namespace AMDiS {

  using namespace std;

30
31
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
32
      schurPrimalSolver(0),
33
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
34
      subdomain(NULL),
35
      massMatrixSolver(NULL),
36
37
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
38
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
39
      printTimings(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
40
41
      stokesMode(false),
      pressureComponent(-1)
42
43
44
45
46
47
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
48
      MSG("Create FETI-DP solver with no preconditioner!\n");
49
50
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
51
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
52
53
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
54
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
55
56
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
57
58
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
59
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
60
61
62
63
64

    Parameters::get("parallel->feti->schur primal solver", schurPrimalSolver);
    TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
      ("Wrong solver \"%d\"for the Schur primal complement!\n", 
       schurPrimalSolver);
65
66

    Parameters::get("parallel->multi level test", multiLevelTest);
67
68
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
69
70

    Parameters::get("parallel->print timings", printTimings);
71
72
73
  }


74
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
75
  {
76
77
    FUNCNAME("PetscSolverFeti::initialize()");

78
79
80
    TEST_EXIT_DBG(meshLevel + 1 == meshDistributor->getMeshLevelData().getLevelNumber())
      ("Mesh hierarchy does not contain %d levels!\n", meshLevel + 1);

81
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
82
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
83

Thomas Witkowski's avatar
Thomas Witkowski committed
84

Thomas Witkowski's avatar
Thomas Witkowski committed
85
86
87
88
89
90
91
92
    stokesMode = false;
    Parameters::get("parallel->feti->stokes mode", stokesMode);
    if (stokesMode) {
      Parameters::get("parallel->feti->pressure component", pressureComponent);
      TEST_EXIT(pressureComponent >= 0)
	("FETI-DP in Stokes mode, no pressure component defined!\n");

      pressureFeSpace = feSpaces[pressureComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
93
94
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
95
96
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
97

98
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
99
	subdomain->setMeshDistributor(meshDistributor, 
100
				      mpiCommGlobal, mpiCommLocal);
101
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
102
	subdomain->setMeshDistributor(meshDistributor, 
103
104
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
105
	subdomain->setLevel(meshLevel);
106
107
      }
    }
108

109
110
111
112
    primalDofMap.init(levelData, feSpaces, uniqueFe);   
    dualDofMap.init(levelData, feSpaces, uniqueFe, false);
    localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0);
    lagrangeMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
113

Thomas Witkowski's avatar
Thomas Witkowski committed
114
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
115
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
116

117
    if (fetiPreconditioner == FETI_DIRICHLET) {
118
119
120
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

121
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
122
    }
123
124
125
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
126
127
128
129
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
130
131
132
133
134
    bool removeDirichletRows = false;
    Parameters::get("parallel->feti->remove dirichlet", removeDirichletRows);
    if (!removeDirichletRows)
      return;

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
    int nComponents = mat.getSize();
    for (int i = 0; i < nComponents; i++) {
      DOFMatrix* dofMat = mat[i][i];
      if (!dofMat)
	continue;
      
      const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
      std::set<DegreeOfFreedom>& dRows = dofMat->getDirichletRows();
      if (dirichletRows.count(feSpace)) {
	WARNING("Implement test that for all components dirichlet rows are equal!\n");
      } else {
	dirichletRows[feSpace] = dRows;
      }
    }
  }


152
153
154
155
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
156
157
    double timeCounter = MPI::Wtime();

158
159
160
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
161

162
163
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

164
165
166
167
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
168
    if (fetiPreconditioner == FETI_DIRICHLET)
169
170
      interiorDofMap.clear();

171
172
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
173

174
175
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
176
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
177
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
178
    if (fetiPreconditioner == FETI_DIRICHLET)
179
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
180

181
182
183
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
184
185
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
186
    if (stokesMode) {
187
188
189
190
191
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

192
193
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
194
    
195
196
      createPrimals(feSpace);  

197
198
199
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
200

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
201
      createIndexB(feSpace);     
202
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
203
204
205
206

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
207
    if (fetiPreconditioner == FETI_DIRICHLET)
208
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
209

Thomas Witkowski's avatar
Thomas Witkowski committed
210
    if (stokesMode)
211
212
      interfaceDofMap.update();

213
214
215
216
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
217

218
219
220
    lagrangeMap.update();


221
222
223
224
225
226
227
228
229
230
231
232
    // === ===

    if (meshLevel == 0) {
      rStartInterior = 0;
      nGlobalOverallInterior = localDofMap.getOverallDofs();
    } else {
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

      int groupRowsInterior = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	groupRowsInterior = localDofMap.getOverallDofs();

233
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
234
235
236
237
238
239
240
241
242
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	tmp = rStartInterior;

      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }

243
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
Thomas Witkowski's avatar
Thomas Witkowski committed
244
245
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
246
      MSG("FETI-DP data for %d-ith FE space: %p\n", i, feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
247

Thomas Witkowski's avatar
Thomas Witkowski committed
248
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
	MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
	    primalDofMap[feSpace].nRankDofs, 
	    primalDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
	    dualDofMap[feSpace].nRankDofs, 
	    dualDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
	    lagrangeMap[feSpace].nRankDofs, 
	    lagrangeMap[feSpace].nOverallDofs);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
265
266
267
// 	TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
// 		      static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
// 	  ("Should not happen!\n");	
Thomas Witkowski's avatar
Thomas Witkowski committed
268
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    }
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
274
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
275
276

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
277
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
278
      timeCounter = MPI::Wtime() - timeCounter;
279
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
	  timeCounter);
    }
282
283
284
  }


285
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
286
  {
287
    FUNCNAME("PetscSolverFeti::createPrimals()");  
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
    if (feSpace == pressureFeSpace)
290
291
      return;

292
293
294
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
295
    // Set of DOF indices that are considered to be primal variables.
296
    DofContainerSet& vertices = 
297
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
298
299

    DofIndexSet primals;
300
    for (DofContainerSet::iterator it = vertices.begin(); 
301
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
302
303
      if (dirichletRows[feSpace].count(**it))
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
304

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
305
306
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
307
      } else {
308
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
309
310
311
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
312
313
314
315
316
317
318
	if (fabs(c[0]) < e || fabs(c[1]) < e ||
	    fabs(c[0] - 25.0) < e || fabs(c[1] - 25.0) < e ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 12.5) < e)) {
	  MSG("PRIMAL COORD %f %f\n", c[0], c[1]);
	  primals.insert(**it);
	}
      }
319
    }
320

321
322
323
324

    // === Calculate the number of primals that are owned by the rank and ===
    // === create local indices of the primals starting at zero.          ===

325
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
326
327
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
328
      else
329
  	primalDofMap[feSpace].insertNonRankDof(*it);
330
331
332
  }


333
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
334
335
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
336

Thomas Witkowski's avatar
Thomas Witkowski committed
337
    if (feSpace == pressureFeSpace)
338
339
      return;

340
341
342
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
343
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
344
    
345
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
346
347
348
349
350
351
352
353
354
355
356
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;

      if (isPrimal(feSpace, **it))
	continue;

      if (meshLevel == 0) {
	dualDofMap[feSpace].insertRankDof(**it);
      } else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
357
	  dualDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
358
359
      }	  
    }
360
361
362
  }

  
363
364
365
366
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
367
    if (feSpace != pressureFeSpace)
368
369
370
371
372
373
      return;

    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
374
375
376
377
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;      
      
378
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
379
	interfaceDofMap[feSpace].insertRankDof(**it);
380
      else
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
381
382
	interfaceDofMap[feSpace].insertNonRankDof(**it);      
    }
383
384
385
  }


386
387
388
389
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    if (feSpace == pressureFeSpace)
391
392
      return;

393
394
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
397
    // Stores for all rank owned communication DOFs, if the counterpart is
    // a rank owned DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFs that fulfill this requirenment.
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
398
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
399

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
400
    if (meshLevel > 0) {
401
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
402
403
404
405
406
407
408
409
410

      for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank()) {

	vector<int> subdomainRankDofs;
	subdomainRankDofs.reserve(it.getDofs().size());

	for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
411
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
412
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
413
	  else
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
	    subdomainRankDofs.push_back(0);
	}

	stdMpi.send(it.getRank(), subdomainRankDofs);
      }	     

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace); 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
429
430
431
432
433
434
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex()))
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
435

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
436
437
438
439
440
441
442
    if (dualDofMap[feSpace].nLocalDofs == 0)
      return;


    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===

443
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
444
445
446
447
448
449
450
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(feSpace, it.getDofIndex())) {
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

451
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
452
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
453
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
454
455
456
457
	}
      }
    }

458
459
460
461

    // === Communicate these sets for all rank owned dual nodes to other ===
    // === ranks that also have this node.                               ===

462
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
463

464
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
465
466
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
467
	if (!isPrimal(feSpace, it.getDofIndex()))
468
469
470
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
471
472
473

    stdMpi.updateSendDataSize();

474
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
475
	 !it.end(); it.nextRank()) {
476
      bool recvFromRank = false;
477
      for (; !it.endDofIter(); it.nextDof()) {
478
	if (!isPrimal(feSpace, it.getDofIndex())) {
479
480
481
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
482
483
484
	    recvFromRank = true;
	    break;
	  }
485
	}
486
      }
487
488

      if (recvFromRank)
489
	stdMpi.recv(it.getRank());
490
    }
491

492
493
    stdMpi.startCommunication();

494
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
495
	 !it.end(); it.nextRank()) {
496
      int i = 0;
497
      for (; !it.endDofIter(); it.nextDof())
498
	if (!isPrimal(feSpace, it.getDofIndex()))
499
500
501
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
502
503
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
504
505
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
506
507
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
508

509
510
511
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

512
    int nRankLagrange = 0;
513
    DofMap& dualMap = dualDofMap[feSpace].getMap();
514
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
515
516
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
517
	int degree = boundaryDofRanks[feSpace][it->first].size();
518
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
519
      } else {
520
	lagrangeMap[feSpace].insertNonRankDof(it->first);
521
522
      }
    }
523
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
524
525
526
  }


527
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
528
  {
529
    FUNCNAME("PetscSolverFeti::createIndexB()");
530

531
    DOFAdmin* admin = feSpace->getAdmin();
532
533
534
535

    // === To ensure that all interior node on each rank are listen first in ===
    // === the global index of all B nodes, insert all interior nodes first, ===
    // === without defining a correct index.                                 ===
536

537
    int nLocalInterior = 0;    
538
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
539
540
541
542
543
544
      if (admin->isDofFree(i) ||
	  isPrimal(feSpace, i) ||
	  isDual(feSpace, i) ||
	  isInterface(feSpace, i) ||
	  dirichletRows[feSpace].count(i))
	continue;      
545

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
      if (meshLevel == 0) {
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);
	
	if (fetiPreconditioner == FETI_DIRICHLET)
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	
	nLocalInterior++;	
      } else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	  localDofMap[feSpace].insertRankDof(i);
	else
	  localDofMap[feSpace].insertNonRankDof(i);
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
561
      }
562
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
563
    
564
565
    // === And finally, add the global indicies of all dual nodes. ===

566
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
567
568
	 it != dualDofMap[feSpace].getMap().end(); ++it) {
      if (meshLevel == 0) {
569
	localDofMap[feSpace].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
570
      } else {
571
572
573
574
575
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
576
    }
577
578
579
  }


580
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
581
582
583
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
584
    double wtime = MPI::Wtime();
585
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
586

587
588
    // === Create distributed matrix for Lagrange constraints. ===

589
590
591
592
593
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
594
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
595

596
597
598
599
600
601
602
    // === Create for all duals the corresponding Lagrange constraints. On ===
    // === each rank we traverse all pairs (n, m) of ranks, with n < m,    ===
    // === that contain this node. If the current rank number is r, and    ===
    // === n == r, the rank sets 1.0 for the corresponding constraint, if  ===
    // === m == r, than the rank sets -1.0 for the corresponding           ===
    // === constraint.                                                     ===

Thomas Witkowski's avatar
Thomas Witkowski committed
603
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
604
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
605

606
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
607
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
608
609
610
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
611
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
612

Thomas Witkowski's avatar
Thomas Witkowski committed
613
	// Copy set of all ranks that contain this dual node.
614
615
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
616
617
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
618
619

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
620
621
622
623
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
624
625
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
626

627
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
628
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
629
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
630
	    index++;	      
631
632
633
634
635
636
637
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
638

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
639
640
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
641
642
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
643
    }
644
645
646
  }


647
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
648
  {
649
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
650

Thomas Witkowski's avatar
Thomas Witkowski committed
651
    if (schurPrimalSolver == 0) {
652
653
      MSG("Create iterative schur primal solver!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
654
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
655
      
656
657
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
658

659
      MatCreateShell(mpiCommGlobal,
660
661
662
663
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
666
667
668
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
669
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
670
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
671
672
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
673
674
      KSPSetFromOptions(ksp_schur_primal);
    } else {
675
676
      MSG("Create direct schur primal solver!\n");

677
678
      double wtime = MPI::Wtime();

679
680
681
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

682
683
684
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
685

Thomas Witkowski's avatar
Thomas Witkowski committed
686
      Mat matBPi;
687
688
689
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
690
691
692
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

693

694
695
696
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
697

698
699
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

700
      for (int i = 0; i < nRowsRankB; i++) {
701
	PetscInt row = localDofMap.getStartDofs() + i;
702
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
703
704
705
706
707

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

708
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
709
710
      }

711
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
712
		primalDofMap.getLocalDofs())
713
	("Should not happen!\n");
714

715
716
717
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
718
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
719
720
721
722
723
724
725

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
726

Thomas Witkowski's avatar
Thomas Witkowski committed
727
	subdomain->solve(tmpVec, tmpVec);
728

Thomas Witkowski's avatar
Thomas Witkowski committed
729
	PetscScalar *tmpValues;
730
	VecGetArray(tmpVec, &tmpValues);
731
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
732
	  MatSetValue(matBPi, 
733
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
734
735
736
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
737
738
739
740
741
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
742
743
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
744

745
746
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
747
748

      Mat matPrimal;
749
750
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
751
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
752
753

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
754
      MatDestroy(&matBPi);
755

756
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
757
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
758
		      SAME_NONZERO_PATTERN);
759
760
761
762
763
764
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPPREONLY);
      PC pc_schur_primal;      
      KSPGetPC(ksp_schur_primal, &pc_schur_primal);
      PCSetType(pc_schur_primal, PCLU);
      PCFactorSetMatSolverPackage(pc_schur_primal, MATSOLVERMUMPS);
765
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
766

Thomas Witkowski's avatar
Thomas Witkowski committed
767
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
768
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
769
770
771
772
773
774
775
776
777
778
	MatInfo minfo;
	MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
	MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);
	
	MSG("FETI-DP timing 06: %.5f seconds (creation of schur primal matrix)\n",
	    MPI::Wtime() - wtime);

	wtime = MPI::Wtime();
	KSPSetUp(ksp_schur_primal);
	KSPSetUpOnBlocks(ksp_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
779
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
780
781
782
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
783
    }
784
785
786
787
788
789
790
  }


  void PetscSolverFeti::destroySchurPrimalKsp()
  {
    FUNCNAME("PetscSolverFeti::destroySchurPrimal()");

791
    if (schurPrimalSolver == 0) {
792
      schurPrimalData.subSolver = NULL;
793

794
795
796
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
797
798
799
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
800
801
802
  }


803
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
804
805
806
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

807
808
    // === Create FETI-DP solver object. ===

809
810
811
812
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
813
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
814
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
815
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
816

Thomas Witkowski's avatar
Thomas Witkowski committed
817
    if (stokesMode == false) {      
818
819
820
821
822
823
824
825
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
826
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
827
828
829
830
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

831
832
833
834
835
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
836
		     &fetiData, &mat_feti);
837
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
838
			   (void(*)(void))petscMultMatFetiInterface);      
839
    }
840

841
    KSPCreate(mpiCommGlobal, &ksp_feti);
842
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
843
844
845
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
846
    KSPSetFromOptions(ksp_feti);
847
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
848
      KSPMonitorSet(ksp_feti, KSPMonitorFetiStokes, &fetiKspData, PETSC_NULL);
849

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
850
851


852
    // === Create FETI-DP preconditioner object. ===
853

Thomas Witkowski's avatar
Thomas Witkowski committed
854
    if (fetiPreconditioner != FETI_NONE || stokesMode) {
855
856
857
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
858

859

Thomas Witkowski's avatar
Thomas Witkowski committed
860
861
862
863
864
    // === Create null space objects. ===
    
    createNullSpace();


865
    switch (fetiPreconditioner) {
866
867
    case FETI_DIRICHLET:
      TEST_EXIT(meshLevel == 0)
868
869
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
870
      TEST_EXIT(!stokesMode)
871
872
	("Stokes mode does not yet support the Dirichlet precondition!\n");

873
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
874
875
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
876
877
878
879
880
881
      KSPSetOptionsPrefix(ksp_interior, "precon_interior_");
      KSPSetType(ksp_interior, KSPPREONLY);
      PC pc_interior;
      KSPGetPC(ksp_interior, &pc_interior);
      PCSetType(pc_interior, PCLU);
      PCFactorSetMatSolverPackage(pc_interior, MATSOLVERUMFPACK);
882
883
884
885
886
887
888
889
890
      KSPSetFromOptions(ksp_interior);
            
      fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior;
      fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals;
      fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals;
      fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
      fetiDirichletPreconData.ksp_interior = &ksp_interior;
      
891
892
      localDofMap.createVec(fetiDirichletPreconData.tmp_vec_b, 
			    nGlobalOverallInterior);
893
894
895
896
897
898
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals1));
      MatGetVecs(mat_interior_interior, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_interior));
899

900
901
902
      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

903
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
904
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
905
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
906
	  DegreeOfFreedom d = it->first;	  
907
908
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
909
910
911
912
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

913
914
915
916
      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
Thomas Witkowski's avatar
Thomas Witkowski committed
917
918
919
920
921
922
923
924


      // For the case, that we want to print the timings, we force the LU 
      // factorization of the local problems to be done here explicitly.
      if (printTimings) {
	double wtime = MPI::Wtime();
	KSPSetUp(ksp_interior);
	KSPSetUpOnBlocks(ksp_interior);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
925
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
926
927
928
	MSG("FETI-DP timing 08: %.5f seconds (factorization of Dirichlet preconditoner matrices)\n",
	    MPI::Wtime() - wtime);
      }
929
930
931
932
      
      break;

    case FETI_LUMPED:
Thomas Witkowski's avatar
Thomas Witkowski committed
933
934
935
      {
	FetiLumpedPreconData *lumpedData = 
	  (stokesMode ? &fetiInterfaceLumpedPreconData : &fetiLumpedPreconData);
936

Thomas Witkowski's avatar
Thomas Witkowski committed
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
	lumpedData->mat_lagrange_scaled = &mat_lagrange_scaled;
	lumpedData->mat_duals_duals = &mat_duals_duals;
	
	localDofMap.createVec(lumpedData->tmp_vec_b0);
	MatGetVecs(mat_duals_duals, PETSC_NULL, 
		   &(lumpedData->tmp_vec_duals0));
	MatGetVecs(mat_duals_duals, PETSC_NULL, 
		   &(lumpedData->tmp_vec_duals1));
	
	for (unsigned int i = 0; i < feSpaces.size(); i++) {
	  if (stokesMode && i == pressureComponent)
	    continue;</