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


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

namespace AMDiS {

  using namespace std;

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
  PetscErrorCode KSPMonitorFeti(KSP ksp, PetscInt n, PetscReal rnorm, void *dummy)
  {    
    Vec            Br,v,w;
    Mat            A;

    KSPGetOperators(ksp, &A, PETSC_NULL, PETSC_NULL);
    MatGetVecs(A, &w, &v);
    KSPBuildResidual(ksp, v, w, &Br);

    Vec nest0, nest1;
    VecNestGetSubVec(Br, 0, &nest0);
    VecNestGetSubVec(Br, 1, &nest1);

    PetscScalar norm0, norm1;
    VecNorm(nest0, NORM_2, &norm0);
    VecNorm(nest1, NORM_2, &norm1);

    VecDestroy(&v);
    VecDestroy(&w);

    PetscPrintf(PETSC_COMM_WORLD, "%3D KSP Component residual norm [%1.12e %1.12e] and relative [%1.12e]\n",
		norm0, norm1, rnorm);

    return 0;
  }

54
55
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
56
      schurPrimalSolver(0),
57
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
58
      subdomain(NULL),
59
60
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
61
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
62
      printTimings(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
63
64
      stokesMode(false),
      pressureComponent(-1)
65
66
67
68
69
70
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
71
      MSG("Create FETI-DP solver with no preconditioner!\n");
72
73
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
74
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
75
76
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
77
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
78
79
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
80
81
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
82
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
83
84
85
86
87

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

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

    Parameters::get("parallel->print timings", printTimings);
94
95
96
  }


97
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
98
  {
99
100
    FUNCNAME("PetscSolverFeti::initialize()");

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

104
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
105
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
106

Thomas Witkowski's avatar
Thomas Witkowski committed
107

Thomas Witkowski's avatar
Thomas Witkowski committed
108
109
110
111
112
113
114
    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");

Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
      MSG("========= !!!! SUPER WARNING !!! ========\n");
      MSG("Use make use of stokes mode which is work in progress! We guarantee for nothing and even less than this!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
117
      pressureFeSpace = feSpaces[pressureComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
118
119
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
122

123
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
124
	subdomain->setMeshDistributor(meshDistributor, 
125
				      mpiCommGlobal, mpiCommLocal);
126
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
127
	subdomain->setMeshDistributor(meshDistributor, 
128
129
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
130
	subdomain->setLevel(meshLevel);
131
132
      }
    }
133

134
135
136
137
    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
138

Thomas Witkowski's avatar
Thomas Witkowski committed
139
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
140
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
141

142
    if (fetiPreconditioner == FETI_DIRICHLET) {
143
144
145
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

146
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
147
    }
148
149
150
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

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


172
173
174
175
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
176
177
    double timeCounter = MPI::Wtime();

178
179
180
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
181

182
183
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

184
185
186
187
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
188
    if (fetiPreconditioner == FETI_DIRICHLET)
189
190
      interiorDofMap.clear();

191
192
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
193

194
195
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
196
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
197
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
198
    if (fetiPreconditioner == FETI_DIRICHLET)
199
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
200

201
202
203
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
204
205
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
206
    if (stokesMode) {
207
208
209
210
211
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

212
213
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
214
    
215
216
      createPrimals(feSpace);  

217
218
219
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
220

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
221
      createIndexB(feSpace);     
222
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
223
224
225
226

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
227
    if (fetiPreconditioner == FETI_DIRICHLET)
228
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
229

Thomas Witkowski's avatar
Thomas Witkowski committed
230
    if (stokesMode)
231
232
      interfaceDofMap.update();

233
234
235
236
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
237

238
239
240
    lagrangeMap.update();


241
242
243
244
245
246
247
248
249
250
251
252
    // === ===

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

253
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
254
255
256
257
258
259
260
261
262
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
268
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
	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
285
286
287
// 	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
288
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
289
    }
290

Thomas Witkowski's avatar
Thomas Witkowski committed
291
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
292
293
294
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
295
296

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
297
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
298
      timeCounter = MPI::Wtime() - timeCounter;
299
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
	  timeCounter);
    }
302
303
304
  }


305
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
306
  {
307
    FUNCNAME("PetscSolverFeti::createPrimals()");  
308

Thomas Witkowski's avatar
Thomas Witkowski committed
309
    if (feSpace == pressureFeSpace)
310
311
      return;

312
313
314
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
315
    // Set of DOF indices that are considered to be primal variables.
316
    DofContainerSet& vertices = 
317
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
318
319

    DofIndexSet primals;
320
    for (DofContainerSet::iterator it = vertices.begin(); 
321
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
322
323
      if (dirichletRows[feSpace].count(**it))
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
324

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
325
326
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
327
      } else {
328
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
331
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
332
333
334
335
336
337
338
	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);
	}
      }
339
    }
340

341
342
343
344

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

345
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
346
347
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
348
      else
349
  	primalDofMap[feSpace].insertNonRankDof(*it);
350
351
352
  }


353
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
354
355
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
356

Thomas Witkowski's avatar
Thomas Witkowski committed
357
    if (feSpace == pressureFeSpace)
358
359
      return;

360
361
362
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
363
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
364
    
365
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
366
367
368
369
370
371
372
373
374
375
376
	 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))
377
	  dualDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
378
379
      }	  
    }
380
381
382
  }

  
383
384
385
386
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
387
    if (feSpace != pressureFeSpace)
388
389
390
391
392
393
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
394
395
396
397
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;      
      
398
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
399
	interfaceDofMap[feSpace].insertRankDof(**it);
400
      else
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
401
402
	interfaceDofMap[feSpace].insertNonRankDof(**it);      
    }
403
404
405
  }


406
407
408
409
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
410
    if (feSpace == pressureFeSpace)
411
412
      return;

413
414
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
    // 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
418
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
419

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
420
    if (meshLevel > 0) {
421
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
422
423
424
425
426
427
428
429
430

      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
431
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
432
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
433
	  else
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
	    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
449
450
451
452
453
454
	   !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());
    }
455

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
456
457
458
459
460
461
462
    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)).         ===

463
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
464
465
466
467
468
469
470
    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);

471
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
472
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
473
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
474
475
476
477
	}
      }
    }

478
479
480
481

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

482
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
483

484
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
485
486
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
487
	if (!isPrimal(feSpace, it.getDofIndex()))
488
489
490
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
491
492
493

    stdMpi.updateSendDataSize();

494
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
495
	 !it.end(); it.nextRank()) {
496
      bool recvFromRank = false;
497
      for (; !it.endDofIter(); it.nextDof()) {
498
	if (!isPrimal(feSpace, it.getDofIndex())) {
499
500
501
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
502
503
504
	    recvFromRank = true;
	    break;
	  }
505
	}
506
      }
507
508

      if (recvFromRank)
509
	stdMpi.recv(it.getRank());
510
    }
511

512
513
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
528

529
530
531
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


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

551
    DOFAdmin* admin = feSpace->getAdmin();
552
553
554
555

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

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

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

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


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

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

607
608
    // === Create distributed matrix for Lagrange constraints. ===

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

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

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
640
641
642
643
	
	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
644
645
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
646

647
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
648
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
649
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
650
	    index++;	      
651
652
653
654
655
656
657
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
659
660
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
661
662
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
663
    }
664
665
666
  }


667
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
668
  {
669
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
670

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

Thomas Witkowski's avatar
Thomas Witkowski committed
674
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
675
      
676
677
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
678

679
      MatCreateShell(mpiCommGlobal,
680
681
682
683
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
684
685
686
687
688
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
689
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
690
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
691
692
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
693
694
      KSPSetFromOptions(ksp_schur_primal);
    } else {
695
696
      MSG("Create direct schur primal solver!\n");

697
698
      double wtime = MPI::Wtime();

699
700
701
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

702
703
704
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
705

Thomas Witkowski's avatar
Thomas Witkowski committed
706
      Mat matBPi;
707
708
709
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
710
711
712
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

713

714
715
716
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
717

718
719
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

720
      for (int i = 0; i < nRowsRankB; i++) {
721
	PetscInt row = localDofMap.getStartDofs() + i;
722
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
723
724
725
726
727

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

728
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
729
730
      }

731
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
732
		primalDofMap.getLocalDofs())
733
	("Should not happen!\n");
734

735
736
737
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
738
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
739
740
741
742
743
744
745

 	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
746

Thomas Witkowski's avatar
Thomas Witkowski committed
747
	subdomain->solve(tmpVec, tmpVec);
748

Thomas Witkowski's avatar
Thomas Witkowski committed
749
	PetscScalar *tmpValues;
750
	VecGetArray(tmpVec, &tmpValues);
751
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
752
	  MatSetValue(matBPi, 
753
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
754
755
756
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
757
758
759
760
761
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
762
763
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
764

765
766
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
767
768

      Mat matPrimal;
769
770
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
771
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
772
773

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
774
      MatDestroy(&matBPi);
775

776
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
777
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
778
		      SAME_NONZERO_PATTERN);
779
780
781
782
783
784
      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);
785
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
786

Thomas Witkowski's avatar
Thomas Witkowski committed
787
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
788
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
789
790
791
792
793
794
795
796
797
798
	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
799
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
800
801
802
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
803
    }
804
805
806
807
808
809
810
  }


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

811
    if (schurPrimalSolver == 0) {
812
      schurPrimalData.subSolver = NULL;
813

814
815
816
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
817
818
819
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
820
821
822
  }


823
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
824
825
826
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

827
828
    // === Create FETI-DP solver object. ===

829
830
831
832
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
833
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
834
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
835
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
836

Thomas Witkowski's avatar
Thomas Witkowski committed
837
    if (stokesMode == false) {      
838
839
840
841
842
843
844
845
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
846
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
847
848
849
850
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

851
852
853
854
855
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
856
		     &fetiData, &mat_feti);
857
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
858
			   (void(*)(void))petscMultMatFetiInterface);      
859
    }
860

861
    KSPCreate(mpiCommGlobal, &ksp_feti);
862
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
863
864
865
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
866
    KSPSetFromOptions(ksp_feti);
867
868


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
869
870


871
    // === Create FETI-DP preconditioner object. ===
872

Thomas Witkowski's avatar
Thomas Witkowski committed
873
    if (fetiPreconditioner != FETI_NONE || stokesMode) {
874
875
876
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
877

878

Thomas Witkowski's avatar
Thomas Witkowski committed
879
880
881
882
883
    // === Create null space objects. ===
    
    createNullSpace();


884
    switch (fetiPreconditioner) {
885
886
    case FETI_DIRICHLET:
      TEST_EXIT(meshLevel == 0)
887
888
	("Check for localDofMap.getLocalMatIndex, which should not work for multilevel FETI-DP!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
889
      TEST_EXIT(!stokesMode)
890
891
	("Stokes mode does not yet support the Dirichlet precondition!\n");

892
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
893
894
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
895
896
897
898
899
900
      KSPSetOptionsPrefix(ksp_interior, "precon_interior_");
      KSPSetType(ksp_interior, KSPPREONLY);
      PC pc_interior;
      KSPGetPC(ksp_interior, &pc_interior);
      PCSetType(pc_interior, PCLU);
      PCFactorSetMatSolverPackage(pc_interior, MATSOLVERUMFPACK);
901
902
903
904
905
906
907
908
909
      KSPSetFromOptions(ksp_interior);
            
      fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior;
      fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals;
      fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals;
      fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
      fetiDirichletPreconData.ksp_interior = &ksp_interior;
      
910
911
      localDofMap.createVec(fetiDirichletPreconData.tmp_vec_b, 
			    nGlobalOverallInterior);
912
913
914
915
916
917
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals1));
      MatGetVecs(mat_interior_interior, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_interior));
918

919
920
921
      TEST_EXIT_DBG(meshLevel == 0)
	("Should not happen, check usage of localDofMap!\n");

922
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
923
	DofMap &dualMap = dualDofMap[feSpaces[i]].getMap();
924
	for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
925
	  DegreeOfFreedom d = it->first;	  
926
927
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
928