PetscSolverFeti.cc 75.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
Thomas Witkowski committed
46
      stokesMode(false),
47
      augmentedLagrange(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
48
      pressureComponent(-1)
49
50
51
52
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

68
69
70
71
72
73
74
75
    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
76
77
78
    TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
      ("Wrong solver \"%d\"for the Schur primal complement!\n", 
       schurPrimalSolver);
79

80
81
82
83
84
85
86
87
88
89
90
    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);

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

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


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
108
    if (subdomain == NULL) {
109
      subdomain = new PetscSolverGlobalMatrix("");
110
      subdomain->setSymmetric(isSymmetric);
111
      subdomain->setHandleDirichletRows(false);
112

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

124
125
126
127
    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
128

Thomas Witkowski's avatar
Thomas Witkowski committed
129
    if (stokesMode)
130
      interfaceDofMap.init(levelData, componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
131

132
    if (fetiPreconditioner == FETI_DIRICHLET) {
133
134
135
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

136
      interiorDofMap.init(levelData, componentSpaces, feSpaces, false);
137
    }
138
139
140
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
141
142
143
144
145
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

    int nComponents = mat.getSize();
146
147
    for (int component = 0; component < nComponents; component++) {
      DOFMatrix* dofMat = mat[component][component];
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
148
149
      if (!dofMat)
	continue;
150
151

      dirichletRows[component] = dofMat->getDirichletRows();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
152
153
154
155
    }
  }


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

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

162
163
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

164
165
166
167
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
168
    if (fetiPreconditioner == FETI_DIRICHLET)
169
170
      interiorDofMap.clear();

171
172
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
173

174
175
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
176
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
177
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
178
    if (fetiPreconditioner == FETI_DIRICHLET)
179
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
180

181
182
183
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
184
185
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

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

192
193
194
195
196
197
    int nComponents = componentSpaces.size();
    for (int component = 0; component < nComponents; component++) {
      createPrimals(component);  
      createDuals(component);
      createInterfaceNodes(component);
      createIndexB(component);
198
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
201
202

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
203
    if (fetiPreconditioner == FETI_DIRICHLET)
204
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
205

Thomas Witkowski's avatar
Thomas Witkowski committed
206
    if (stokesMode)
207
208
      interfaceDofMap.update();

209
210
211
    for (int component = 0; component < nComponents; component++) {
      createLagrange(component);
      createAugmentedLagrange(component);
212
    }
213

214
215
216
    lagrangeMap.update();


217
218
219
220
221
222
223
224
225
226
227
228
    // === ===

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

229
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
230
231
232
233
234
235
236
237
238
			   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);
    }

239
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
240
241
    for (unsigned int i = 0; i < componentSpaces.size(); i++) {
      const FiniteElemSpace *feSpace = componentSpaces[i];
242
243

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
265
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
268
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
269
270

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
271
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
272
      timeCounter = MPI::Wtime() - timeCounter;
273
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
	  timeCounter);
    }
276
277
278
  }


279
  void PetscSolverFeti::createPrimals(int component)
Thomas Witkowski's avatar
Thomas Witkowski committed
280
  {
281
    FUNCNAME("PetscSolverFeti::createPrimals()");  
282

283
    if (component == pressureComponent)
284
285
      return;

286
287
    const FiniteElemSpace *feSpace = componentSpaces[component];

288
289
290
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
291
    // Set of DOF indices that are considered to be primal variables.
292
    DofContainerSet& vertices = 
293
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
294
295

    DofIndexSet primals;
296
    for (DofContainerSet::iterator it = vertices.begin(); 
297
	 it != vertices.end(); ++it) {
298
      
299
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
300
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
301

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
302
303
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
304
      } else {
305
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
306
307
308
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
309
310
311
312
313
314
315
	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);
	}
      }
316
    }
317

318
319
320
321

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

322
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
323
      if (dofMap[feSpace].isRankDof(*it)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
324
	primalDofMap[component].insertRankDof(*it);
325
      } else
Thomas Witkowski's avatar
Thomas Witkowski committed
326
  	primalDofMap[component].insertNonRankDof(*it);
327
328
329
  }


330
  void PetscSolverFeti::createDuals(int component)
331
332
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
333

334
    if (component == pressureComponent)
335
336
      return;

337
338
    const FiniteElemSpace *feSpace = componentSpaces[component];

339
340
341
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
342
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
343
    
344
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
345
	 it != allBoundaryDofs.end(); ++it) {
346
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
347
348
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
349
      if (isPrimal(component, **it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
350
351
352
	continue;

      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
353
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
354
      } else {
355
	if (dofMapSd[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
356
	  dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
357
358
      }	  
    }
359
360
361
  }

  
362
  void PetscSolverFeti::createInterfaceNodes(int component)
363
364
365
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

366
    if (component != pressureComponent)
367
368
      return;

369
    const FiniteElemSpace *feSpace = componentSpaces[component];
370
371
372
373
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
374
	 it != allBoundaryDofs.end(); ++it) {
375
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
376
377
	continue;      
      
378
      if (dofMap[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
379
	interfaceDofMap[component].insertRankDof(**it);
380
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
381
	interfaceDofMap[component].insertNonRankDof(**it);      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
382
    }
383
384
385
  }


386
  void PetscSolverFeti::createLagrange(int component)
387
388
389
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

390
    if (component == pressureComponent)
391
392
      return;

393
    const FiniteElemSpace *feSpace = componentSpaces[component];
394
395
    boundaryDofRanks[feSpace].clear();

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

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

      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()) {
412
	  if (dofMapSd[feSpace].isRankDof(it.getDofIndex()))
413
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
414
	  else
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
	    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
430
431
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
432
	  if (!isPrimal(component, it.getDofIndex()))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
433
434
435
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
436

Thomas Witkowski's avatar
Thomas Witkowski committed
437
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
438
439
440
441
442
443
      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)).         ===

444
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
445
446
447
448
    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
449
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
450
451
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

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

459
460
461
462

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

463
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
464

465
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
466
467
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
468
	if (!isPrimal(component, it.getDofIndex()))
469
470
471
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
472
473
474

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
490
	stdMpi.recv(it.getRank());
491
    }
492

493
494
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
509

510
511
512
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


528
  void PetscSolverFeti::createAugmentedLagrange(int component)
529
530
531
532
533
534
535
536
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


537
  void PetscSolverFeti::createIndexB(int component)
538
  {
539
    FUNCNAME("PetscSolverFeti::createIndexB()");
540

541
542

    const FiniteElemSpace *feSpace = componentSpaces[component];
543
    DOFAdmin* admin = feSpace->getAdmin();
544
545
546
547

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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
558
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
559
	localDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
560
561
	
	if (fetiPreconditioner == FETI_DIRICHLET)
562
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
563
564
565
	
	nLocalInterior++;	
      } else {
566
	if (dofMapSd[feSpace].isRankDof(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
567
	  localDofMap[component].insertRankDof(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
568
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
569
	  localDofMap[component].insertNonRankDof(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
570
571
572
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
573
      }
574
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
575
    
576
577
    // === And finally, add the global indicies of all dual nodes. ===

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


592
  void PetscSolverFeti::createMatLagrange()
593
594
595
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

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

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

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

608
609
610
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

611
612
613
614
615
616
617
    // === 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.                                                     ===

618
    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
619
      DofMap &dualMap = dualDofMap[component].getMap();
620

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

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

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

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
653
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
654
655
656
657
658
659
	  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);
660
661
662
663
664
665
666
	  }
	}
      }
    }

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

668

Thomas Witkowski's avatar
So on  
Thomas Witkowski committed
669
670
671
672
673
674
675
676
677
678
679
680
    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);


681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
    // === 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
696
697
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
698
699
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
700
    }
701
702
  }

703

Thomas Witkowski's avatar
Thomas Witkowski committed
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
  bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
  {
    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);
  }
720

721

722
  void PetscSolverFeti::createMatAugmentedLagrange()
723
724
725
726
727
728
729
730
731
732
  {
    FUNCNAME("PetscSolverFeti::createMatAugmentedLagrange()");

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

    nOverallEdges = 0;
    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
Thomas Witkowski's avatar
Thomas Witkowski committed
733
    std::set<BoundaryObject> allEdges;
734
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
735
736
737
      if ((it->rankObj.subObj == FACE ||
	   (it->rankObj.subObj == EDGE && 
	    testWirebasketEdge(it->rankObj, feSpaces[0]))) && 
738
739
740
741
742
	  allEdges.count(it->rankObj) == 0) {
	bool dirichletOnlyEdge = true;

	DofContainer edgeDofs;
	it->rankObj.el->getAllDofs(feSpaces[0], it->rankObj, edgeDofs);
743
	
744
745
	for (DofContainer::iterator dit = edgeDofs.begin();
	     dit != edgeDofs.end(); ++dit) {
746
	  if (dirichletRows[0].count(**dit) == 0) {
747
748
749
750
751
752
753
754
755
	    dirichletOnlyEdge = false;
	    break;
	  }
	}

	if (!dirichletOnlyEdge)
	  allEdges.insert(it->rankObj);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
756
757

    nRankEdges = allEdges.size();    
758
759
760
761
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

    MSG("nRankEdges = %d, nOverallEdges = %d\n", nRankEdges, nOverallEdges);
762
    
763
764
765
    nRankEdges *= componentSpaces.size();
    rStartEdges *= componentSpaces.size();
    nOverallEdges *= componentSpaces.size();
766
767

    MatCreateAIJ(mpiCommGlobal,
768
769
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
770
		 2, PETSC_NULL, 2, PETSC_NULL, 
771
772
773
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

774
    int rowCounter = rStartEdges;
Thomas Witkowski's avatar
Thomas Witkowski committed
775
776
    for (std::set<BoundaryObject>::iterator edgeIt = allEdges.begin(); 
	 edgeIt != allEdges.end(); ++edgeIt) {
777
      for (int component = 0; component < componentSpaces.size(); component++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
778
	DofContainer edgeDofs;
779
	edgeIt->el->getAllDofs(componentSpaces[component], *edgeIt, edgeDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
780
781
782
	
	for (DofContainer::iterator it = edgeDofs.begin();
	     it != edgeDofs.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
783
	  TEST_EXIT_DBG(isPrimal(component, **it) == false)
Thomas Witkowski's avatar
Thomas Witkowski committed
784
785
	    ("Should not be primal!\n");
	  
786
	  if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
787
788
	    continue;
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
789
	  int col = lagrangeMap.getMatIndex(component, **it);
Thomas Witkowski's avatar
Thomas Witkowski committed
790
791
	  double value = 1.0;
	  MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
792
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
793
794
795
	
	rowCounter++;
      }      
796
797
798
799
800
    }

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

Thomas Witkowski's avatar
So on  
Thomas Witkowski committed
801
802
803
804
805
    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);

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


814
  void PetscSolverFeti::createSchurPrimalKsp()
815
  {
816
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
817

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

821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
      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;
838
	schurPrimalAugmentedData.nestedVec = true;
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857

	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
858

Thomas Witkowski's avatar