PetscSolverFeti.cc 71.1 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
    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)) {
147
148
149
150
151
152
153
154
155
	// === Run test if Dirichlet rows are all the same in all components ===
	// === of one FE space.                                              ===
	TEST_EXIT(dirichletRows[feSpace].size() == dRows.size())
	  ("Wrong number of dirichlet rows!\n");
	
	for (std::set<DegreeOfFreedom>::iterator it0 = dirichletRows[feSpace].begin(),
	       it1 = dRows.begin(); it1 != dRows.end(); ++it0, ++it1) {
	  TEST_EXIT(*it0 == *it1)("Wrong DOFs %d %d!\n", *it0, *it1);
	}
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
156
157
158
159
160
161
162
      } else {
	dirichletRows[feSpace] = dRows;
      }
    }
  }


163
164
165
166
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
    double timeCounter = MPI::Wtime();

169
170
171
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
172

173
174
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

175
176
177
178
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
179
    if (fetiPreconditioner == FETI_DIRICHLET)
180
181
      interiorDofMap.clear();

182
183
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
184

185
186
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
187
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
188
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
189
    if (fetiPreconditioner == FETI_DIRICHLET)
190
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
191

192
193
194
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
195
196
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
197
    if (stokesMode) {
198
199
200
201
202
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

203
204
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
205
    
206
207
      createPrimals(feSpace);  

208
209
210
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
211

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
212
      createIndexB(feSpace);     
213
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
216
217

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
218
    if (fetiPreconditioner == FETI_DIRICHLET)
219
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
220

Thomas Witkowski's avatar
Thomas Witkowski committed
221
    if (stokesMode)
222
223
      interfaceDofMap.update();

224
225
226
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
227
      createAugmentedLagrange(feSpace);
228
    }
229

230
231
232
    lagrangeMap.update();


233
234
235
236
237
238
239
240
241
242
243
244
    // === ===

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

245
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
246
247
248
249
250
251
252
253
254
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
260
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
261
262
263
264
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
265
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
266
	    primalDofMap[feSpace].nRankDofs, 
267
	    primalDofMap[feSpace].nLocalDofs,
Thomas Witkowski's avatar
Thomas Witkowski committed
268
269
270
271
272
273
274
275
276
277
	    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
278
279
280
// 	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
281
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
282
    }
283

Thomas Witkowski's avatar
Thomas Witkowski committed
284
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
287
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
288
289

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
290
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
291
      timeCounter = MPI::Wtime() - timeCounter;
292
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
	  timeCounter);
    }
295
296
297
  }


298
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
299
  {
300
    FUNCNAME("PetscSolverFeti::createPrimals()");  
301

Thomas Witkowski's avatar
Thomas Witkowski committed
302
    if (feSpace == pressureFeSpace)
303
304
      return;

305
306
307
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
308
    // Set of DOF indices that are considered to be primal variables.
309
    DofContainerSet& vertices = 
310
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
311
312

    DofIndexSet primals;
313
    for (DofContainerSet::iterator it = vertices.begin(); 
314
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
315
316
      if (dirichletRows[feSpace].count(**it))
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
317

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
318
319
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
320
      } else {
321
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
324
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
325
326
327
328
329
330
331
	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);
	}
      }
332
    }
333

334
335
336
337

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

338
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
339
340
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
341
      else
342
  	primalDofMap[feSpace].insertNonRankDof(*it);
343
344
345
  }


346
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
347
348
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
349

Thomas Witkowski's avatar
Thomas Witkowski committed
350
    if (feSpace == pressureFeSpace)
351
352
      return;

353
354
355
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
356
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
357
    
358
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
359
360
361
362
363
364
365
366
367
368
369
	 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))
370
	  dualDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
371
372
      }	  
    }
373
374
375
  }

  
376
377
378
379
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
380
    if (feSpace != pressureFeSpace)
381
382
383
384
385
386
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
387
388
389
390
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;      
      
391
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
392
	interfaceDofMap[feSpace].insertRankDof(**it);
393
      else
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
394
395
	interfaceDofMap[feSpace].insertNonRankDof(**it);      
    }
396
397
398
  }


399
400
401
402
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
403
    if (feSpace == pressureFeSpace)
404
405
      return;

406
407
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
410
    // 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
411
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
412

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
413
    if (meshLevel > 0) {
414
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
415
416
417
418
419
420
421
422
423

      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
424
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
425
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
426
	  else
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
	    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
442
443
444
445
446
447
	   !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());
    }
448

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
449
450
451
452
453
454
455
    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)).         ===

456
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
457
458
459
460
461
462
463
    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);

464
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
465
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
466
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
467
468
469
470
	}
      }
    }

471
472
473
474

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

475
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
476

477
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
478
479
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
480
	if (!isPrimal(feSpace, it.getDofIndex()))
481
482
483
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
484
485
486

    stdMpi.updateSendDataSize();

487
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
488
	 !it.end(); it.nextRank()) {
489
      bool recvFromRank = false;
490
      for (; !it.endDofIter(); it.nextDof()) {
491
	if (!isPrimal(feSpace, it.getDofIndex())) {
492
493
494
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
495
496
497
	    recvFromRank = true;
	    break;
	  }
498
	}
499
      }
500
501

      if (recvFromRank)
502
	stdMpi.recv(it.getRank());
503
    }
504

505
506
    stdMpi.startCommunication();

507
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
508
	 !it.end(); it.nextRank()) {
509
      int i = 0;
510
      for (; !it.endDofIter(); it.nextDof())
511
	if (!isPrimal(feSpace, it.getDofIndex()))
512
513
514
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
515
516
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
517
518
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
519
520
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
521

522
523
524
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

525
    int nRankLagrange = 0;
526
    DofMap& dualMap = dualDofMap[feSpace].getMap();
527
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
528
529
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
530
	int degree = boundaryDofRanks[feSpace][it->first].size();
531
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
532
      } else {
533
	lagrangeMap[feSpace].insertNonRankDof(it->first);
534
535
      }
    }
536
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
537
538
539
  }


540
541
542
543
544
545
546
547
548
  void PetscSolverFeti::createAugmentedLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


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

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

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

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

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

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


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

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

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

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

618
619
620
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

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

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
645
	
646
	int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
647
648
649
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
650
651
652
653
654
	      MatSetValue(mat_lagrange, 
			  index + counter, 
			  localDofMap.getMatIndex(k, it->first) + rStartInterior,
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
655
	    }
656
657
658
659
660
661
662
663
664
665
666
667
668
669
	    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);
670
671
672
673
674
675
676
	  }
	}
      }
    }

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

678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693

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

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

    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
694
695
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
696
697
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
698
    }
699
700
701
  }


702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
  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,
728
729
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
730
731
732
733
		 1, PETSC_NULL, 1, PETSC_NULL, 
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

734
    int rowCounter = rStartEdges;
735
736
737
738
739
740
741
742
743
744
745
746
747
    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");

748
	    int col = lagrangeMap.getMatIndex(i, **it);
749
	    double value = 1.0;
750
	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
751
752
	  }

753
	  rowCounter++;
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
	}
      }
    }

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


769
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
770
  {
771
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
772

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

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
806
807
808
809
810
811
      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
812

813
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
814
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
815
816
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
817
818
      KSPSetFromOptions(ksp_schur_primal);
    } else {
819
820
      MSG("Create direct schur primal solver!\n");

821
822
      TEST_EXIT(!augmentedLagrange)("Not yet supported!\n");

823
824
      double wtime = MPI::Wtime();

825
826
827
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

828
829
830
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
831

Thomas Witkowski's avatar
Thomas Witkowski committed
832
      Mat matBPi;
833
834
835
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
836
837
838
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

839

840
841
842
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
843

844
845
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

846
      for (int i = 0; i < nRowsRankB; i++) {
847
	PetscInt row = localDofMap.getStartDofs() + i;
848
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
849
850
851
852
853

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

854
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
855
856
      }

857
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
858
		primalDofMap.getLocalDofs())
859
	("Should not happen!\n");
860

861
862
863
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
864
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
865
866
867
868
869
870
871

 	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
872

Thomas Witkowski's avatar
Thomas Witkowski committed
873
	subdomain->solve(tmpVec, tmpVec);
874

Thomas Witkowski's avatar
Thomas Witkowski committed
875
	PetscScalar *tmpValues;
876
	VecGetArray(tmpVec, &tmpValues);
877
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
878
	  MatSetValue(matBPi, 
879
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
880
881
882
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
883
884
885
886
887
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
888
889
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
890

891
892
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
893
894

      Mat matPrimal;
895
896
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
897
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
898
899

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
900
      MatDestroy(&matBPi);
901

902
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
903
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
904
		      SAME_NONZERO_PATTERN);
905
906
907
908
909
910
      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);
911
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
912

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


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

937
    if (schurPrimalSolver == 0) {
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
      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);
      }
953
    }
954
955
956
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
957
958
959
  }


960
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
961
962
963
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

964
965
    // === Create FETI-DP solver object. ===

966
967
968
969
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
970
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
971
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
972
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
973

Thomas Witkowski's avatar
Thomas Witkowski committed
974
    if (stokesMode == false) {      
975
976
977
978
979
980
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
981
982
983
984
985
986
987
988
989
      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);	
      }
990
    }  else {
991
992
      TEST_EXIT_DBG(!augmentedLagrange)("Not yet implemented!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
993
994
995
996
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

997
998
999
1000
1001
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
1002
		     &fetiData, &mat_feti);
1003
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
1004
			   (void(*)(void))petscMultMatFetiInterface);      
1005
    }
1006

1007
    KSPCreate(mpiCommGlobal, &ksp_feti);
1008
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
1009
1010
1011
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
1012
    KSPSetFromOptions(ksp_feti);
1013
    if (stokesMode)