PetscSolverFeti.cc 77.8 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.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include "MatrixVector.h"
14
#include "parallel/PetscHelper.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
  PetscSolverFeti::PetscSolverFeti(string name)
31
    : PetscSolver(name),
32
33
34
35
36
37
      primalDofMap(COMPONENT_WISE),
      dualDofMap(COMPONENT_WISE),
      interfaceDofMap(COMPONENT_WISE),
      localDofMap(COMPONENT_WISE),
      lagrangeMap(COMPONENT_WISE),
      interiorDofMap(COMPONENT_WISE),
38
      schurPrimalSolver(0),
39
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
40
      subdomain(NULL),
41
      massMatrixSolver(NULL),
42
43
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
44
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
45
      printTimings(false),
Thomas Witkowski's avatar
d    
Thomas Witkowski committed
46
      dirichletMode(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
47
      stokesMode(false),
48
      augmentedLagrange(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
49
      pressureComponent(-1)
50
51
52
53
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

69
70
71
72
73
74
75
76
    preconditionerName = "";
    Parameters::get(initFileStr + "->right precon", preconditionerName);
    if (preconditionerName != "" && preconditionerName != "no") {
      ERROR_EXIT("FETI-DP does not support right preconditioning! (parameter \"%s->right precon\" has value \"%s\")\n",
		 initFileStr.c_str(), preconditionerName.c_str());
    }

    Parameters::get(initFileStr + "->feti->schur primal solver", schurPrimalSolver);
Thomas Witkowski's avatar
Thomas Witkowski committed
77
78
79
    TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
      ("Wrong solver \"%d\"for the Schur primal complement!\n", 
       schurPrimalSolver);
80

81
82
83
84
85
86
87
88
89
90
91
    Parameters::get(initFileStr + "->feti->stokes mode", stokesMode);
    if (stokesMode) {
      Parameters::get(initFileStr + "->feti->pressure component", pressureComponent);
      TEST_EXIT(pressureComponent >= 0)
	("FETI-DP in Stokes mode, no pressure component defined!\n");
    }
			   
    Parameters::get(initFileStr + "->feti->augmented lagrange", augmentedLagrange);

    Parameters::get(initFileStr + "->feti->symmetric", isSymmetric);

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

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


100
  void PetscSolverFeti::initialize()
101
  {
102
103
    FUNCNAME("PetscSolverFeti::initialize()");

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

107
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
Thomas Witkowski's avatar
Thomas Witkowski committed
108

Thomas Witkowski's avatar
Thomas Witkowski committed
109
    if (subdomain == NULL) {
110
      subdomain = new PetscSolverGlobalMatrix("");
111
      subdomain->setSymmetric(isSymmetric);
Thomas Witkowski's avatar
d    
Thomas Witkowski committed
112
113
114
115
      if (dirichletMode == 0)
	subdomain->setHandleDirichletRows(true);
      else
	subdomain->setHandleDirichletRows(false);
116

117
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
118
	subdomain->setMeshDistributor(meshDistributor, 
119
				      mpiCommGlobal, mpiCommLocal);
120
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
121
	subdomain->setMeshDistributor(meshDistributor, 
122
123
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
124
	subdomain->setLevel(meshLevel);
125
126
      }
    }
127

128
129
130
131
    primalDofMap.init(levelData, componentSpaces, feSpaces);   
    dualDofMap.init(levelData, componentSpaces, feSpaces, false);
    localDofMap.init(levelData, componentSpaces, feSpaces, meshLevel != 0);
    lagrangeMap.init(levelData, componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
132

Thomas Witkowski's avatar
Thomas Witkowski committed
133
    if (stokesMode)
134
      interfaceDofMap.init(levelData, componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
135

136
    if (fetiPreconditioner == FETI_DIRICHLET) {
137
138
139
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

140
      interiorDofMap.init(levelData, componentSpaces, feSpaces, false);
141
    }
142
143
144
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
145
146
147
148
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
d    
Thomas Witkowski committed
149
150
151
152
153
154
155
156
157
    if (dirichletMode == 1) {
      int nComponents = mat.getSize();
      for (int component = 0; component < nComponents; component++) {
	DOFMatrix* dofMat = mat[component][component];
	if (!dofMat)
	  continue;
	
	dirichletRows[component] = dofMat->getDirichletRows();
      }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
158
159
160
161
    }
  }


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

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

168
169
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

170
171
172
173
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
174
    if (fetiPreconditioner == FETI_DIRICHLET)
175
176
      interiorDofMap.clear();

177
178
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
179

180
181
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
182
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
183
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
184
    if (fetiPreconditioner == FETI_DIRICHLET)
185
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
186

187
188
189
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
190
191
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

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

198
199
200
201
202
203
    int nComponents = componentSpaces.size();
    for (int component = 0; component < nComponents; component++) {
      createPrimals(component);  
      createDuals(component);
      createInterfaceNodes(component);
      createIndexB(component);
204
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
205
206
207
208

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

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

215
216
217
    for (int component = 0; component < nComponents; component++) {
      createLagrange(component);
      createAugmentedLagrange(component);
218
    }
219

220
221
222
    lagrangeMap.update();


223
224
225
226
227
228
229
230
231
232
233
234
    // === ===

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

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

245
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
246
247
    for (unsigned int i = 0; i < componentSpaces.size(); i++) {
      const FiniteElemSpace *feSpace = componentSpaces[i];
248
249

      MSG("FETI-DP data for %d-ith component (FE space %p):\n", i, feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
250

251
      if (i == pressureComponent) {
Thomas Witkowski's avatar
Thomas Witkowski committed
252
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
253
254
	    interfaceDofMap[i].nRankDofs, 
	    interfaceDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
255
      } else {
256
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
257
258
259
	    primalDofMap[i].nRankDofs, 
	    primalDofMap[i].nLocalDofs,
	    primalDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
260
261
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
262
263
	    dualDofMap[i].nRankDofs, 
	    dualDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
266
267
	    lagrangeMap[i].nRankDofs, 
	    lagrangeMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
268
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    }
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
274
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
275
276

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


    bool writePrimals = false;
    Parameters::get("parallel->debug->write primals", writePrimals);
    if (writePrimals)
      PetscSolverFetiDebug::writePrimalFiles(*this);
288
289
290
  }


291
  void PetscSolverFeti::createPrimals(int component)
Thomas Witkowski's avatar
Thomas Witkowski committed
292
  {
293
    FUNCNAME("PetscSolverFeti::createPrimals()");  
294

295
    if (component == pressureComponent)
296
297
      return;

298
299
    const FiniteElemSpace *feSpace = componentSpaces[component];

300
301
302
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

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

    DofIndexSet primals;
308
    for (DofContainerSet::iterator it = vertices.begin(); 
309
	 it != vertices.end(); ++it) {
310
      
311
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
312
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
313

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
321
322
323
324
325
326
327
	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);
	}
      }
328
    }
329

330
331
332
333

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

334
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
335
      if (dofMap[feSpace].isRankDof(*it)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
336
	primalDofMap[component].insertRankDof(*it);
337
      } else
Thomas Witkowski's avatar
Thomas Witkowski committed
338
  	primalDofMap[component].insertNonRankDof(*it);
339
340
341
  }


342
  void PetscSolverFeti::createDuals(int component)
343
344
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
345

346
    if (component == pressureComponent)
347
348
      return;

349
350
    const FiniteElemSpace *feSpace = componentSpaces[component];

351
352
353
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
354
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
355
    
356
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
357
	 it != allBoundaryDofs.end(); ++it) {
358
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
359
360
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
361
      if (isPrimal(component, **it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
362
363
364
	continue;

      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
365
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
366
      } else {
367
	if (dofMapSd[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
368
	  dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
369
370
      }	  
    }
371
372
373
  }

  
374
  void PetscSolverFeti::createInterfaceNodes(int component)
375
376
377
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

378
    if (component != pressureComponent)
379
380
      return;

381
    const FiniteElemSpace *feSpace = componentSpaces[component];
382
383
384
385
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

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


398
  void PetscSolverFeti::createLagrange(int component)
399
400
401
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

402
    if (component == pressureComponent)
403
404
      return;

405
    const FiniteElemSpace *feSpace = componentSpaces[component];
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()) {
424
	  if (dofMapSd[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
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
444
	  if (!isPrimal(component, it.getDofIndex()))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
445
446
447
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
448

Thomas Witkowski's avatar
Thomas Witkowski committed
449
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
450
451
452
453
454
455
      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
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
461
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
462
463
	  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())
Thomas Witkowski's avatar
Thomas Witkowski committed
480
	if (!isPrimal(component, 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()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
491
	if (!isPrimal(component, it.getDofIndex())) {
492
493
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
494
	       dofMapSd[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())
Thomas Witkowski's avatar
Thomas Witkowski committed
511
	if (!isPrimal(component, it.getDofIndex()))
512
513
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
514
	       dofMapSd[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
515
516
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
517
	  else
Thomas Witkowski's avatar
Thomas Witkowski committed
518
	    lagrangeMap[component].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[component].getMap();
527
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
528
      if (dofMap[feSpace].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
529
	lagrangeMap[component].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 {
Thomas Witkowski's avatar
Thomas Witkowski committed
533
	lagrangeMap[component].insertNonRankDof(it->first);
534
535
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
536
    lagrangeMap[component].nRankDofs = nRankLagrange;
537
538
539
  }


540
  void PetscSolverFeti::createAugmentedLagrange(int component)
541
542
543
544
545
546
547
548
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


549
  void PetscSolverFeti::createIndexB(int component)
550
  {
551
    FUNCNAME("PetscSolverFeti::createIndexB()");
552

553
554

    const FiniteElemSpace *feSpace = componentSpaces[component];
555
    DOFAdmin* admin = feSpace->getAdmin();
556
557
558
559

    // === To ensure that all interior node on each rank are listen first in ===
    // === the global index of all B nodes, insert all interior nodes first, ===
    // === without defining a correct index.                                 ===
560

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
570
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
571
	localDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
572
573
	
	if (fetiPreconditioner == FETI_DIRICHLET)
574
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
575
576
577
	
	nLocalInterior++;	
      } else {
578
	if (dofMapSd[feSpace].isRankDof(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
579
	  localDofMap[component].insertRankDof(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
580
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
581
	  localDofMap[component].insertNonRankDof(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
582
583
584
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
585
      }
586
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
587
    
588
589
    // === And finally, add the global indicies of all dual nodes. ===

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


604
  void PetscSolverFeti::createMatLagrange()
605
606
607
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

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

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

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

620
621
622
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

623
624
625
626
627
628
629
    // === 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.                                                     ===

630
    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
631
      DofMap &dualMap = dualDofMap[component].getMap();
632

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

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

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

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
665
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
666
667
668
669
670
671
	  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);
672
673
674
675
676
677
678
	  }
	}
      }
    }

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

680

Thomas Witkowski's avatar
So on    
Thomas Witkowski committed
681
682
683
684
685
686
687
688
689
690
691
692
    int nZeroRows = PetscSolverFetiDebug::testZeroRows(mat_lagrange);
    int m,n;
    MatGetSize(mat_lagrange, &m ,&n);
    MSG("Lagrange matrix has %d zero rows and global size of %d %d!\n", nZeroRows, m, n);

    PetscViewer petscView;
    PetscViewerBinaryOpen(PETSC_COMM_WORLD, "lagrange.mat", 
			  FILE_MODE_WRITE, &petscView);
    MatView(mat_lagrange, petscView);
    PetscViewerDestroy(&petscView);


693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
    // === 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
708
709
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
710
711
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
712
    }
713
714
  }

715

Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
  bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
720
    if (meshDistributor->getIntBoundary().getDegreeOwn(edge) != 3)
      return false;

Thomas Witkowski's avatar
Thomas Witkowski committed
721
722
723
724
725
726
727
728
729
730
731
732
733
734
    Element *el = edge.el;
    int i0 = el->getVertexOfEdge(edge.ithObj, 0);
    int i1 = el->getVertexOfEdge(edge.ithObj, 1);
    DegreeOfFreedom d0 = el->getDof(i0, 0);
    DegreeOfFreedom d1 = el->getDof(i1, 0);
    WorldVector<double> c0, c1;
    el->getMesh()->getDofIndexCoords(d0, feSpace, c0);
    el->getMesh()->getDofIndexCoords(d1, feSpace, c1);
    bool xe = fabs(c0[0] - c1[0]) < 1e-8;
    bool ye = fabs(c0[1] - c1[1]) < 1e-8;
    bool ze = fabs(c0[2] - c1[2]) < 1e-8;
    int counter = static_cast<int>(xe) + static_cast<int>(ye) + static_cast<int>(ze);
    return (counter == 2);
  }
735

736

737
  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseEdges()
738
  {
739
    FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");
740
741

    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
Thomas Witkowski's avatar
Thomas Witkowski committed
742
    std::set<BoundaryObject> allEdges;
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
      if (it->rankObj.subObj == EDGE && 
	  testWirebasketEdge(it->rankObj, feSpaces[0]) &&
	  allEdges.count(it->rankObj) == 0)
	allEdges.insert(it->rankObj);

    int nEmptyEdges = 0;
    vector<vector<BoundaryObject> > data;
    for (std::set<BoundaryObject>::iterator it = allEdges.begin();
	 it != allEdges.end(); ++it) {
      DofContainer edgeDofs;
      it->el->getAllDofs(feSpaces[0], *it, edgeDofs);
      if (edgeDofs.size() == 0) {
	nEmptyEdges++;
      } else {
	vector<BoundaryObject> oneBoundary;
	oneBoundary.push_back(*it);      
	data.push_back(oneBoundary);
      }
    }
763

764
765
766
    int nEdges = allEdges.size();
    mpi::globalAdd(nEdges);
    mpi::globalAdd(nEmptyEdges);
Thomas Witkowski's avatar
Thomas Witkowski committed
767

768
769
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
    MSG("Coarse space edges: %d (empty: %d)\n", nEdges, nEmptyEdges);
    
    return data;
  }


  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseFaces()
  {
    FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");

    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
    map<int, std::set<BoundaryObject> > allFaces;
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
      if (it->rankObj.subObj == FACE)
	allFaces[it.getRank()].insert(it->rankObj);

    int nEmptyFaces = 0;
    vector<vector<BoundaryObject> > data;
    for (map<int, std::set<BoundaryObject> >::iterator it = allFaces.begin();
	 it != allFaces.end(); ++it) {
      vector<BoundaryObject> oneBoundary;
      for (std::set<BoundaryObject>::iterator bIt = it->second.begin();
	   bIt != it->second.end(); ++bIt) {
	DofContainer faceDofs;
	bIt->el->getAllDofs(feSpaces[0], *bIt, faceDofs);
	if (faceDofs.size())
	  oneBoundary.push_back(*bIt);
795
      }
796
797
798
799
      if (oneBoundary.size()) 
	data.push_back(oneBoundary);
      else 
	nEmptyFaces++;
800
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
801

802
    int nFaces = allFaces.size();
803
    mpi::globalAdd(nFaces);
804
    mpi::globalAdd(nEmptyFaces);
805

806
807
808
    MSG("Coarse space faces: %d (empty: %d)\n", nFaces, nEmptyFaces);
   
    return data;
809
810
811
812
813
814
815
816
817
818
819
820
  }


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

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

821
822
823
824
825
    vector<vector<BoundaryObject> > allEdges = getCoarseEdges();
    vector<vector<BoundaryObject> > allFaces = getCoarseFaces(); 
    allEdges.insert(allEdges.end(), allFaces.begin(), allFaces.end());

    nRankEdges = allEdges.size();
826
827
828
829
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

    MSG("nRankEdges = %d, nOverallEdges = %d\n", nRankEdges, nOverallEdges);
Thomas Witkowski's avatar
Thomas Witkowski committed
830

831
832
833
    nRankEdges *= componentSpaces.size();
    rStartEdges *= componentSpaces.size();
    nOverallEdges *= componentSpaces.size();
834
835

    MatCreateAIJ(mpiCommGlobal,
836
837
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
838
		 2, PETSC_NULL, 2, PETSC_NULL, 
839
840
841
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

842
    int rowCounter = rStartEdges;
843
844
    for (vector<vector<BoundaryObject> >::iterator it = allEdges.begin(); 
	 it != allEdges.end(); ++it) {
845
      for (int component = 0; component < componentSpaces.size(); component++) {
846
847
848
849
850
	for (vector<BoundaryObject>::iterator edgeIt = it->begin();
	     edgeIt != it->end(); ++edgeIt) {
	  
	  DofContainer edgeDofs;
	  edgeIt->el->getAllDofs(componentSpaces[component], *edgeIt, edgeDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
851
	  
852
	  TEST_EXIT(edgeDofs.size())("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
853
	  
854
855
856
857
858
859
860
861
862
	  for (DofContainer::iterator it = edgeDofs.begin();
	       it != edgeDofs.end(); ++it) {
	    TEST_EXIT(isPrimal(component, **it) == false)
	      ("Should not be primal!\n");
	    
	    int col = lagrangeMap.getMatIndex(component, **it);
	    double value = 1.0;
	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
	  }
863
	}
864
865
	rowCounter++;	
      }
866
867
868
869
870
    }

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

Thomas Witkowski's avatar
So on    
Thomas Witkowski committed
871
872
873
874
875
    int nZeroRows = PetscSolverFetiDebug::testZeroRows(mat_augmented_lagrange);
    int m,n;
    MatGetSize(mat_augmented_lagrange, &m ,&n);
    MSG("Augmented lagrange matrix has %d zero rows and global size of %d %d!\n", nZeroRows, m, n);

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


884
  void PetscSolverFeti::createSchurPrimalKsp()
885
  {
886
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
887

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

891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
      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;
908
	schurPrimalAugmentedData.nestedVec = true;
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927

	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
928

929
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
930
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
931
932
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
933
934
      KSPSetFromOptions(ksp_schur_primal);
    } else {
935
936
      MSG("Create direct schur primal solver!\n");

937
938
      double wtime = MPI::Wtime();

939
940
941
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

942

943
      // === Create explicit matrix representation of the Schur primal system. ===
944

945
946
947
948
      if (!augmentedLagrange)
	createMatExplicitSchurPrimal();
      else
	createMatExplicitAugmentedSchurPrimal();
949
950


Thomas Witkowski's avatar
Thomas Witkowski committed
951
      // === Create KSP solver object and set appropriate solver options. ===
952

953
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
<