PetscSolverFeti.cc 84.9 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
      fetiSolverType(EXACT),
33
34
35
36
37
38
      primalDofMap(COMPONENT_WISE),
      dualDofMap(COMPONENT_WISE),
      interfaceDofMap(COMPONENT_WISE),
      localDofMap(COMPONENT_WISE),
      lagrangeMap(COMPONENT_WISE),
      interiorDofMap(COMPONENT_WISE),
39
      schurPrimalSolver(0),
40
      subDomainIsLocal(true),
Thomas Witkowski's avatar
Thomas Witkowski committed
41
      subdomain(NULL),
42
      massMatrixSolver(NULL),
Thomas Witkowski's avatar
Thomas Witkowski committed
43
      printTimings(false),
Thomas Witkowski's avatar
d    
Thomas Witkowski committed
44
      dirichletMode(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
45
      stokesMode(false),
46
      augmentedLagrange(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
47
      pressureComponent(-1)
48
49
50
51
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

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

90
91
92
93
    {
      MSG("WARNING: CHECK THIS HERE BEFORE GOING INTO RUNNING MULTILEVEL FETI-DP!\n");
      Parameters::get("parallel->level mode", levelMode);
    }
94
95
96
97
98
99
100
101
102
    
    {
      int tmp = 0;
      Parameters::get(initFileStr + "->feti->inexact", tmp);
      if (tmp == 1)
	fetiSolverType = INEXACT;
      if (tmp == 2)
	fetiSolverType = INEXACT_REDUCED;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
103
104

    Parameters::get("parallel->print timings", printTimings);
105
106
107
  }


108
  void PetscSolverFeti::initialize()
109
  {
110
111
    FUNCNAME("PetscSolverFeti::initialize()");

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

114
115
116
    TEST_EXIT_DBG(meshLevel + 2 <= 
		  meshDistributor->getMeshLevelData().getNumberOfLevels())
      ("Mesh hierarchy does not contain at least %d levels!\n", meshLevel + 2);
117

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

120
121
    subDomainIsLocal = (levelData.getMpiComm(meshLevel + 1) == MPI::COMM_SELF);

Thomas Witkowski's avatar
Thomas Witkowski committed
122
    if (subdomain == NULL) {
123
124
125
126
127
128
129
130
131
132
133
      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();
134
      subdomain->setSymmetric(isSymmetric);
135
136
      subdomain->setHandleDirichletRows(dirichletMode == 0);
      subdomain->setMeshDistributor(meshDistributor, meshLevel + 1);
137
      subdomain->init(componentSpaces, feSpaces);
138
139

      delete solverCreator;
140
    }
141

142
143
144
145
    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
146

Thomas Witkowski's avatar
Thomas Witkowski committed
147
    if (stokesMode)
148
      interfaceDofMap.init(componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
149

150
    if (fetiPreconditioner == FETI_DIRICHLET) {
151
      TEST_EXIT(levelMode == 1)
152
153
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

154
      interiorDofMap.init(componentSpaces, feSpaces, false);
155
    }
156
157
158
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
159
160
161
162
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
d    
Thomas Witkowski committed
163
164
165
166
167
168
169
170
171
    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
172
173
174
175
    }
  }


176
177
178
179
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
180
181
    double timeCounter = MPI::Wtime();

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

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

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

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

201
202
    localDofMap.setDofComm(meshDistributor->getDofComm(meshLevel + 1));
    
Thomas Witkowski's avatar
Thomas Witkowski committed
203
    if (stokesMode) {
204
      interfaceDofMap.clear();
205
206
      interfaceDofMap.setDofComm(meshDistributor->getDofComm(meshLevel));
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0));
207
208
    }

209
210
211
212
213
214
    int nComponents = componentSpaces.size();
    for (int component = 0; component < nComponents; component++) {
      createPrimals(component);  
      createDuals(component);
      createInterfaceNodes(component);
      createIndexB(component);
215
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
216
217
218
219

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

221
    if (fetiPreconditioner == FETI_DIRICHLET)
222
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
223

Thomas Witkowski's avatar
Thomas Witkowski committed
224
    if (stokesMode)
225
226
      interfaceDofMap.update();

227
228
229
    for (int component = 0; component < nComponents; component++) {
      createLagrange(component);
      createAugmentedLagrange(component);
230
    }
231

232
233
234
    lagrangeMap.update();


235
236
    // === ===

237
238
239
    if (subDomainIsLocal) {
      MSG("WARNING: MAKE GENERAL!\n");

240
      rStartInterior = 0;
241
242
243
244
      int localDofs = localDofMap.getOverallDofs();

      mpi::getDofNumbering(domainComm, localDofs,
			   rStartInterior, nGlobalOverallInterior);
245
    } else {
246
247
      MSG("WARNING: MAKE GENERAL!\n");

248
249
250
251
252
253
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

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

254
      mpi::getDofNumbering(domainComm, groupRowsInterior,
255
256
257
258
259
260
			   rStartInterior, nGlobalOverallInterior);

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

261
262
      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1, 
					MPI_INT, MPI_SUM);
263
264
    }

265

266
267
    for (unsigned int i = 0; i < componentSpaces.size(); i++) {
      const FiniteElemSpace *feSpace = componentSpaces[i];
268
269

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
285
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
286
287
	    lagrangeMap[i].nRankDofs, 
	    lagrangeMap[i].nOverallDofs);
288
289
290
291

	MSG("  nRankLocal = %d  nOverallLocal = %d\n",
	    localDofMap[i].nRankDofs,
	    localDofMap[i].nOverallDofs);       
Thomas Witkowski's avatar
Thomas Witkowski committed
292
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
293
    }
294

Thomas Witkowski's avatar
Thomas Witkowski committed
295
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
298
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300

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


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


315
  void PetscSolverFeti::createPrimals(int component)
Thomas Witkowski's avatar
Thomas Witkowski committed
316
  {
317
    FUNCNAME("PetscSolverFeti::createPrimals()");  
318

319
    if (component == pressureComponent)
320
321
      return;

322
323
    const FiniteElemSpace *feSpace = componentSpaces[component];

324
325
326
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
327
    // Set of DOF indices that are considered to be primal variables.
328
    DofContainerSet& vertices = 
329
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
330
331

    DofIndexSet primals;
332
    for (DofContainerSet::iterator it = vertices.begin(); 
333
	 it != vertices.end(); ++it) {
334
      
335
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
336
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
337

338
339
340
      if (meshLevel == 1 && not (*interiorMap)[component].isSet(**it))
	continue;

341
      if (subDomainIsLocal) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
342
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
343
      } else {
344
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
347
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

348
349
350
351
	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
352
353
354
	    (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);
355
356
	} else {
	  MSG("OMMIT SOME PRIMAL!\n");
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
357
358
	}
      }
359
    }
360

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

364
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it) {
365
      if (dofMap[feSpace].isRankDof(*it)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
366
	primalDofMap[component].insertRankDof(*it);
367
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
368
  	primalDofMap[component].insertNonRankDof(*it);
369
370
      }
    }
371
372
373
  }


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

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

381
382
    const FiniteElemSpace *feSpace = componentSpaces[component];

383
384
385
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
386
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
387
    
388
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
389
	 it != allBoundaryDofs.end(); ++it) {
390
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
391
392
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
393
      if (isPrimal(component, **it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
394
	continue;
395
396
397

      if (meshLevel == 1 && not (*interiorMap)[component].isSet(**it))
	continue;
398
      
399
      if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
400
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
401
    }
402
403
404
  }

  
405
  void PetscSolverFeti::createInterfaceNodes(int component)
406
407
408
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

409
    if (component != pressureComponent)
410
411
      return;

412
    const FiniteElemSpace *feSpace = componentSpaces[component];
413
414
415
416
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
417
	 it != allBoundaryDofs.end(); ++it) {
418
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
419
420
	continue;      
      
421
      if (dofMap[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
422
	interfaceDofMap[component].insertRankDof(**it);
423
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
424
	interfaceDofMap[component].insertNonRankDof(**it);      
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
425
    }
426
427
428
  }


429
  void PetscSolverFeti::createLagrange(int component)
430
431
432
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

433
    if (component == pressureComponent)
434
435
      return;

436
    const FiniteElemSpace *feSpace = componentSpaces[component];
437
438
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
439
440
441
    // 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.
442
    map<int, std::set<DegreeOfFreedom> > subDomainRankDofs;
443

444
445
    if (not subDomainIsLocal) {
      StdMpi<vector<int> > stdMpi(domainComm);
446

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

450
451
	vector<int> dofs;
	dofs.reserve(it.getDofs().size());
452
453

	for (; !it.endDofIter(); it.nextDof()) {
454
	  if (dofMapSubDomain[feSpace].isRankDof(it.getDofIndex()))
455
	    dofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
456
	  else
457
	    dofs.push_back(0);
458
459
	}

460
	stdMpi.send(it.getRank(), dofs);
461
462
      }	     

463
      for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getSendDofs(), feSpace);
464
465
466
467
468
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

469
      for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getSendDofs(), feSpace); 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
470
471
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
472
473
474
	  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
475
    }
476

Thomas Witkowski's avatar
Thomas Witkowski committed
477
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
478
479
480
481
482
483
      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)).         ===

484
485
    int mpiRank = domainComm.Get_rank();
    for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getSendDofs(), feSpace); 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
486
487
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
488
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
489
490
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

491
492
493
494
495
	  // 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()))
496
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
497
498
499
500
	}
      }
    }

501
502
503
504

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

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

507
    for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getSendDofs(), feSpace);
508
509
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
510
	if (!isPrimal(component, it.getDofIndex()))
511
 	  if (subDomainIsLocal || subDomainRankDofs[it.getRank()].count(it.getDofIndex()))
512
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
513
514
515

    stdMpi.updateSendDataSize();

516
    for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getRecvDofs(), feSpace); 
517
	 !it.end(); it.nextRank()) {
518
      bool recvFromRank = false;
519
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
520
	if (!isPrimal(component, it.getDofIndex())) {
521
 	  if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
522
523
524
	    recvFromRank = true;
	    break;
	  }
525
	}
526
      }
527
528

      if (recvFromRank)
529
	stdMpi.recv(it.getRank());
530
    }
531

532
533
    stdMpi.startCommunication();

534
    for (DofComm::Iterator it(meshDistributor->getDofComm(meshLevel).getRecvDofs(), feSpace); 
535
	 !it.end(); it.nextRank()) {
536
      int i = 0;
537
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
538
	if (!isPrimal(component, it.getDofIndex()))
539
 	  if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(it.getDofIndex()))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
540
541
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
542
	  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
543
	    lagrangeMap[component].insertNonRankDof(it.getDofIndex());
544
	  }
545
546
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
547

548
549
550
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

551
    int nRankLagrange = 0;
552
    DofMap& dualMap = dualDofMap[component].getMap();
553
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
554
      if (dofMap[feSpace].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
555
	lagrangeMap[component].insertRankDof(it->first, nRankLagrange);
556
	int degree = boundaryDofRanks[feSpace][it->first].size();
557
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
558
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
559
	lagrangeMap[component].insertNonRankDof(it->first);
560
561
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
562
    lagrangeMap[component].nRankDofs = nRankLagrange;
563
564
565
  }


566
  void PetscSolverFeti::createAugmentedLagrange(int component)
567
568
569
570
571
572
573
574
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


575
  void PetscSolverFeti::createIndexB(int component)
576
  {
577
    FUNCNAME("PetscSolverFeti::createIndexB()");
578

579
    const FiniteElemSpace *feSpace = componentSpaces[component];
580
    DOFAdmin* admin = feSpace->getAdmin();
581
582
583
584

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

586
    int nLocalInterior = 0;    
587
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
588
      if (admin->isDofFree(i) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
589
590
591
	  isPrimal(component, i) ||
	  isDual(component, i) ||
	  isInterface(component, i) ||
592
	  dirichletRows[component].count(i))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
593
	continue;      
594

595
596
597
      if (meshLevel == 1 &&  not (*interiorMap)[component].isSet(i))
	continue;

598
      if (subDomainIsLocal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
599
	localDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
600
601
	
	if (fetiPreconditioner == FETI_DIRICHLET)
602
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
603
604
605
	
	nLocalInterior++;	
      } else {
606
	if (dofMapSubDomain[feSpace].isRankDof(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
607
	  localDofMap[component].insertRankDof(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
608
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
609
	  localDofMap[component].insertNonRankDof(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
610
611
612
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
613
      }
614
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
615
    
616
617
    // === And finally, add the global indicies of all dual nodes. ===

618
619
    for (DofMap::iterator it = dualDofMap[component].getMap().begin();
	 it != dualDofMap[component].getMap().end(); ++it) {
620
      if (subDomainIsLocal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
621
	localDofMap[component].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
622
      } else {
623
	if (dofMapSubDomain[feSpace].isRankDof(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
624
	  localDofMap[component].insertRankDof(it->first);
625
	else 
Thomas Witkowski's avatar
Thomas Witkowski committed
626
	  localDofMap[component].insertNonRankDof(it->first);
627
      }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
628
    }
629
630
631
  }


632
  void PetscSolverFeti::createMatLagrange()
633
634
635
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
636
    double wtime = MPI::Wtime();
637
    int mpiRank = domainComm.Get_rank();
Thomas Witkowski's avatar
Thomas Witkowski committed
638

639
640
    // === Create distributed matrix for Lagrange constraints. ===

641
    MatCreateAIJ(domainComm,
642
643
644
645
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
646
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
647

648

649
650
651
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

652

653
654
655
656
657
658
659
    // === 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.                                                     ===

660
    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
661
      DofMap &dualMap = dualDofMap[component].getMap();
662

663
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
664
	TEST_EXIT_DBG(boundaryDofRanks[componentSpaces[component]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
665
666
667
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
668
	int index = lagrangeMap.getMatIndex(component, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
669

Thomas Witkowski's avatar
Thomas Witkowski committed
670
	// Copy set of all ranks that contain this dual node.
671
672
	vector<int> W(boundaryDofRanks[componentSpaces[component]][it->first].begin(), 
		      boundaryDofRanks[componentSpaces[component]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
673
674
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
675
676

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
677
	
678
	int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
681
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
682
683
	      MatSetValue(mat_lagrange, 
			  index + counter, 
684
			  localDofMap.getMatIndex(component, it->first) + rStartInterior,
685
686
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
687
	    }
688
689
690
691
692
693
694
	    counter++;
	  }
	}

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
695
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
696
697
698
699
700
701
	  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);
702
703
704
705
706
707
708
	  }
	}
      }
    }

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

710

711
#if (DEBUG != 0)
712
713
714
715
    {
      int nZeroRows = PetscSolverFetiDebug::testZeroRows(mat_lagrange);
      int m,n;
      MatGetSize(mat_lagrange, &m ,&n);
716
      mpi::globalAdd(domainComm, nZeroRows);
717
      MSG("Lagrange matrix has %d zero rows and global size of %d %d!\n", nZeroRows, m, n);
718
719

      TEST_EXIT(nZeroRows == 0)("Lagrange matrix has zero rows!\n");
720
    }
721
#endif
Thomas Witkowski's avatar
So on    
Thomas Witkowski committed
722
723


724
725
726
727
728
    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

729
730
731
    if (fetiPreconditioner != FETI_NONE || 
	fetiSolverType == INEXACT || 
	stokesMode) {
732
733
734
735
736
737
738
739
740
      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
741
742
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
745
    }
746
747
  }

748

Thomas Witkowski's avatar
Thomas Witkowski committed
749
750
  bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
  {
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
751
752
    FUNCNAME("PetscSolverFeti::testWirebasketEdge()");

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
753
754
    return true;

Thomas Witkowski's avatar
a    
Thomas Witkowski committed
755
756
757
    if (meshDistributor->getMesh()->getDim() == 2)
      return true;

758
    if (meshDistributor->getIntBoundary(meshLevel).getDegreeOwn(edge) != 3)
Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
      return false;

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
761
762
    return true;

Thomas Witkowski's avatar
Thomas Witkowski committed
763
764
765
766
767
768
769
770
771
772
773
774
775
776
    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);
  }
777

778

779
  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseEdges()
780
  {
781
    FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");
782

783
    InteriorBoundary &intBound = meshDistributor->getIntBoundary(meshLevel);
Thomas Witkowski's avatar
Thomas Witkowski committed
784
    std::set<BoundaryObject> allEdges;
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
    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);
      }
    }
805

806
807
808
    int nEdges = allEdges.size();
    mpi::globalAdd(nEdges);
    mpi::globalAdd(nEmptyEdges);
Thomas Witkowski's avatar
Thomas Witkowski committed
809

810
811
812
813
814
815
816
817
818
819
    MSG("Coarse space edges: %d (empty: %d)\n", nEdges, nEmptyEdges);
    
    return data;
  }


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

820
    InteriorBoundary &intBound = meshDistributor->getIntBoundary(meshLevel);
821
822
823
824
825
    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
826
#if 0
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
    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
843
844
845
846
847
848
    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
849
850
	      intBound.getDegreeOwn(bs[i].rankObj) == 1 &&
	      allMyEdges.count(bs[i].rankObj)) {
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
851
852
	    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
853
	    it->second.insert(bs[i].rankObj);
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
854
855
856
857
	  }	      
	}	  
      }
    }
Thomas Witkowski's avatar
da    
Thomas Witkowski committed
858
#endif
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
859
860


861
862
863
864
865
866
867
868
869
870
871
    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);
872
      }
Thomas Witkowski's avatar