PetscSolverFeti.cc 59.4 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
37
      printTimings(false),
      enableStokesMode(false)
38
39
40
41
42
43
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

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

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


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
80
81
82
83
84
85

    enableStokesMode = false;
    Parameters::get("parallel->stokes mode", enableStokesMode);
    if (enableStokesMode) {
      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");
86
      fullInterface = feSpaces[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
87
88
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
89
90
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
91

92
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
93
	subdomain->setMeshDistributor(meshDistributor, 
94
				      mpiCommGlobal, mpiCommLocal);
95
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
96
	subdomain->setMeshDistributor(meshDistributor, 
97
98
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
99
	subdomain->setLevel(meshLevel);
100
101
      }
    }
102

103
104
105
106
    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
107

108
    if (fullInterface != NULL)
Thomas Witkowski's avatar
Thomas Witkowski committed
109
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
110

111
    if (fetiPreconditioner == FETI_DIRICHLET) {
112
113
114
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

115
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
116
    }
117
118
119
120
121
122
123
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
124
125
    double timeCounter = MPI::Wtime();

126
127
128
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
129

130
131
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

132
133
134
135
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
136
    if (fetiPreconditioner == FETI_DIRICHLET)
137
138
      interiorDofMap.clear();

139
140
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
141

142
143
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
144
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
145
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
146
    if (fetiPreconditioner == FETI_DIRICHLET)
147
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
148

149
150
151
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
152
153
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

154
155
156
157
158
159
    if (fullInterface != NULL) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

160
161
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
162
    
163
164
      createPrimals(feSpace);  

165
166
167
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
168

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
169
      createIndexB(feSpace);     
170
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
171
172
173
174

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
175
    if (fetiPreconditioner == FETI_DIRICHLET)
176
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
177

178
179
180
    if (fullInterface != NULL)
      interfaceDofMap.update();

181
182
183
184
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
185

186
187
188
    lagrangeMap.update();


189
190
191
192
193
194
195
196
197
198
199
200
    // === ===

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

201
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
202
203
204
205
206
207
208
209
210
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
      if (feSpace == fullInterface) {
	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
237
    }
238

Thomas Witkowski's avatar
Thomas Witkowski committed
239
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
240
241
242

    if (enableStokesMode) {
      MSG("WARNING: MAKE THIS MORE GENERAL!!!!!\n");
Thomas Witkowski's avatar
So    
Thomas Witkowski committed
243
244
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, 0);
      //      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
245
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 1);
Thomas Witkowski's avatar
So    
Thomas Witkowski committed
246
247
248
249
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 2);
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 3);
      subdomain->setCoarseSpaceDofMapping(&primalDofMap, 4);

Thomas Witkowski's avatar
Thomas Witkowski committed
250
251
252
    } else {	  
      subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
253
254

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


263
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
264
  {
265
    FUNCNAME("PetscSolverFeti::createPrimals()");  
266

267
268
269
    if (feSpace == fullInterface)
      return;

270
271
272
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
273
    // Set of DOF indices that are considered to be primal variables.
274
    DofContainerSet& vertices = 
275
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
276
277

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

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

296
297
298
299

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

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


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

312
313
314
    if (feSpace == fullInterface)
      return;

315
316
317
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
318
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
319
320
321
322

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

  
332
333
334
335
336
337
338
339
340
341
342
343
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

    if (feSpace != fullInterface)
      return;

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

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


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

355
356
357
    if (feSpace == fullInterface)
      return;

358
359
    boundaryDofRanks[feSpace].clear();

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

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

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

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

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

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

423
424
425
426

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

427
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
428

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
454
	stdMpi.recv(it.getRank());
455
    }
456

457
458
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
473

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

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


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

496
    DOFAdmin* admin = feSpace->getAdmin();
497
498
499
500

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

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

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

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

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


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

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

549
550
    // === Create distributed matrix for Lagrange constraints. ===

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

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

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

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

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

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

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

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


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

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

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

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

639
640
      double wtime = MPI::Wtime();

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

644
645
646
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
647

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

655

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

660
661
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

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

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

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

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

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

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

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
704
705
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
706

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

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

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
716
      MatDestroy(&matBPi);
717

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

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


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

753
    if (schurPrimalSolver == 0) {
754
      schurPrimalData.subSolver = NULL;
755

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


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

769
770
    // === Create FETI-DP solver object. ===

771
772
773
774
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
    if (enableStokesMode) {
      Vec nullSpaceInterface;
      Vec nullSpaceLagrange;
      Vec nullSpaceBasis;

      interfaceDofMap.createVec(nullSpaceInterface);
      VecSet(nullSpaceInterface, 1.0);

      lagrangeMap.createVec(nullSpaceLagrange);
      VecSet(nullSpaceLagrange, 0.0);

      Vec vecArray[2] = {nullSpaceInterface, nullSpaceLagrange};
      VecCreateNest(mpiCommGlobal, 2, PETSC_NULL, vecArray, &nullSpaceBasis);

      
      MatNullSpace matNullspace;
      MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullspace);
827
828
      MatSetNullSpace(mat_feti, matNullspace);
      KSPSetNullSpace(ksp_feti, matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
829
830
    }

831

832
    // === Create FETI-DP preconditioner object. ===
833

834
835
836
837
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
838

839

840
    switch (fetiPreconditioner) {
841
842
    case FETI_DIRICHLET:
      TEST_EXIT(meshLevel == 0)
843
844
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

845
846
847
      TEST_EXIT(!enableStokesMode)
	("Stokes mode does not yet support the Dirichlet precondition!\n");

848
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
849
850
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
851
852
853
854
855
856
      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);
857
858
859
860
861
862
863
864
865
      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;
      
866
867
      localDofMap.createVec(fetiDirichletPreconData.tmp_vec_b, 
			    nGlobalOverallInterior);
868
869
870
871
872
873
      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));
874

875
876
877
      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

878
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
879
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
880
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
881
	  DegreeOfFreedom d = it->first;	  
882
883
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
884
885
886
887
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

888
889
890
891
      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
892
893
894
895
896
897
898
899


      // 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
900
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
901
902
903
	MSG("FETI-DP timing 08: %.5f seconds (factorization of Dirichlet preconditoner matrices)\n",
	    MPI::Wtime() - wtime);
      }
904
905
906
907
      
      break;

    case FETI_LUMPED:
908
      fetiLumpedPreconData.enableStokesMode = enableStokesMode;
909
910
911
      fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;

Thomas Witkowski's avatar
Thomas Witkowski committed
912
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
913
914
915
	if (enableStokesMode && feSpaces[i] == fullInterface)
	  continue;

916
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
917
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
918
	  DegreeOfFreedom d = it->first;
919
920
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
Thomas Witkowski's avatar
Thomas Witkowski committed
921
922
923
924
	  fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

925
      localDofMap.createVec(fetiLumpedPreconData.tmp_vec_b);
926
927
928
929
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
930
931
932
933
934
935

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
936
937
      break;
    default:
938
939
      break;
    }
940
941
942
943
944
945
946
  }
  

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

947
948
    // === Destroy FETI-DP solver object. ===

949
    fetiData.mat_lagrange = PETSC_NULL;
950
    fetiData.subSolver = NULL;
951
    fetiData.ksp_schur_primal = PETSC_NULL;