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

namespace AMDiS {

  using namespace std;

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

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

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

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

    Parameters::get("parallel->print timings", printTimings);
68
69
70
  }


71
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
72
  {
73
74
    FUNCNAME("PetscSolverFeti::initialize()");

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

78
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
79
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
80

Thomas Witkowski's avatar
Thomas Witkowski committed
81

Thomas Witkowski's avatar
Thomas Witkowski committed
82
83
84
85
86
87
88
    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");

Thomas Witkowski's avatar
Thomas Witkowski committed
89
90
      MSG("========= !!!! SUPER WARNING !!! ========\n");
      MSG("Use make use of stokes mode which is work in progress! We guarantee for nothing and even less than this!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
91
      pressureFeSpace = feSpaces[pressureComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
92
93
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
94
95
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
96

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

108
109
110
111
    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
112

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

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

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


  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
129
130
    double timeCounter = MPI::Wtime();

131
132
133
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
134

135
136
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

137
138
139
140
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
141
    if (fetiPreconditioner == FETI_DIRICHLET)
142
143
      interiorDofMap.clear();

144
145
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
146

147
148
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
149
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
150
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
151
    if (fetiPreconditioner == FETI_DIRICHLET)
152
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
153

154
155
156
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
157
158
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
159
    if (stokesMode) {
160
161
162
163
164
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

165
166
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
167
    
168
169
      createPrimals(feSpace);  

170
171
172
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
173

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
174
      createIndexB(feSpace);     
175
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
176
177
178
179

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
180
    if (fetiPreconditioner == FETI_DIRICHLET)
181
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
182

Thomas Witkowski's avatar
Thomas Witkowski committed
183
    if (stokesMode)
184
185
      interfaceDofMap.update();

186
187
188
189
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
190

191
192
193
    lagrangeMap.update();


194
195
196
197
198
199
200
201
202
203
204
205
    // === ===

    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();

206
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
207
208
209
210
211
212
213
214
215
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
221
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
	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);

	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
242
    }
243

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
247
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
248
249

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
250
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
251
      timeCounter = MPI::Wtime() - timeCounter;
252
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
253
254
	  timeCounter);
    }
255
256
257
  }


258
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
259
  {
260
    FUNCNAME("PetscSolverFeti::createPrimals()");  
261

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    if (feSpace == pressureFeSpace)
263
264
      return;

265
266
267
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    // Set of DOF indices that are considered to be primal variables.
269
    DofContainerSet& vertices = 
270
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
271
272

    DofIndexSet primals;
273
    for (DofContainerSet::iterator it = vertices.begin(); 
274
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
275
      if (meshLevel == 0) {
276
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
277
      } else {
278
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
279
280
281
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
282
283
284
285
286
287
288
	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);
	}
      }
289
    }
290

291
292
293
294

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

295
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
296
297
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
298
      else
299
  	primalDofMap[feSpace].insertNonRankDof(*it);
300
301
302
  }


303
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
304
305
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
306

Thomas Witkowski's avatar
Thomas Witkowski committed
307
    if (feSpace == pressureFeSpace)
308
309
      return;

310
311
312
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
313
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
314
315
316
317

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
318
319
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
320
321
322
323
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
324
325
326
  }

  
327
328
329
330
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
331
    if (feSpace != pressureFeSpace)
332
333
334
335
336
337
338
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
339
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
340
	interfaceDofMap[feSpace].insertRankDof(**it);
341
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
342
	interfaceDofMap[feSpace].insertNonRankDof(**it);
343
344
345
  }


346
347
348
349
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
350
    if (feSpace == pressureFeSpace)
351
352
      return;

353
354
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
355
356
357
    // 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
358
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
359

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
360
    if (meshLevel > 0) {
361
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
362
363
364
365
366
367
368
369
370

      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
371
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
372
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
373
	  else
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
	    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
389
390
391
392
393
394
	   !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());
    }
395

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
396
397
398
399
400
401
402
    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)).         ===

403
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
404
405
406
407
408
409
410
    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);

411
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
412
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
413
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
414
415
416
417
	}
      }
    }

418
419
420
421

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

422
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
423

424
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
425
426
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
427
	if (!isPrimal(feSpace, it.getDofIndex()))
428
429
430
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
431
432
433

    stdMpi.updateSendDataSize();

434
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
435
	 !it.end(); it.nextRank()) {
436
      bool recvFromRank = false;
437
      for (; !it.endDofIter(); it.nextDof()) {
438
	if (!isPrimal(feSpace, it.getDofIndex())) {
439
440
441
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
442
443
444
	    recvFromRank = true;
	    break;
	  }
445
	}
446
      }
447
448

      if (recvFromRank)
449
	stdMpi.recv(it.getRank());
450
    }
451

452
453
    stdMpi.startCommunication();

454
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
455
	 !it.end(); it.nextRank()) {
456
      int i = 0;
457
      for (; !it.endDofIter(); it.nextDof())
458
	if (!isPrimal(feSpace, it.getDofIndex()))
459
460
461
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
462
463
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
464
465
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
466
467
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
468

469
470
471
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

472
    int nRankLagrange = 0;
473
    DofMap& dualMap = dualDofMap[feSpace].getMap();
474
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
475
476
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
477
	int degree = boundaryDofRanks[feSpace][it->first].size();
478
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
479
      } else {
480
	lagrangeMap[feSpace].insertNonRankDof(it->first);
481
482
      }
    }
483
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
484
485
486
  }


487
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
488
  {
489
    FUNCNAME("PetscSolverFeti::createIndexB()");
490

491
    DOFAdmin* admin = feSpace->getAdmin();
492
493
494
495

    // === 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.                                 ===
496

497
    int nLocalInterior = 0;    
498
    for (int i = 0; i < admin->getUsedSize(); i++) {
499
      if (admin->isDofFree(i) == false && 
500
	  isPrimal(feSpace, i) == false &&
501
502
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
503
504
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
505

506
	  if (fetiPreconditioner == FETI_DIRICHLET)
507
508
509
510
511
512
513
514
	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	  
	  nLocalInterior++;	  
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	    localDofMap[feSpace].insertRankDof(i);
	  else
	    localDofMap[feSpace].insertNonRankDof(i);
515
516
517

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
518
	}
519
      }
520
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
521
    
522
523
    // === And finally, add the global indicies of all dual nodes. ===

524
525
526
527
528
529
530
531
532
533
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
      if (meshLevel == 0)
	localDofMap[feSpace].insertRankDof(it->first);
      else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
534
535
536
  }


537
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
538
539
540
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
541
    double wtime = MPI::Wtime();
542
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
543

544
545
    // === Create distributed matrix for Lagrange constraints. ===

546
547
548
549
550
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
551
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
552

553
554
555
556
557
558
559
    // === 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
560
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
561
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
562

563
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
564
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
565
566
567
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
568
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
569

Thomas Witkowski's avatar
Thomas Witkowski committed
570
	// Copy set of all ranks that contain this dual node.
571
572
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
573
574
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
575
576

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
577
578
579
580
	
	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
581
582
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
583

584
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
585
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
586
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
587
	    index++;	      
588
589
590
591
592
593
594
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
596
597
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
598
599
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
600
    }
601
602
603
  }


604
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
605
  {
606
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
607

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

Thomas Witkowski's avatar
Thomas Witkowski committed
611
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
612
      
613
614
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
615

616
      MatCreateShell(mpiCommGlobal,
617
618
619
620
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
621
622
623
624
625
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
626
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
627
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
628
629
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
630
631
      KSPSetFromOptions(ksp_schur_primal);
    } else {
632
633
      MSG("Create direct schur primal solver!\n");

634
635
      double wtime = MPI::Wtime();

636
637
638
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

639
640
641
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
642

Thomas Witkowski's avatar
Thomas Witkowski committed
643
      Mat matBPi;
644
645
646
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
647
648
649
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

650

651
652
653
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
654

655
656
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

657
      for (int i = 0; i < nRowsRankB; i++) {
658
	PetscInt row = localDofMap.getStartDofs() + i;
659
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
660
661
662
663
664

	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]));

665
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
666
667
      }

668
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
669
		primalDofMap.getLocalDofs())
670
	("Should not happen!\n");
671

672
673
674
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
675
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
676
677
678
679
680
681
682

 	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);
683
	
Thomas Witkowski's avatar
Thomas Witkowski committed
684
	subdomain->solve(tmpVec, tmpVec);
685

Thomas Witkowski's avatar
Thomas Witkowski committed
686
	PetscScalar *tmpValues;
687
	VecGetArray(tmpVec, &tmpValues);
688
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
689
	  MatSetValue(matBPi, 
690
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
691
692
693
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
694
695
696
697
698
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
701

702
703
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
704
705

      Mat matPrimal;
706
707
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
708
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
709
710

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
711
      MatDestroy(&matBPi);
712

713
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
714
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
715
		      SAME_NONZERO_PATTERN);
716
717
718
719
720
721
      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);
722
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
723

Thomas Witkowski's avatar
Thomas Witkowski committed
724
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
725
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
726
727
728
729
730
731
732
733
734
735
	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
736
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
737
738
739
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
740
    }
741
742
743
744
745
746
747
  }


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

748
    if (schurPrimalSolver == 0) {
749
      schurPrimalData.subSolver = NULL;
750

751
752
753
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
754
755
756
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
757
758
759
  }


760
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
761
762
763
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

764
765
    // === Create FETI-DP solver object. ===

766
767
768
769
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
770
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
771
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
772
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
773

Thomas Witkowski's avatar
Thomas Witkowski committed
774
    if (stokesMode == false) {      
775
776
777
778
779
780
781
782
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
783
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
784
785
786
787
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

788
789
790
791
792
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
793
		     &fetiData, &mat_feti);
794
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
795
			   (void(*)(void))petscMultMatFetiInterface);      
796
    }
797

798
    KSPCreate(mpiCommGlobal, &ksp_feti);
799
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
800
801
802
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
803
    KSPSetFromOptions(ksp_feti);
804
805


806
    // === Create FETI-DP preconditioner object. ===
807

Thomas Witkowski's avatar
Thomas Witkowski committed
808
    if (fetiPreconditioner != FETI_NONE || stokesMode) {
809
810
811
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
812

813

Thomas Witkowski's avatar
Thomas Witkowski committed
814
815
816
817
818
    // === Create null space objects. ===
    
    createNullSpace();


819
    switch (fetiPreconditioner) {
820
821
    case FETI_DIRICHLET:
      TEST_EXIT(meshLevel == 0)
822
823
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

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

827
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
828
829
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
830
831
832
833
834
835
      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);
836
837
838
839
840
841
842
843
844
      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;
      
845
846
      localDofMap.createVec(fetiDirichletPreconData.tmp_vec_b, 
			    nGlobalOverallInterior);
847
848
849
850
851
852
      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));
853

854
855
856
      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

857
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
858
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
859
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
860
	  DegreeOfFreedom d = it->first;	  
861
862
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
863
864
865
866
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

867
868
869
870
      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
871
872
873
874
875
876
877
878


      // 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
879
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
880
881
882
	MSG("FETI-DP timing 08: %.5f seconds (factorization of Dirichlet preconditoner matrices)\n",
	    MPI::Wtime() - wtime);
      }
883
884
885
886
      
      break;

    case FETI_LUMPED:
Thomas Witkowski's avatar
Thomas Witkowski committed
887
888
889
      {
	FetiLumpedPreconData *lumpedData = 
	  (stokesMode ? &fetiInterfaceLumpedPreconData : &fetiLumpedPreconData);
890

Thomas Witkowski's avatar
Thomas Witkowski committed
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
	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;
	  
	  DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
	  for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
	    DegreeOfFreedom d = it->first;
	    int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	    int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	    lumpedData->localToDualMap[matIndexLocal] = matIndexDual;
	  }
	}
	
	if (stokesMode) {
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
	  DOFVector<double> tmpvec(pressureFeSpace, "tmpvec");
	  createH2vec(tmpvec);
	  interfaceDofMap.createVec(fetiInterfaceLumpedPreconData.h2vec);

	  DofMap& m = interfaceDofMap[pressureFeSpace].getMap();
	  for (DofMap::iterator it = m.begin(); it != m.end(); ++it) {
	    if (meshDistributor->getDofMap()[pressureFeSpace].isRankDof(it->first)) {
	      int index = interfaceDofMap.getMatIndex(pressureComponent, it->first);
	      VecSetValue(fetiInterfaceLumpedPreconData.h2vec, 
			  index, tmpvec[it->first], INSERT_VALUES);
	    }
	  }  

	  VecAssemblyBegin(fetiInterfaceLumpedPreconData.h2vec);
	  VecAssemblyEnd(fetiInterfaceLumpedPreconData.h2vec);
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
930
931
932
933
934
935
936
937
938
939
940
941
	  localDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_vec_b1);
	  fetiInterfaceLumpedPreconData.subSolver = subdomain;	
	}
	
	KSPGetPC(ksp_feti, &precon_feti);
	PCSetType(precon_feti, PCSHELL);
	if (stokesMode) {
	  PCShellSetContext(precon_feti, static_cast<void*>(&fetiInterfaceLumpedPreconData));
	  PCShellSetApply(precon_feti, petscApplyFetiInterfaceLumpedPrecon);
	} else {
	  PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
	  PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
Thomas Witkowski's avatar
Thomas Witkowski committed
942
943
944
	}
      }

945
946
      break;
    default:
947
948
      break;
    }
949
950
951
952
953
954
955
  }
  

  void PetscSolverFeti::destroyFetiKsp()
  {
    FUNCNAME("PetscSolverFeti::destroyFetiKsp()");

956
957
    // === Destroy FETI-DP solver object. ===

958
    fetiData