PetscSolverFeti.cc 69.6 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) {
Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
278
279
280
281
282
283
	MSG("WARNING: Modified primal detection algorithm!\n");

	double e = 1e-3;
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

	if (fabs(c[0]) > e && fabs(c[1]) > e && fabs(c[0] - 1.0) > e && fabs(c[1] - 1.0) > e)
	  primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
284
      } else {
285
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
286
287
288
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
289
290
291
292
293
294
295
	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);
	}
      }
296
    }
297

298
299
300
301

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

302
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
303
304
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
305
      else
306
  	primalDofMap[feSpace].insertNonRankDof(*it);
307
308
309
  }


310
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
311
312
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
313

Thomas Witkowski's avatar
Thomas Witkowski committed
314
    if (feSpace == pressureFeSpace)
315
316
      return;

317
318
319
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
320
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
321
322
323
324

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
325
326
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
327
328
329
330
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
331
332
333
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
338
    if (feSpace != pressureFeSpace)
339
340
341
342
343
344
345
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
346
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
347
	interfaceDofMap[feSpace].insertRankDof(**it);
348
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
349
	interfaceDofMap[feSpace].insertNonRankDof(**it);
350
351
352
  }


353
354
355
356
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
357
    if (feSpace == pressureFeSpace)
358
359
      return;

360
361
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
364
    // 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
365
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
366

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
367
    if (meshLevel > 0) {
368
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
369
370
371
372
373
374
375
376
377

      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
378
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
379
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
380
	  else
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
	    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
396
397
398
399
400
401
	   !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());
    }
402

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
403
404
405
406
407
408
409
    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)).         ===

410
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
411
412
413
414
415
416
417
    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);

418
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
419
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
420
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
421
422
423
424
	}
      }
    }

425
426
427
428

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

429
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
430

431
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
432
433
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
434
	if (!isPrimal(feSpace, it.getDofIndex()))
435
436
437
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
438
439
440

    stdMpi.updateSendDataSize();

441
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
442
	 !it.end(); it.nextRank()) {
443
      bool recvFromRank = false;
444
      for (; !it.endDofIter(); it.nextDof()) {
445
	if (!isPrimal(feSpace, it.getDofIndex())) {
446
447
448
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
449
450
451
	    recvFromRank = true;
	    break;
	  }
452
	}
453
      }
454
455

      if (recvFromRank)
456
	stdMpi.recv(it.getRank());
457
    }
458

459
460
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
475

476
477
478
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


494
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
495
  {
496
    FUNCNAME("PetscSolverFeti::createIndexB()");
497

498
    DOFAdmin* admin = feSpace->getAdmin();
499
500
501
502

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

504
    int nLocalInterior = 0;    
505
    for (int i = 0; i < admin->getUsedSize(); i++) {
506
      if (admin->isDofFree(i) == false && 
507
	  isPrimal(feSpace, i) == false &&
508
509
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
510
511
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
512

513
	  if (fetiPreconditioner == FETI_DIRICHLET)
514
515
516
517
518
519
520
521
	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	  
	  nLocalInterior++;	  
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	    localDofMap[feSpace].insertRankDof(i);
	  else
	    localDofMap[feSpace].insertNonRankDof(i);
522
523
524

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
525
	}
526
      }
527
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
528
    
529
530
    // === And finally, add the global indicies of all dual nodes. ===

531
532
533
534
535
536
537
538
539
540
    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);
      }
541
542
543
  }


544
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
545
546
547
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
548
    double wtime = MPI::Wtime();
549
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
550

551
552
    // === Create distributed matrix for Lagrange constraints. ===

553
554
555
556
557
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
558
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
559

560
561
562
563
564
565
566
    // === 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
567
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
568
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
569

570
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
571
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
572
573
574
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
575
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
576

Thomas Witkowski's avatar
Thomas Witkowski committed
577
	// Copy set of all ranks that contain this dual node.
578
579
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
580
581
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
582
583

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
584
585
586
587
	
	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
588
589
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
590

591
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
592
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
593
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
594
	    index++;	      
595
596
597
598
599
600
601
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
603
604
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
605
606
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
607
    }
608
609
610
  }


611
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
612
  {
613
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
614

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

Thomas Witkowski's avatar
Thomas Witkowski committed
618
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
619
      
620
621
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
622

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

641
642
      double wtime = MPI::Wtime();

643
644
645
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

646
647
648
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
649

Thomas Witkowski's avatar
Thomas Witkowski committed
650
      Mat matBPi;
651
652
653
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
654
655
656
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

657

658
659
660
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
661

662
663
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

664
      for (int i = 0; i < nRowsRankB; i++) {
665
	PetscInt row = localDofMap.getStartDofs() + i;
666
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
667
668
669
670
671

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

672
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
673
674
      }

675
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
676
		primalDofMap.getLocalDofs())
677
	("Should not happen!\n");
678

679
680
681
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
682
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
683
684
685
686
687
688
689

 	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);
690
	
Thomas Witkowski's avatar
Thomas Witkowski committed
691
	subdomain->solve(tmpVec, tmpVec);
692

Thomas Witkowski's avatar
Thomas Witkowski committed
693
	PetscScalar *tmpValues;
694
	VecGetArray(tmpVec, &tmpValues);
695
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
696
	  MatSetValue(matBPi, 
697
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
698
699
700
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
701
702
703
704
705
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
706
707
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
708

709
710
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
711
712

      Mat matPrimal;
713
714
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
715
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
716
717

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
718
      MatDestroy(&matBPi);
719

720
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
721
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
722
		      SAME_NONZERO_PATTERN);
723
724
725
726
727
728
      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);
729
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
730

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


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

755
    if (schurPrimalSolver == 0) {
756
      schurPrimalData.subSolver = NULL;
757

758
759
760
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
761
762
763
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
764
765
766
  }


767
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
768
769
770
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

771
772
    // === Create FETI-DP solver object. ===

773
774
775
776
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
777
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
778
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
779
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
780

Thomas Witkowski's avatar
Thomas Witkowski committed
781
    if (stokesMode == false) {      
782
783
784
785
786
787
788
789
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
790
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
791
792
793
794
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

795
796
797
798
799
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
800
		     &fetiData, &mat_feti);
801
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
802
			   (void(*)(void))petscMultMatFetiInterface);      
803
    }
804

805
    KSPCreate(mpiCommGlobal, &ksp_feti);
806
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
807
808
809
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
810
    KSPSetFromOptions(ksp_feti);
811
812


813
    // === Create FETI-DP preconditioner object. ===
814

Thomas Witkowski's avatar
Thomas Witkowski committed
815
    if (fetiPreconditioner != FETI_NONE || stokesMode) {
816
817
818
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
819

820

Thomas Witkowski's avatar
Thomas Witkowski committed
821
822
823
824
825
    // === Create null space objects. ===
    
    createNullSpace();


826
    switch (fetiPreconditioner) {
827
828
    case FETI_DIRICHLET:
      TEST_EXIT(meshLevel == 0)
829
830
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

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

834
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
835
836
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
837
838
839
840
841
842
      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);
843
844
845
846
847
848
849
850
851
      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;
      
852
853
      localDofMap.createVec(fetiDirichletPreconData.tmp_vec_b, 
			    nGlobalOverallInterior);
854
855
856
857
858
859
      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));
860

861
862
863
      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

864
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
865
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
866
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
867
	  DegreeOfFreedom d = it->first;	  
868
869
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
870
871
872
873
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

874
875
876
877
      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
878
879
880
881
882
883
884
885


      // 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
886
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
887
888
889
	MSG("FETI-DP timing 08: %.5f seconds (factorization of Dirichlet preconditoner matrices)\n",
	    MPI::Wtime() - wtime);
      }
890
891
892
893
      
      break;

    case FETI_LUMPED:
Thomas Witkowski's avatar
Thomas Witkowski committed
894
895
896
      {
	FetiLumpedPreconData *lumpedData = 
	  (stokesMode ? &fetiInterfaceLumpedPreconData : &fetiLumpedPreconData);
897

Thomas Witkowski's avatar
Thomas Witkowski committed
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
	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;
	  }
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
919

Thomas Witkowski's avatar
Thomas Witkowski committed
920
	if (stokesMode) {
921
922
	  DOFVector<double> tmpvec(pressureFeSpace, "tmpvec");
	  createH2vec(tmpvec);
Thomas Witkowski's avatar
Thomas Witkowski committed
923
	  VtkWriter::writeFile(&tmpvec, "tmpvec.vtu");
924
925
926
927
928
929
930
931
932
933
934
935
936
	  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
937
  
Thomas Witkowski's avatar
Thomas Witkowski committed
938
	  localDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_vec_b1);
Thomas Witkowski's avatar
Thomas Witkowski committed
939
	  primalDofMap.createVec(fetiInterfaceLumpedPreconData.tmp_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
940
941
942
943
944
945
	  fetiInterfaceLumpedPreconData.subSolver = subdomain;	
	}
	
	KSPGetPC(ksp_feti, &precon_feti);
	PCSetType(precon_feti, PCSHELL);
	if (stokesMode) {