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
355
356
	    (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);
	}
      }
357
    }
358

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

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


372
  void PetscSolverFeti::createDuals(int component)
373
374
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
375

376
    if (component == pressureComponent)
377
378
      return;

379
380
    const FiniteElemSpace *feSpace = componentSpaces[component];

381
382
383
    // === Create global index of the dual nodes on each rank. ===

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

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

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

  
403
  void PetscSolverFeti::createInterfaceNodes(int component)
404
405
406
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

407
    if (component != pressureComponent)
408
409
      return;

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

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


427
  void PetscSolverFeti::createLagrange(int component)
428
429
430
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

431
    if (component == pressureComponent)
432
433
      return;

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

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

442
443
    if (not subDomainIsLocal) {
      StdMpi<vector<int> > stdMpi(domainComm);
444

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

448
449
	vector<int> dofs;
	dofs.reserve(it.getDofs().size());
450
451

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

458
	stdMpi.send(it.getRank(), dofs);
459
460
      }	     

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

      stdMpi.startCommunication();

467
      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())
470
471
472
	  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
473
    }
474

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

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

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

499
500
501
502

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

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

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
527
	stdMpi.recv(it.getRank());
528
    }
529

530
531
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
545

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

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


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

    if (!augmentedLagrange)
      return;
  }


573
  void PetscSolverFeti::createIndexB(int component)
574
  {
575
    FUNCNAME("PetscSolverFeti::createIndexB()");
576

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

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

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

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

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

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


630
  void PetscSolverFeti::createMatLagrange()
631
632
633
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

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

637
638
    // === Create distributed matrix for Lagrange constraints. ===

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

646

647
648
649
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

650

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

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

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

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

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

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

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

708

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

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


722
723
724
725
726
    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

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

746

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

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
751
752
    return true;

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

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

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
759
760
    return true;

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

776

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

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

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

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


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

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


859
860
861
862
863
864
865
866
867
868
869
    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);
870
      }
871
872
873
874
      if (oneBoundary.size()) 
	data.push_back(oneBoundary);
      else 
	nEmptyFaces++;