PetscSolverFeti.cc 76.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;

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
  PetscErrorCode KSPMonitorFeti(KSP ksp, PetscInt n, PetscReal rnorm, void *dummy)
  {    
    Vec            Br,v,w;
    Mat            A;

    KSPGetOperators(ksp, &A, PETSC_NULL, PETSC_NULL);
    MatGetVecs(A, &w, &v);
    KSPBuildResidual(ksp, v, w, &Br);

    Vec nest0, nest1;
    VecNestGetSubVec(Br, 0, &nest0);
    VecNestGetSubVec(Br, 1, &nest1);

    PetscScalar norm0, norm1;
    VecNorm(nest0, NORM_2, &norm0);
    VecNorm(nest1, NORM_2, &norm1);

    VecDestroy(&v);
    VecDestroy(&w);

    PetscPrintf(PETSC_COMM_WORLD, "%3D KSP Component residual norm [%1.12e %1.12e] and relative [%1.12e]\n",
		norm0, norm1, rnorm);

    return 0;
  }

54
55
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
56
      schurPrimalSolver(0),
57
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
58
      subdomain(NULL),
59
60
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
61
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
62
      printTimings(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
63
64
      stokesMode(false),
      pressureComponent(-1)
65
66
67
68
69
70
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
71
      MSG("Create FETI-DP solver with no preconditioner!\n");
72
73
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
74
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
75
76
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
77
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
78
79
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
80
81
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
82
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
83
84
85
86
87

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

    Parameters::get("parallel->multi level test", multiLevelTest);
90
91
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
92
93

    Parameters::get("parallel->print timings", printTimings);
94
95
96
  }


97
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
98
  {
99
100
    FUNCNAME("PetscSolverFeti::initialize()");

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

104
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
105
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
106

Thomas Witkowski's avatar
Thomas Witkowski committed
107

Thomas Witkowski's avatar
Thomas Witkowski committed
108
109
110
111
112
113
114
115
    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
116
117
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
118
119
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
120

121
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
122
	subdomain->setMeshDistributor(meshDistributor, 
123
				      mpiCommGlobal, mpiCommLocal);
124
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
125
	subdomain->setMeshDistributor(meshDistributor, 
126
127
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
128
	subdomain->setLevel(meshLevel);
129
130
      }
    }
131

132
133
134
135
    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
136

Thomas Witkowski's avatar
Thomas Witkowski committed
137
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
138
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
139

140
    if (fetiPreconditioner == FETI_DIRICHLET) {
141
142
143
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

144
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
145
    }
146
147
148
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
149
150
151
152
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
153
154
155
156
157
    bool removeDirichletRows = false;
    Parameters::get("parallel->feti->remove dirichlet", removeDirichletRows);
    if (!removeDirichletRows)
      return;

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
    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;
      }
    }
  }


175
176
177
178
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
    double timeCounter = MPI::Wtime();

181
182
183
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
184

185
186
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

187
188
189
190
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
191
    if (fetiPreconditioner == FETI_DIRICHLET)
192
193
      interiorDofMap.clear();

194
195
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
196

197
198
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
199
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
200
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
201
    if (fetiPreconditioner == FETI_DIRICHLET)
202
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
203

204
205
206
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
207
208
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
209
    if (stokesMode) {
210
211
212
213
214
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

215
216
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
217
    
218
219
      createPrimals(feSpace);  

220
221
222
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
223

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
224
      createIndexB(feSpace);     
225
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
226
227
228
229

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
230
    if (fetiPreconditioner == FETI_DIRICHLET)
231
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
232

Thomas Witkowski's avatar
Thomas Witkowski committed
233
    if (stokesMode)
234
235
      interfaceDofMap.update();

236
237
238
239
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
240

241
242
243
    lagrangeMap.update();


244
245
246
247
248
249
250
251
252
253
254
255
    // === ===

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

256
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
257
258
259
260
261
262
263
264
265
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
271
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
	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
288
289
290
// 	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
291
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
292
    }
293

Thomas Witkowski's avatar
Thomas Witkowski committed
294
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
295
296
297
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
300
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
301
      timeCounter = MPI::Wtime() - timeCounter;
302
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
	  timeCounter);
    }
305
306
307
  }


308
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
309
  {
310
    FUNCNAME("PetscSolverFeti::createPrimals()");  
311

Thomas Witkowski's avatar
Thomas Witkowski committed
312
    if (feSpace == pressureFeSpace)
313
314
      return;

315
316
317
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
318
    // Set of DOF indices that are considered to be primal variables.
319
    DofContainerSet& vertices = 
320
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
321
322

    DofIndexSet primals;
323
    for (DofContainerSet::iterator it = vertices.begin(); 
324
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
325
326
      if (dirichletRows[feSpace].count(**it))
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
327

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
328
329
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
330
      } else {
331
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
334
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
335
336
337
338
339
340
341
	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);
	}
      }
342
    }
343

344
345
346
347

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

348
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
349
350
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
351
      else
352
  	primalDofMap[feSpace].insertNonRankDof(*it);
353
354
355
  }


356
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
357
358
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
359

Thomas Witkowski's avatar
Thomas Witkowski committed
360
    if (feSpace == pressureFeSpace)
361
362
      return;

363
364
365
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
366
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
367
    
368
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
369
370
371
372
373
374
375
376
377
378
379
	 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))
380
	  dualDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
381
382
      }	  
    }
383
384
385
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    if (feSpace != pressureFeSpace)
391
392
393
394
395
396
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
397
398
399
400
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;      
      
401
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
402
	interfaceDofMap[feSpace].insertRankDof(**it);
403
      else
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
404
405
	interfaceDofMap[feSpace].insertNonRankDof(**it);      
    }
406
407
408
  }


409
410
411
412
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
413
    if (feSpace == pressureFeSpace)
414
415
      return;

416
417
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
418
419
420
    // 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
421
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
422

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
423
    if (meshLevel > 0) {
424
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
425
426
427
428
429
430
431
432
433

      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
434
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
435
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
436
	  else
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
	    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
452
453
454
455
456
457
	   !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());
    }
458

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
459
460
461
462
463
464
465
    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)).         ===

466
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
467
468
469
470
471
472
473
    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);

474
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
475
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
476
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
477
478
479
480
	}
      }
    }

481
482
483
484

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

485
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
486

487
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
488
489
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
490
	if (!isPrimal(feSpace, it.getDofIndex()))
491
492
493
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
494
495
496

    stdMpi.updateSendDataSize();

497
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
498
	 !it.end(); it.nextRank()) {
499
      bool recvFromRank = false;
500
      for (; !it.endDofIter(); it.nextDof()) {
501
	if (!isPrimal(feSpace, it.getDofIndex())) {
502
503
504
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
505
506
507
	    recvFromRank = true;
	    break;
	  }
508
	}
509
      }
510
511

      if (recvFromRank)
512
	stdMpi.recv(it.getRank());
513
    }
514

515
516
    stdMpi.startCommunication();

517
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
518
	 !it.end(); it.nextRank()) {
519
      int i = 0;
520
      for (; !it.endDofIter(); it.nextDof())
521
	if (!isPrimal(feSpace, it.getDofIndex()))
522
523
524
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
525
526
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
527
528
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
529
530
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
531

532
533
534
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

535
    int nRankLagrange = 0;
536
    DofMap& dualMap = dualDofMap[feSpace].getMap();
537
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
538
539
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
540
	int degree = boundaryDofRanks[feSpace][it->first].size();
541
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
542
      } else {
543
	lagrangeMap[feSpace].insertNonRankDof(it->first);
544
545
      }
    }
546
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
547
548
549
  }


550
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
551
  {
552
    FUNCNAME("PetscSolverFeti::createIndexB()");
553

554
    DOFAdmin* admin = feSpace->getAdmin();
555
556
557
558

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

560
    int nLocalInterior = 0;    
561
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
562
563
564
565
566
567
      if (admin->isDofFree(i) ||
	  isPrimal(feSpace, i) ||
	  isDual(feSpace, i) ||
	  isInterface(feSpace, i) ||
	  dirichletRows[feSpace].count(i))
	continue;      
568

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
      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");	
584
      }
585
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
586
    
587
588
    // === And finally, add the global indicies of all dual nodes. ===

589
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
590
591
	 it != dualDofMap[feSpace].getMap().end(); ++it) {
      if (meshLevel == 0) {
592
	localDofMap[feSpace].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
593
      } else {
594
595
596
597
598
	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
599
    }
600
601
602
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
607
    double wtime = MPI::Wtime();
608
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
609

610
611
    // === Create distributed matrix for Lagrange constraints. ===

612
613
614
615
616
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
617
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
618

619
620
621
622
623
624
625
    // === 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
626
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
627
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
628

629
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
630
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
631
632
633
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
634
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
635

Thomas Witkowski's avatar
Thomas Witkowski committed
636
	// Copy set of all ranks that contain this dual node.
637
638
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
639
640
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
641
642

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
643
644
645
646
	
	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
647
648
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
649

650
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
651
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
652
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
653
	    index++;	      
654
655
656
657
658
659
660
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
662
663
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
666
    }
667
668
669
  }


670
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
671
  {
672
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
673

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

Thomas Witkowski's avatar
Thomas Witkowski committed
677
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
678
      
679
680
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
681

682
      MatCreateShell(mpiCommGlobal,
683
684
685
686
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
687
688
689
690
691
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
692
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
693
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
694
695
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
696
697
      KSPSetFromOptions(ksp_schur_primal);
    } else {
698
699
      MSG("Create direct schur primal solver!\n");

700
701
      double wtime = MPI::Wtime();

702
703
704
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

705
706
707
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
708

Thomas Witkowski's avatar
Thomas Witkowski committed
709
      Mat matBPi;
710
711
712
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
713
714
715
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

716

717
718
719
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
720

721
722
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

723
      for (int i = 0; i < nRowsRankB; i++) {
724
	PetscInt row = localDofMap.getStartDofs() + i;
725
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
726
727
728
729
730

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

731
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
732
733
      }

734
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
735
		primalDofMap.getLocalDofs())
736
	("Should not happen!\n");
737

738
739
740
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
741
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
742
743
744
745
746
747
748

 	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
749

Thomas Witkowski's avatar
Thomas Witkowski committed
750
	subdomain->solve(tmpVec, tmpVec);
751

Thomas Witkowski's avatar
Thomas Witkowski committed
752
	PetscScalar *tmpValues;
753
	VecGetArray(tmpVec, &tmpValues);
754
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
755
	  MatSetValue(matBPi, 
756
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
757
758
759
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
760
761
762
763
764
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
765
766
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
767

768
769
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
770
771

      Mat matPrimal;
772
773
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
774
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
775
776

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
777
      MatDestroy(&matBPi);
778

779
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
780
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
781
		      SAME_NONZERO_PATTERN);
782
783
784
785
786
787
      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);
788
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
789

Thomas Witkowski's avatar
Thomas Witkowski committed
790
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
791
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
792
793
794
795
796
797
798
799
800
801
	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
802
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
803
804
805
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
806
    }
807
808
809
810
811
812
813
  }


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

814
    if (schurPrimalSolver == 0) {
815
      schurPrimalData.subSolver = NULL;
816

817
818
819
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
820
821
822
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
823
824
825
  }


826
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
827
828
829
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

830
831
    // === Create FETI-DP solver object. ===

832
833
834
835
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
836
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
837
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
838
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
839

Thomas Witkowski's avatar
Thomas Witkowski committed
840
    if (stokesMode == false) {      
841
842
843
844
845
846
847
848
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
849
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
850
851
852
853
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

854
855
856
857
858
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
859
		     &fetiData, &mat_feti);
860
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
861
			   (void(*)(void))petscMultMatFetiInterface);      
862
    }
863

864
    KSPCreate(mpiCommGlobal, &ksp_feti);
865
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
866
867
868
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
869
    KSPSetFromOptions(ksp_feti);
870
871


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
872
873


874
    // === Create FETI-DP preconditioner object. ===
875

Thomas Witkowski's avatar
Thomas Witkowski committed
876
    if (fetiPreconditioner != FETI_NONE || stokesMode) {
877
878
879
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
880

881

Thomas Witkowski's avatar
Thomas Witkowski committed
882
883
884
885
886
    // === Create null space objects. ===
    
    createNullSpace();


887
    switch (fetiPreconditioner) {
888
889
    case FETI_DIRICHLET:
      TEST_EXIT(meshLevel == 0)
890
891
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

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

895
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
896
897
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
898
899
900
901
902
903
      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);
904
905
906
907
908
909
910
911
912
      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;
      
913
914
      localDofMap.createVec(fetiDirichletPreconData.tmp_vec_b, 
			    nGlobalOverallInterior);
915
916
917
918
919
920
      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));
921

922
923
924
      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

925
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
926
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
927
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
928
	  DegreeOfFreedom d = it->first;	  
929
930
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
931
932
933
934
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

935
936
937
938
      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
939
940
941
942
943
944
945
946


      // 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
947
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
948
949
950
	MSG("FETI-DP timing 08: %.5f seconds (factorization of Dirichlet preconditoner matrices)\n",
	    MPI::Wtime() - wtime);
      }
951
952
953
954
      
      break;

    case FETI_LUMPED:
Thomas Witkowski's avatar
Thomas Witkowski committed
955
956
957
      {
	FetiLumpedPreconData *lumpedData = 
	  (stokesMode ? &fetiInterfaceLumpedPreconData : &fetiLumpedPreconData);
958

Thomas Witkowski's avatar
Thomas Witkowski committed
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976