PetscSolverFeti.cc 79.6 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
      subDomainIsLocal(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
40
      subdomain(NULL),
41
      massMatrixSolver(NULL),
42
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
43
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
44
      printTimings(false),
Thomas Witkowski's avatar
d    
Thomas Witkowski committed
45
      dirichletMode(0),
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
92
93
94
95
    {
      MSG("WARNING: CHECK THIS HERE BEFORE GOING INTO RUNNING MULTILEVEL FETI-DP!\n");
      Parameters::get("parallel->level mode", levelMode);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
96
97

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

99
100
101
  }


102
  void PetscSolverFeti::initialize()
103
  {
104
105
    FUNCNAME("PetscSolverFeti::initialize()");

106
    MSG_DBG("Init FETI-DP on mesh level %d\n", meshLevel);
107

108
109
110
    TEST_EXIT_DBG(meshLevel + 2 <= 
		  meshDistributor->getMeshLevelData().getNumberOfLevels())
      ("Mesh hierarchy does not contain at least %d levels!\n", meshLevel + 2);
111

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

114
115
    subDomainIsLocal = (levelData.getMpiComm(meshLevel + 1) == MPI::COMM_SELF);

Thomas Witkowski's avatar
Thomas Witkowski committed
116
    if (subdomain == NULL) {
117
118
119
120
121
122
123
124
125
126
127
      string subSolverInitStr = initFileStr + "->subsolver";
      string solverType = "petsc";
      Parameters::get(subSolverInitStr, solverType);
      OEMSolverCreator *solverCreator =
	dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType, initFileStr));
      TEST_EXIT(solverCreator)
	("No valid solver type found in parameter \"%s\"\n", 
	 initFileStr.c_str());
      solverCreator->setName(subSolverInitStr);
      subdomain = dynamic_cast<PetscSolver*>(solverCreator->create());
      subdomain->initParameters();
128
      subdomain->setSymmetric(isSymmetric);
129
130
      subdomain->setHandleDirichletRows(dirichletMode == 0);
      subdomain->setMeshDistributor(meshDistributor, meshLevel + 1);
131
132

      delete solverCreator;
133
    }
134

135
136
137
138
    primalDofMap.init(componentSpaces, feSpaces);   
    dualDofMap.init(componentSpaces, feSpaces, false);
    localDofMap.init(componentSpaces, feSpaces, !subDomainIsLocal);
    lagrangeMap.init(componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
139

Thomas Witkowski's avatar
Thomas Witkowski committed
140
    if (stokesMode)
141
      interfaceDofMap.init(componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
142

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

147
      interiorDofMap.init(componentSpaces, feSpaces, false);
148
    }
149
150
151
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
152
153
154
155
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
d    
Thomas Witkowski committed
156
157
158
159
160
161
162
163
164
    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
165
166
167
168
    }
  }


169
170
171
172
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
    double timeCounter = MPI::Wtime();

175
176
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

177
178
179
180
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
181
    if (fetiPreconditioner == FETI_DIRICHLET)
182
183
      interiorDofMap.clear();

184
185
    primalDofMap.setDofComm(meshDistributor->getDofComm(meshLevel));
    lagrangeMap.setDofComm(meshDistributor->getDofComm(meshLevel));
186

187
188
189
190
    primalDofMap.setMpiComm(levelData.getMpiComm(meshLevel));
    dualDofMap.setMpiComm(levelData.getMpiComm(meshLevel));
    lagrangeMap.setMpiComm(levelData.getMpiComm(meshLevel));
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel + 1));
191
    if (fetiPreconditioner == FETI_DIRICHLET)
192
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel + 1));
193

194
195
    localDofMap.setDofComm(meshDistributor->getDofComm(meshLevel + 1));
    
Thomas Witkowski's avatar
Thomas Witkowski committed
196
    if (stokesMode) {
197
      interfaceDofMap.clear();
198
199
      interfaceDofMap.setDofComm(meshDistributor->getDofComm(meshLevel));
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0));
200
201
    }

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

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
213

214
    if (fetiPreconditioner == FETI_DIRICHLET)
215
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
216

Thomas Witkowski's avatar
Thomas Witkowski committed
217
    if (stokesMode)
218
219
      interfaceDofMap.update();

220
221
222
    for (int component = 0; component < nComponents; component++) {
      createLagrange(component);
      createAugmentedLagrange(component);
223
    }
224

225
226
227
    lagrangeMap.update();


228
229
    // === ===

230
231
232
    if (subDomainIsLocal) {
      MSG("WARNING: MAKE GENERAL!\n");

233
      rStartInterior = 0;
234
235
236
237
      int localDofs = localDofMap.getOverallDofs();

      mpi::getDofNumbering(domainComm, localDofs,
			   rStartInterior, nGlobalOverallInterior);
238
    } else {
239
240
      MSG("WARNING: MAKE GENERAL!\n");

241
242
243
244
245
246
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

      int groupRowsInterior = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	groupRowsInterior = localDofMap.getOverallDofs();

247
      mpi::getDofNumbering(domainComm, groupRowsInterior,
248
249
250
251
252
253
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	tmp = rStartInterior;

254
255
      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1, 
					MPI_INT, MPI_SUM);
256
257
    }

258

259
260
    for (unsigned int i = 0; i < componentSpaces.size(); i++) {
      const FiniteElemSpace *feSpace = componentSpaces[i];
261
262

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

264
      if (i == pressureComponent) {
Thomas Witkowski's avatar
Thomas Witkowski committed
265
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
266
267
	    interfaceDofMap[i].nRankDofs, 
	    interfaceDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
268
      } else {
269
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
270
271
272
	    primalDofMap[i].nRankDofs, 
	    primalDofMap[i].nLocalDofs,
	    primalDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
273
274
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
275
276
	    dualDofMap[i].nRankDofs, 
	    dualDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
279
280
	    lagrangeMap[i].nRankDofs, 
	    lagrangeMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
281
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
282
    }
283

Thomas Witkowski's avatar
Thomas Witkowski committed
284
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
287
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
288
289

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


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


304
  void PetscSolverFeti::createPrimals(int component)
Thomas Witkowski's avatar
Thomas Witkowski committed
305
  {
306
    FUNCNAME("PetscSolverFeti::createPrimals()");  
307

308
    if (component == pressureComponent)
309
310
      return;

311
312
    const FiniteElemSpace *feSpace = componentSpaces[component];

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

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

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

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

334
335
336
337
	if ((fabs(c[0]) < e && fabs(c[1] - 12.5) < e) ||
	    (fabs(c[0] - 25.0) < e && fabs(c[1] - 12.5) < e) ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1]) < e) ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 25.0) < e) ||
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
338
339
340
341
342
	    (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);
	}
      }
343
    }
344

345
346
347
348

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

349
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it) {
350
      if (dofMap[feSpace].isRankDof(*it)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
351
	primalDofMap[component].insertRankDof(*it);
352
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
353
  	primalDofMap[component].insertNonRankDof(*it);
354
355
      }
    }
356
357
358
  }


359
  void PetscSolverFeti::createDuals(int component)
360
361
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
362

363
    if (component == pressureComponent)
364
365
      return;

366
367
    const FiniteElemSpace *feSpace = componentSpaces[component];

368
369
370
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
371
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
372
    
373
    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;

Thomas Witkowski's avatar
Thomas Witkowski committed
378
      if (isPrimal(component, **it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
379
	continue;
380
381
      
      if (subDomainIsLocal || dofMapSd[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
382
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
383
    }
384
385
386
  }

  
387
  void PetscSolverFeti::createInterfaceNodes(int component)
388
389
390
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

391
    if (component != pressureComponent)
392
393
      return;

394
    const FiniteElemSpace *feSpace = componentSpaces[component];
395
396
397
398
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
399
	 it != allBoundaryDofs.end(); ++it) {
400
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
401
402
	continue;      
      
403
      if (dofMap[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
404
	interfaceDofMap[component].insertRankDof(**it);
405
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
406
	interfaceDofMap[component].insertNonRankDof(**it);      
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
407
    }
408
409
410
  }


411
  void PetscSolverFeti::createLagrange(int component)
412
413
414
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

415
    if (component == pressureComponent)
416
417
      return;

418
    const FiniteElemSpace *feSpace = componentSpaces[component];
419
420
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
423
    // 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.
424
    map<int, std::set<DegreeOfFreedom> > subDomainRankDofs;
425

426
427
    if (not subDomainIsLocal) {
      StdMpi<vector<int> > stdMpi(domainComm);
428

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

432
433
	vector<int> dofs;
	dofs.reserve(it.getDofs().size());
434
435

	for (; !it.endDofIter(); it.nextDof()) {
436
	  if (dofMapSd[feSpace].isRankDof(it.getDofIndex()))
437
	    dofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
438
	  else
439
	    dofs.push_back(0);
440
441
	}

442
	stdMpi.send(it.getRank(), dofs);
443
444
      }	     

445
      for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getSendDofs(), feSpace);
446
447
448
449
450
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

451
      for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getSendDofs(), feSpace); 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
452
453
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
454
455
456
	  if (!isPrimal(component, it.getDofIndex()) &&
	      stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	    subDomainRankDofs[it.getRank()].insert(it.getDofIndex());
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
457
    }
458

Thomas Witkowski's avatar
Thomas Witkowski committed
459
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
460
461
462
463
464
465
      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)).         ===

466
467
    int mpiRank = domainComm.Get_rank();
    for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getSendDofs(), feSpace); 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
468
469
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
470
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
471
472
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

473
474
475
476
477
	  // If the subdomain is local, always add the counterpart rank, 
	  // otherwise check if the other rank is the owner of the DOF in 
	  // its subdomain.
 	  if (subDomainIsLocal || 
	      subDomainRankDofs[it.getRank()].count(it.getDofIndex()))
478
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
479
480
481
482
	}
      }
    }

483
484
485
486

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

487
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm(meshLevel));
488

489
    for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getSendDofs(), feSpace);
490
491
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
492
	if (!isPrimal(component, it.getDofIndex()))
493
 	  if (subDomainIsLocal || subDomainRankDofs[it.getRank()].count(it.getDofIndex()))
494
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
495
496
497

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
511
	stdMpi.recv(it.getRank());
512
    }
513

514
515
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
529

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

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


548
  void PetscSolverFeti::createAugmentedLagrange(int component)
549
550
551
552
553
554
555
556
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


557
  void PetscSolverFeti::createIndexB(int component)
558
  {
559
    FUNCNAME("PetscSolverFeti::createIndexB()");
560

561
    const FiniteElemSpace *feSpace = componentSpaces[component];
562
    DOFAdmin* admin = feSpace->getAdmin();
563
564
565
566

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

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

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

597
598
    for (DofMap::iterator it = dualDofMap[component].getMap().begin();
	 it != dualDofMap[component].getMap().end(); ++it) {
599
      if (subDomainIsLocal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
600
	localDofMap[component].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
601
      } else {
602
	if (dofMapSd[feSpace].isRankDof(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
603
	  localDofMap[component].insertRankDof(it->first);
604
	else 
Thomas Witkowski's avatar
Thomas Witkowski committed
605
	  localDofMap[component].insertNonRankDof(it->first);
606
      }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
607
    }
608
609
610
  }


611
  void PetscSolverFeti::createMatLagrange()
612
613
614
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
615
    double wtime = MPI::Wtime();
616
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
617

618
619
    // === Create distributed matrix for Lagrange constraints. ===

620
    MatCreateAIJ(domainComm,
621
622
623
624
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
625
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
626

627

628
629
630
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

631

632
633
634
635
636
637
638
    // === 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.                                                     ===

639
    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
640
      DofMap &dualMap = dualDofMap[component].getMap();
641

642
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
643
	TEST_EXIT_DBG(boundaryDofRanks[componentSpaces[component]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
646
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
647
	int index = lagrangeMap.getMatIndex(component, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
648

Thomas Witkowski's avatar
Thomas Witkowski committed
649
	// Copy set of all ranks that contain this dual node.
650
651
	vector<int> W(boundaryDofRanks[componentSpaces[component]][it->first].begin(), 
		      boundaryDofRanks[componentSpaces[component]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
652
653
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
654
655

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

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
674
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
675
676
677
678
679
680
	  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);
681
682
683
684
685
686
687
	  }
	}
      }
    }

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

689

690
691
692
693
694
695
    {
      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);
    }
Thomas Witkowski's avatar
So on    
Thomas Witkowski committed
696
697


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

720

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

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
725
726
    return true;

Thomas Witkowski's avatar
a    
Thomas Witkowski committed
727
728
729
    if (meshDistributor->getMesh()->getDim() == 2)
      return true;

730
    if (meshDistributor->getIntBoundary(meshLevel).getDegreeOwn(edge) != 3)
Thomas Witkowski's avatar
Thomas Witkowski committed
731
732
      return false;

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
733
734
    return true;

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

750

751
  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseEdges()
752
  {
753
    FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");
754

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

778
779
780
    int nEdges = allEdges.size();
    mpi::globalAdd(nEdges);
    mpi::globalAdd(nEmptyEdges);
Thomas Witkowski's avatar
Thomas Witkowski committed
781

782
783
784
785
786
787
788
789
790
791
    MSG("Coarse space edges: %d (empty: %d)\n", nEdges, nEmptyEdges);
    
    return data;
  }


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

792
    InteriorBoundary &intBound = meshDistributor->getIntBoundary(meshLevel);
793
794
795
796
797
    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
798
#if 0
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
    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
815
816
817
818
819
820
    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
821
822
	      intBound.getDegreeOwn(bs[i].rankObj) == 1 &&
	      allMyEdges.count(bs[i].rankObj)) {
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
823
824
	    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
825
	    it->second.insert(bs[i].rankObj);
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
826
827
828
829
	  }	      
	}	  
      }
    }
Thomas Witkowski's avatar
da    
Thomas Witkowski committed
830
#endif
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
831
832


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

851
    int nFaces = allFaces.size();
852
    mpi::globalAdd(nFaces);
853
    mpi::globalAdd(nEmptyFaces);
854

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


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

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

870
871
872
873
874
    vector<vector<BoundaryObject> > allEdges = getCoarseEdges();
    vector<vector<BoundaryObject> > allFaces = getCoarseFaces(); 
    allEdges.insert(allEdges.end(), allFaces.begin(), allFaces.end());

    nRankEdges = allEdges.size();
875
    int rStartEdges = 0;
876
    mpi::getDofNumbering(domainComm, nRankEdges, rStartEdges, nOverallEdges);
877
878

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

880
881
882
    nRankEdges *= componentSpaces.size();
    rStartEdges *= componentSpaces.size();
    nOverallEdges *= componentSpaces.size();
883

Thomas Witkowski's avatar