PetscSolverFeti.cc 79 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
343
  }


344
  void PetscSolverFeti::createDuals(int component)
345
346
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
347

348
    if (component == pressureComponent)
349
350
      return;

351
352
    const FiniteElemSpace *feSpace = componentSpaces[component];

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

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

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

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

  
376
  void PetscSolverFeti::createInterfaceNodes(int component)
377
378
379
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

380
    if (component != pressureComponent)
381
382
      return;

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

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


400
  void PetscSolverFeti::createLagrange(int component)
401
402
403
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

404
    if (component == pressureComponent)
405
406
      return;

407
    const FiniteElemSpace *feSpace = componentSpaces[component];
408
409
    boundaryDofRanks[feSpace].clear();

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

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

      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()) {
426
	  if (dofMapSd[feSpace].isRankDof(it.getDofIndex()))
427
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
428
	  else
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
	    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
444
445
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
446
	  if (!isPrimal(component, it.getDofIndex()))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
447
448
449
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
450

Thomas Witkowski's avatar
Thomas Witkowski committed
451
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
452
453
454
455
456
457
      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)).         ===

458
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
459
460
461
462
    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
463
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
464
465
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

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

473
474
475
476

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

477
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
478

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
504
	stdMpi.recv(it.getRank());
505
    }
506

507
508
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
523

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

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


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

    if (!augmentedLagrange)
      return;
  }


551
  void PetscSolverFeti::createIndexB(int component)
552
  {
553
    FUNCNAME("PetscSolverFeti::createIndexB()");
554

555
556

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

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

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

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

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


606
  void PetscSolverFeti::createMatLagrange()
607
608
609
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

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

613
614
    // === Create distributed matrix for Lagrange constraints. ===

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

622
623
624
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

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

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

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

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

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

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

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

682

Thomas Witkowski's avatar
So on    
Thomas Witkowski committed
683
684
685
686
687
688
689
690
691
692
693
694
    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);


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

717

Thomas Witkowski's avatar
Thomas Witkowski committed
718
719
  bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
  {
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
720
721
    FUNCNAME("PetscSolverFeti::testWirebasketEdge()");

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
722
723
    return true;

Thomas Witkowski's avatar
a    
Thomas Witkowski committed
724
725
726
    if (meshDistributor->getMesh()->getDim() == 2)
      return true;

Thomas Witkowski's avatar
Thomas Witkowski committed
727
728
729
    if (meshDistributor->getIntBoundary().getDegreeOwn(edge) != 3)
      return false;

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
730
731
    return true;

Thomas Witkowski's avatar
Thomas Witkowski committed
732
733
734
735
736
737
738
739
740
741
742
743
744
745
    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);
  }
746

747

748
  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseEdges()
749
  {
750
    FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");
751
752

    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
Thomas Witkowski's avatar
Thomas Witkowski committed
753
    std::set<BoundaryObject> allEdges;
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
    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);
      }
    }
774

775
776
777
    int nEdges = allEdges.size();
    mpi::globalAdd(nEdges);
    mpi::globalAdd(nEmptyEdges);
Thomas Witkowski's avatar
Thomas Witkowski committed
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);

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
795
#if 0
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
    std::set<BoundaryObject> allMyEdges;
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(meshDistributor->getMesh(), 0, Mesh::CALL_EL_LEVEL | Mesh::FILL_BOUND);
    while (elInfo) {      
      Element *el = elInfo->getElement();
      for (int i = 0; i < el->getGeo(EDGE); i++) {
	BoundaryObject bobj(el, elInfo->getType(), EDGE, i);
	if (intBound.getDegreeOwn(bobj) == 1 && elInfo->getBoundary(EDGE, i) == INTERIOR) {
	  allMyEdges.insert(bobj);
	}
      }
      elInfo = stack.traverseNext(elInfo);
    }



Thomas Witkowski's avatar
a    
Thomas Witkowski committed
812
813
814
815
816
817
    for (map<int, std::set<BoundaryObject> >::iterator it = allFaces.begin();
	 it != allFaces.end(); ++it) {
      if (it->second.size() == 2) {
	vector<AtomicBoundary> &bs = intBound.getOwn()[it->first];
	for (int i = 0; i < static_cast<int>(bs.size()); i++) {
	  if (bs[i].rankObj.subObj == EDGE &&
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
818
819
	      intBound.getDegreeOwn(bs[i].rankObj) == 1 &&
	      allMyEdges.count(bs[i].rankObj)) {
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
820
821
	    MSG("FOUND AN EDGE: %d %d %d\n", 
		bs[i].rankObj.elIndex, bs[i].rankObj.subObj, bs[i].rankObj.ithObj);
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
822
	    it->second.insert(bs[i].rankObj);
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
823
824
825
826
	  }	      
	}	  
      }
    }
Thomas Witkowski's avatar
da    
Thomas Witkowski committed
827
#endif
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
828
829


830
831
832
833
834
835
836
837
838
839
840
    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);
841
      }
842
843
844
845
      if (oneBoundary.size()) 
	data.push_back(oneBoundary);
      else 
	nEmptyFaces++;
846
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
847

848
    int nFaces = allFaces.size();
849
    mpi::globalAdd(nFaces);
850
    mpi::globalAdd(nEmptyFaces);
851

852
853
854
    MSG("Coarse space faces: %d (empty: %d)\n", nFaces, nEmptyFaces);
   
    return data;
855
856
857
858
859
860
861
862
863
864
865
866
  }


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

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

867
868
869
870
871
    vector<vector<BoundaryObject> > allEdges = getCoarseEdges();
    vector<vector<BoundaryObject> > allFaces = getCoarseFaces(); 
    allEdges.insert(allEdges.end(), allFaces.begin(), allFaces.end());

    nRankEdges = allEdges.size();
872
873
874
875
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

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

877
878
879
    nRankEdges *= componentSpaces.size();
    rStartEdges *= componentSpaces.size();
    nOverallEdges *= componentSpaces.size();
880
881

    MatCreateAIJ(mpiCommGlobal,
882
883
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
884
		 2, PETSC_NULL, 2, PETSC_NULL, 
885
886
887
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

888
    int rowCounter = rStartEdges;
889
890
    for (vector<vector<BoundaryObject> >::iterator it = allEdges.begin(); 
	 it != allEdges.end(); ++it) {
891
      for (int component = 0; component < componentSpaces.size(); component++) {
892
893
894
895
896
	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
897
	  
898
	  TEST_EXIT(edgeDofs.size())("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
899
	  
900
901
902
903
904
905
906
907
908
	  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);
	  }
909
	}
910
911
	rowCounter++;	
      }
912
913
914
915
916
    }

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

Thomas Witkowski's avatar
So on    
Thomas Witkowski committed
917
918
919
920
921
    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);

922
923
924
925
926
927
928
929
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 05a: %.5f seconds (creation of augmented lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
    }
  }


930
  void PetscSolverFeti::createSchurPrimalKsp()