PetscSolverFeti.cc 77.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/PetscSolverFeti.h"
16
#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
  PetscErrorCode KSPMonitorFeti(KSP ksp, PetscInt n, PetscReal rnorm, void *data)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
29
  {    
30
31
32
    Vec Br,v,w;
    VecDuplicate(static_cast<FetiKspData*>(data)->draft, &v);
    VecDuplicate(static_cast<FetiKspData*>(data)->draft, &w);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
33
34
35
36
37
38
39

    KSPBuildResidual(ksp, v, w, &Br);

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

40
41
    PetscScalar norm, norm0, norm1;
    VecNorm(Br, NORM_2, &norm);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
42
43
44
45
46
47
    VecNorm(nest0, NORM_2, &norm0);
    VecNorm(nest1, NORM_2, &norm1);

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

48
49
    PetscPrintf(PETSC_COMM_WORLD, "%3D KSP residual norm %1.12e [%1.12e %1.12e] and preconditioned norm [%1.12e]\n",
		n, norm, norm0, norm1, rnorm);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
50
51
52
53

    return 0;
  }

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

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

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

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

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


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
108

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

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

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

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

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

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


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

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

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


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

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

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

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

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

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

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

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

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

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

221
222
223
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
224

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

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

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

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

242
243
244
    lagrangeMap.update();


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

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

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

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

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

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

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


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

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

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

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

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

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

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

345
346
347
348

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

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


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

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

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

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

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

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

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

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


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

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

417
418
    boundaryDofRanks[feSpace].clear();

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

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

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

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

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

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

482
483
484
485

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

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

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

    stdMpi.updateSendDataSize();

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

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

516
517
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
532

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

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


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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

717

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

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

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

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

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

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

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

 	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
750

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

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

	VecDestroy(&tmpVec);
      }

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

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

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

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

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

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


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

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

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


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

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

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

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

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

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

865
    KSPCreate(mpiCommGlobal, &ksp_feti);
866
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
867
868
869
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
870
    KSPSetFromOptions(ksp_feti);
871
872
    if (stokesMode)
      KSPMonitorSet(ksp_feti, KSPMonitorFeti, &fetiKspData, PETSC_NULL);
873

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
874
875


876
    // === Create FETI-DP preconditioner object. ===
877

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

883

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


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

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

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

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

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

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


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