PetscSolverFeti.cc 68.9 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
257
258
259
260
261
262
263
264
265
266
267
268
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
	MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
	    primalDofMap[feSpace].nRankDofs, 
	    primalDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
	    dualDofMap[feSpace].nRankDofs, 
	    dualDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
	    lagrangeMap[feSpace].nRankDofs, 
	    lagrangeMap[feSpace].nOverallDofs);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
269
270
271
// 	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
272
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
273
    }
274

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

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


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

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

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

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

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

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

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

325
326
327
328

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

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


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

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

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

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

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

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

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

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


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

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

397
398
    boundaryDofRanks[feSpace].clear();

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

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

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

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

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

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

462
463
464
465

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

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

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

    stdMpi.updateSendDataSize();

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

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

496
497
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
512

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

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


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

    if (!augmentedLagrange)
      return;
  }


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

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

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

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

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

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


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

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

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

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

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

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
633
634
635
636
	
	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
637
638
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
639

640
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
641
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
642
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
643
	    index++;	      
644
645
646
647
648
649
650
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
652
653
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
654
655
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
656
    }
657
658
659
  }


660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
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
722
723
724
725
726
  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,
		 lagrangeMap.getRankDofs(), nRankEdges,
		 lagrangeMap.getOverallDofs(), nOverallEdges,
		 1, PETSC_NULL, 1, PETSC_NULL, 
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

    int colCounter = rStartEdges;
    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");

	    int row = lagrangeMap.getMatIndex(i, **it);
	    double value = 1.0;
	    MatSetValue(mat_augmented_lagrange, row, colCounter, value, INSERT_VALUES);
	  }

	  colCounter++;
	}
      }
    }

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


727
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
728
  {
729
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
730

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

734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
      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
770

771
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
772
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
773
774
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
      KSPSetFromOptions(ksp_schur_primal);
    } else {
777
778
      MSG("Create direct schur primal solver!\n");

779
780
      TEST_EXIT(!augmentedLagrange)("Not yet supported!\n");

781
782
      double wtime = MPI::Wtime();

783
784
785
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

786
787
788
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
789

Thomas Witkowski's avatar
Thomas Witkowski committed
790
      Mat matBPi;
791
792
793
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
794
795
796
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

797

798
799
800
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
801

802
803
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

804
      for (int i = 0; i < nRowsRankB; i++) {
805
	PetscInt row = localDofMap.getStartDofs() + i;
806
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
807
808
809
810
811

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

812
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
813
814
      }

815
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
816
		primalDofMap.getLocalDofs())
817
	("Should not happen!\n");
818

819
820
821
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
822
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
823
824
825
826
827
828
829

 	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
830

Thomas Witkowski's avatar
Thomas Witkowski committed
831
	subdomain->solve(tmpVec, tmpVec);
832

Thomas Witkowski's avatar
Thomas Witkowski committed
833
	PetscScalar *tmpValues;
834
	VecGetArray(tmpVec, &tmpValues);
835
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
836
	  MatSetValue(matBPi, 
837
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
838
839
840
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
841
842
843
844
845
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
846
847
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
848

849
850
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
851
852

      Mat matPrimal;
853
854
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
855
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
856
857

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
858
      MatDestroy(&matBPi);
859

860
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
861
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
862
		      SAME_NONZERO_PATTERN);
863
864
865
866
867
868
      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);
869
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
870

Thomas Witkowski's avatar
Thomas Witkowski committed
871
      if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
872
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
873
874
875
876
877
878
879
880
881
882
	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
883
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
884
885
886
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
887
    }
888
889
890
891
892
893
894
  }


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

895
    if (schurPrimalSolver == 0) {
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
      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);
      }
911
    }
912
913
914
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
915
916
917
  }


918
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
919
920
921
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

922
923
    // === Create FETI-DP solver object. ===

924
925
926
927
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
928
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
929
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
930
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
931

Thomas Witkowski's avatar
Thomas Witkowski committed
932
    if (stokesMode == false) {      
933
934
935
936
937
938
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
939
940
941
942
943
944
945
946
947
      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);	
      }
948
    }  else {
949
950
      TEST_EXIT_DBG(!augmentedLagrange)("Not yet implemented!\n");

Thomas Witkowski's avatar