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
  void PetscSolverFeti::createDuals(int component)
343
344
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
345

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

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

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

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

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

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

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

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

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

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


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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
410
    // Stores for all rank owned communication DOFs, if the counterpart is
    // a rank owned DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFs that fulfill this requirenment.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
411
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
412

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

      for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank()) {

	vector<int> subdomainRankDofs;
	subdomainRankDofs.reserve(it.getDofs().size());

	for (; !it.endDofIter(); it.nextDof()) {
424
	  if (dofMapSd[feSpace].isRankDof(it.getDofIndex()))
425
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
426
	  else
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
	    subdomainRankDofs.push_back(0);
	}

	stdMpi.send(it.getRank(), subdomainRankDofs);
      }	     

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace); 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
442
443
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
444
	  if (!isPrimal(component, it.getDofIndex()))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
445
446
447
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
448

Thomas Witkowski's avatar
Thomas Witkowski committed
449
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
450
451
452
453
454
455
      return;


    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===

456
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
457
458
459
460
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
461
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
462
463
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

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

471
472
473
474

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

475
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
476

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
502
	stdMpi.recv(it.getRank());
503
    }
504

505
506
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
521

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

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


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

    if (!augmentedLagrange)
      return;
  }


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

553
554

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

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

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

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

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


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

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

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

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

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

623
624
625
626
627
628
629
    // === Create for all duals the corresponding Lagrange constraints. On ===
    // === each rank we traverse all pairs (n, m) of ranks, with n < m,    ===
    // === that contain this node. If the current rank number is r, and    ===
    // === n == r, the rank sets 1.0 for the corresponding constraint, if  ===
    // === m == r, than the rank sets -1.0 for the corresponding           ===
    // === constraint.                                                     ===

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

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

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

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

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
665
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
666
667
668
669
670
671
	  int nConstraints = (degree * (degree - 1)) / 2;
	  for (int i = 0; i < nConstraints; i++) {
	    VecSetValue(vec_scale_lagrange,
			index + i,
			1.0 / static_cast<double>(degree),
			INSERT_VALUES);
672
673
674
675
676
677
678
	  }
	}
      }
    }

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

680

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

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


693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

    if (fetiPreconditioner != FETI_NONE || stokesMode) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatDiagonalScale(mat_lagrange_scaled, vec_scale_lagrange, PETSC_NULL);
    }

    VecDestroy(&vec_scale_lagrange);


    // === Print final timings. ===

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
708
709
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
710
711
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
712
    }
713
714
  }

715

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

Thomas Witkowski's avatar
da  
Thomas Witkowski committed
720
721
    return true;

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

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

Thomas Witkowski's avatar
da  
Thomas Witkowski committed
728
729
    return true;

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

745

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

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

773
774
775
    int nEdges = allEdges.size();
    mpi::globalAdd(nEdges);
    mpi::globalAdd(nEmptyEdges);
Thomas Witkowski's avatar
Thomas Witkowski committed
776

777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
    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
793
#if 0
Thomas Witkowski's avatar
a  
Thomas Witkowski committed
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
    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
810
811
812
813
814
815
    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
816
817
	      intBound.getDegreeOwn(bs[i].rankObj) == 1 &&
	      allMyEdges.count(bs[i].rankObj)) {
Thomas Witkowski's avatar
a  
Thomas Witkowski committed
818
819
	    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
820
	    it->second.insert(bs[i].rankObj);
Thomas Witkowski's avatar
a  
Thomas Witkowski committed
821
822
823
824
	  }	      
	}	  
      }
    }
Thomas Witkowski's avatar
da  
Thomas Witkowski committed
825
#endif
Thomas Witkowski's avatar
a  
Thomas Witkowski committed
826
827


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

846
    int nFaces = allFaces.size();
847
    mpi::globalAdd(nFaces);
848
    mpi::globalAdd(nEmptyFaces);
849

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


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

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

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

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

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

875
876
877
    nRankEdges *= componentSpaces.size();
    rStartEdges *= componentSpaces.size();
    nOverallEdges *= componentSpaces.size();
878
879

    MatCreateAIJ(mpiCommGlobal,
880
881
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),