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


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/PetscSolverFeti.h"
16
17
#include "parallel/PetscSolverFetiDebug.h"
#include "parallel/PetscSolverFetiMonitor.h"
18
#include "parallel/PetscSolverFetiStructs.h"
19
20
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
21
22
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
23
#include "parallel/PetscSolverGlobalMatrix.h"
24
#include "io/VtkWriter.h"
25
26
27
28
29

namespace AMDiS {

  using namespace std;

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

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

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

    Parameters::get("parallel->multi level test", multiLevelTest);
68
69
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
70
71

    Parameters::get("parallel->print timings", printTimings);
72
73

    Parameters::get("parallel->feti->augmented lagrange", augmentedLagrange);
74
75
76
  }


77
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
78
  {
79
80
    FUNCNAME("PetscSolverFeti::initialize()");

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

84
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
85
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
86

Thomas Witkowski's avatar
Thomas Witkowski committed
87

Thomas Witkowski's avatar
Thomas Witkowski committed
88
89
90
91
92
93
94
95
    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
96
97
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
98
99
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
100

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

112
113
114
115
    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
116

Thomas Witkowski's avatar
Thomas Witkowski committed
117
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
118
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
119

120
    if (fetiPreconditioner == FETI_DIRICHLET) {
121
122
123
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

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


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
129
130
131
132
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
133
134
135
136
137
    bool removeDirichletRows = false;
    Parameters::get("parallel->feti->remove dirichlet", removeDirichletRows);
    if (!removeDirichletRows)
      return;

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    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;
      }
    }
  }


155
156
157
158
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
    double timeCounter = MPI::Wtime();

161
162
163
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
164

165
166
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

167
168
169
170
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
171
    if (fetiPreconditioner == FETI_DIRICHLET)
172
173
      interiorDofMap.clear();

174
175
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
176

177
178
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
179
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
180
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
181
    if (fetiPreconditioner == FETI_DIRICHLET)
182
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
183

184
185
186
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
187
188
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
189
    if (stokesMode) {
190
191
192
193
194
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

195
196
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
197
    
198
199
      createPrimals(feSpace);  

200
201
202
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
203

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
204
      createIndexB(feSpace);     
205
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
206
207
208
209

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
210
    if (fetiPreconditioner == FETI_DIRICHLET)
211
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
212

Thomas Witkowski's avatar
Thomas Witkowski committed
213
    if (stokesMode)
214
215
      interfaceDofMap.update();

216
217
218
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
219
      createAugmentedLagrange(feSpace);
220
    }
221

222
223
224
    lagrangeMap.update();


225
226
227
228
229
230
231
232
233
234
235
236
    // === ===

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

237
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
238
239
240
241
242
243
244
245
246
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
252
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
253
254
255
256
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
257
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
258
	    primalDofMap[feSpace].nRankDofs, 
259
	    primalDofMap[feSpace].nLocalDofs,
Thomas Witkowski's avatar
Thomas Witkowski committed
260
261
262
263
264
265
266
267
268
269
	    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
270
271
272
// 	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
273
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
274
    }
275

Thomas Witkowski's avatar
Thomas Witkowski committed
276
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
279
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
282
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
283
      timeCounter = MPI::Wtime() - timeCounter;
284
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
	  timeCounter);
    }
287
288
289
  }


290
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
291
  {
292
    FUNCNAME("PetscSolverFeti::createPrimals()");  
293

Thomas Witkowski's avatar
Thomas Witkowski committed
294
    if (feSpace == pressureFeSpace)
295
296
      return;

297
298
299
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
300
    // Set of DOF indices that are considered to be primal variables.
301
    DofContainerSet& vertices = 
302
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
303
304

    DofIndexSet primals;
305
    for (DofContainerSet::iterator it = vertices.begin(); 
306
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
307
308
      if (dirichletRows[feSpace].count(**it))
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
309

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
310
311
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
312
      } else {
313
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
316
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
317
318
319
320
321
322
323
	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);
	}
      }
324
    }
325

326
327
328
329

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

330
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
331
332
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
333
      else
334
  	primalDofMap[feSpace].insertNonRankDof(*it);
335
336
337
  }


338
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
339
340
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
341

Thomas Witkowski's avatar
Thomas Witkowski committed
342
    if (feSpace == pressureFeSpace)
343
344
      return;

345
346
347
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
348
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
349
    
350
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
351
352
353
354
355
356
357
358
359
360
361
	 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))
362
	  dualDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
363
364
      }	  
    }
365
366
367
  }

  
368
369
370
371
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
372
    if (feSpace != pressureFeSpace)
373
374
375
376
377
378
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
379
380
381
382
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;      
      
383
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
384
	interfaceDofMap[feSpace].insertRankDof(**it);
385
      else
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
386
387
	interfaceDofMap[feSpace].insertNonRankDof(**it);      
    }
388
389
390
  }


391
392
393
394
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
395
    if (feSpace == pressureFeSpace)
396
397
      return;

398
399
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
402
    // 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
403
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
404

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
405
    if (meshLevel > 0) {
406
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
407
408
409
410
411
412
413
414
415

      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
416
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
417
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
418
	  else
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
	    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
434
435
436
437
438
439
	   !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());
    }
440

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
441
442
443
444
445
446
447
    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)).         ===

448
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
449
450
451
452
453
454
455
    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);

456
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
457
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
458
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
459
460
461
462
	}
      }
    }

463
464
465
466

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

467
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
468

469
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
470
471
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
472
	if (!isPrimal(feSpace, it.getDofIndex()))
473
474
475
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
476
477
478

    stdMpi.updateSendDataSize();

479
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
480
	 !it.end(); it.nextRank()) {
481
      bool recvFromRank = false;
482
      for (; !it.endDofIter(); it.nextDof()) {
483
	if (!isPrimal(feSpace, it.getDofIndex())) {
484
485
486
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
487
488
489
	    recvFromRank = true;
	    break;
	  }
490
	}
491
      }
492
493

      if (recvFromRank)
494
	stdMpi.recv(it.getRank());
495
    }
496

497
498
    stdMpi.startCommunication();

499
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
500
	 !it.end(); it.nextRank()) {
501
      int i = 0;
502
      for (; !it.endDofIter(); it.nextDof())
503
	if (!isPrimal(feSpace, it.getDofIndex()))
504
505
506
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
507
508
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
509
510
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
511
512
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
513

514
515
516
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

517
    int nRankLagrange = 0;
518
    DofMap& dualMap = dualDofMap[feSpace].getMap();
519
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
520
521
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
522
	int degree = boundaryDofRanks[feSpace][it->first].size();
523
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
524
      } else {
525
	lagrangeMap[feSpace].insertNonRankDof(it->first);
526
527
      }
    }
528
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
529
530
531
  }


532
533
534
535
536
537
538
539
540
  void PetscSolverFeti::createAugmentedLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


541
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
542
  {
543
    FUNCNAME("PetscSolverFeti::createIndexB()");
544

545
    DOFAdmin* admin = feSpace->getAdmin();
546
547
548
549

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

551
    int nLocalInterior = 0;    
552
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
553
554
555
556
557
558
      if (admin->isDofFree(i) ||
	  isPrimal(feSpace, i) ||
	  isDual(feSpace, i) ||
	  isInterface(feSpace, i) ||
	  dirichletRows[feSpace].count(i))
	continue;      
559

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
      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");	
575
      }
576
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
577
    
578
579
    // === And finally, add the global indicies of all dual nodes. ===

580
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
581
582
	 it != dualDofMap[feSpace].getMap().end(); ++it) {
      if (meshLevel == 0) {
583
	localDofMap[feSpace].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
584
      } else {
585
586
587
588
589
	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
590
    }
591
592
593
  }


594
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
595
596
597
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
598
    double wtime = MPI::Wtime();
599
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
600

601
602
    // === Create distributed matrix for Lagrange constraints. ===

603
604
605
606
607
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
608
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
609

610
611
612
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

613
614
615
616
617
618
619
    // === 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
620
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
621
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
622

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

Thomas Witkowski's avatar
Thomas Witkowski committed
630
	// Copy set of all ranks that contain this dual node.
631
632
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
633
634
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
635
636

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
637
	
638
	int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
639
640
641
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
642
643
644
645
646
	      MatSetValue(mat_lagrange, 
			  index + counter, 
			  localDofMap.getMatIndex(k, it->first) + rStartInterior,
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
647
	    }
648
649
650
651
652
653
654
655
656
657
658
659
660
661
	    counter++;
	  }
	}

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
	if (meshDistributor->getDofMap()[feSpaces[k]].isRankDof(it->first)) {
	  int nConstraints = (degree * (degree - 1)) / 2;
	  for (int i = 0; i < nConstraints; i++) {
	    VecSetValue(vec_scale_lagrange,
			index + i,
			1.0 / static_cast<double>(degree),
			INSERT_VALUES);
662
663
664
665
666
667
668
	  }
	}
      }
    }

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

670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687

    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

    VecView(vec_scale_lagrange, PETSC_VIEWER_STDOUT_WORLD);

    if (fetiPreconditioner != FETI_NONE || stokesMode) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatDiagonalScale(mat_lagrange_scaled, vec_scale_lagrange, PETSC_NULL);
    }

    VecDestroy(&vec_scale_lagrange);


    // === Print final timings. ===

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
688
689
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
690
691
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
692
    }
693
694
695
  }


696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
  void PetscSolverFeti::createMatAugmentedLagrange(vector<const FiniteElemSpace*> &feSpaces)
  {
    FUNCNAME("PetscSolverFeti::createMatAugmentedLagrange()");

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

    nRankEdges = 0;
    nOverallEdges = 0;
    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
      if (it->rankObj.subObj == EDGE)
	nRankEdges++;
    
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

    MSG("nRankEdges = %d, nOverallEdges = %d\n", nRankEdges, nOverallEdges);

    nRankEdges *= feSpaces.size();
    rStartEdges *= feSpaces.size();
    nOverallEdges *= feSpaces.size();

    MatCreateAIJ(mpiCommGlobal,
722
723
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
724
725
726
727
		 1, PETSC_NULL, 1, PETSC_NULL, 
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

728
    int rowCounter = rStartEdges;
729
730
731
732
733
734
735
736
737
738
739
740
741
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
      if (it->rankObj.subObj == EDGE) {
	for (int i = 0; i < feSpaces.size(); i++) {
	  DofContainer edgeDofs;
	  it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj, edgeDofs);

	  MSG("SIZE = %d\n", edgeDofs.size());

	  for (DofContainer::iterator it = edgeDofs.begin();
	       it != edgeDofs.end(); it++) {
	    TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false)
	      ("Should not be primal!\n");

742
	    int col = lagrangeMap.getMatIndex(i, **it);
743
	    double value = 1.0;
744
	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
745
746
	  }

747
	  rowCounter++;
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
	}
      }
    }

    MatAssemblyBegin(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);

    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 05a: %.5f seconds (creation of augmented lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
    }
  }


763
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
764
  {
765
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
766

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

770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
      if (augmentedLagrange == false) {
	schurPrimalData.subSolver = subdomain;
	
	localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
	
	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getOverallDofs(), 
		       primalDofMap.getOverallDofs(),
		       &schurPrimalData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimal);	
      } else {
	schurPrimalAugmentedData.subSolver = subdomain;

	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b0, nGlobalOverallInterior);
	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b1, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalAugmentedData.tmp_vec_primal);
	lagrangeMap.createVec(schurPrimalAugmentedData.tmp_vec_lagrange);

	schurPrimalAugmentedData.mat_lagrange = &mat_lagrange;
	schurPrimalAugmentedData.mat_augmented_lagrange = &mat_augmented_lagrange;

	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges,
		       &schurPrimalAugmentedData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimalAugmented);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
806

807
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
808
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
809
810
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
811
812
      KSPSetFromOptions(ksp_schur_primal);
    } else {
813
814
      MSG("Create direct schur primal solver!\n");

815
816
      TEST_EXIT(!augmentedLagrange)("Not yet supported!\n");

817
818
      double wtime = MPI::Wtime();

819
820
821
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

822
823
824
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
825

Thomas Witkowski's avatar
Thomas Witkowski committed
826
      Mat matBPi;
827
828
829
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
830
831
832
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

833

834
835
836
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
837

838
839
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

840
      for (int i = 0; i < nRowsRankB; i++) {
841
	PetscInt row = localDofMap.getStartDofs() + i;
842
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
843
844
845
846
847

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

848
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
849
850
      }

851
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
852
		primalDofMap.getLocalDofs())
853
	("Should not happen!\n");
854

855
856
857
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
858
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
859
860
861
862
863
864
865

 	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
866

Thomas Witkowski's avatar
Thomas Witkowski committed
867
	subdomain->solve(tmpVec, tmpVec);
868

Thomas Witkowski's avatar
Thomas Witkowski committed
869
	PetscScalar *tmpValues;
870
	VecGetArray(tmpVec, &tmpValues);
871
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
872
	  MatSetValue(matBPi, 
873
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
874
875
876
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
877
878
879
880
881
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
882
883
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
884

885
886
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
887
888

      Mat matPrimal;
889
890
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
891
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
892
893

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
894
      MatDestroy(&matBPi);
895

896
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
897
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
898
		      SAME_NONZERO_PATTERN);
899
900
901
902
903
904
      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);
905
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
906

Thomas Witkowski's avatar
Thomas Witkowski committed
907
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
908
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
909
910
911
912
913
914
915
916
917
918
	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
919
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
920
921
922
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
923
    }
924
925
926
927
928
929
930
  }


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

931
    if (schurPrimalSolver == 0) {
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
      if (augmentedLagrange == false) {
	schurPrimalData.subSolver = NULL;
	
	VecDestroy(&schurPrimalData.tmp_vec_b);
	VecDestroy(&schurPrimalData.tmp_vec_primal);
      } else {
	schurPrimalAugmentedData.subSolver = NULL;
	schurPrimalAugmentedData.mat_lagrange = NULL;
	schurPrimalAugmentedData.mat_augmented_lagrange = NULL;

	VecDestroy(&schurPrimalAugmentedData.tmp_vec_b0);
	VecDestroy(&schurPrimalAugmentedData.tmp_vec_b1);
	VecDestroy(&schurPrimalAugmentedData.tmp_vec_primal);
	VecDestroy(&schurPrimalAugmentedData.tmp_vec_lagrange);
      }
947
    }
948
949
950
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
951
952
953
  }


954
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
955
956
957
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

958
959
    // === Create FETI-DP solver object. ===

960
961
962
963
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
964
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
965
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
966
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
967

Thomas Witkowski's avatar
Thomas Witkowski committed
968
    if (stokesMode == false) {      
969
970
971
972
973
974
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
975
976
977
978
979
980
981
982
983
      if (augmentedLagrange == false) {
	MatShellSetOperation(mat_feti, MATOP_MULT, 
			     (void(*)(void))petscMultMatFeti);
      } else {
	fetiData.mat_augmented_lagrange = &mat_augmented_lagrange;
	primalDofMap.createVec(fetiData.tmp_vec_primal1);
	MatShellSetOperation(mat_feti, MATOP_MULT, 
			     (void(*)(void))petscMultMatFetiAugmented);	
      }
984
    }  else {
985
986
      TEST_EXIT_DBG(!augmentedLagrange)("Not yet implemented!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
987
988
989
990
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

991
992
993
994
995
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
996
		     &fetiData, &mat_feti);
997
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
998
			   (void(*)(void))petscMultMatFetiInterface);      
999
    }
1000

1001
    KSPCreate(mpiCommGlobal, &ksp_feti);
1002
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
1003
1004
1005
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
1006
    KSPSetFromOptions(ksp_feti);
1007
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
1008
      KSPMonitorSet(ksp_feti, KSPMonitorFetiStokes, &fetiKspData, PETSC_NULL);
1009

Thomas Witkowski's avatar
Thomas Witkowski committed
1010
1011
1012
1013
1014
    // === Create null space objects. ===
    
    createNullSpace();