PetscSolverFeti.cc 85.8 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20
21


Thomas Witkowski's avatar
Thomas Witkowski committed
22
#include "MatrixVector.h"
23
#include "parallel/PetscHelper.h"
24
#include "parallel/PetscSolverFeti.h"
25
26
#include "parallel/PetscSolverFetiDebug.h"
#include "parallel/PetscSolverFetiMonitor.h"
27
#include "parallel/PetscSolverFetiStructs.h"
28
29
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
30
31
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
32
#include "parallel/PetscSolverGlobalMatrix.h"
33
#include "io/VtkWriter.h"
34
#include "utility/PetscWrapper.h"
35

36
namespace AMDiS { namespace Parallel {
37
38
39

  using namespace std;

40
  PetscSolverFeti::PetscSolverFeti(string name)
41
    : PetscSolver(name),
42
      fetiSolverType(EXACT),
43
      dofMapSubDomain(FESPACE_WISE, true),
44
45
46
47
48
49
      primalDofMap(COMPONENT_WISE),
      dualDofMap(COMPONENT_WISE),
      interfaceDofMap(COMPONENT_WISE),
      localDofMap(COMPONENT_WISE),
      lagrangeMap(COMPONENT_WISE),
      interiorDofMap(COMPONENT_WISE),
50
      schurPrimalSolver(0),
51
      levelMode(1),
52
      subDomainIsLocal(true),
53
54
      subdomain(NULL),
      massMatrixSolver(NULL),
Thomas Witkowski's avatar
Thomas Witkowski committed
55
      printTimings(false),
56
57
58
      augmentedLagrange(false),
      nRankEdges(0),
      nOverallEdges(0),
Thomas Witkowski's avatar
d    
Thomas Witkowski committed
59
      dirichletMode(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
60
      stokesMode(false),
61
      pressureComponent(-1)
62
63
64
65
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
66
    Parameters::get(name + "->left precon", preconditionerName);
67
    if (preconditionerName == "" || preconditionerName == "no") {
68
      MSG("Create FETI-DP solver with no preconditioner!\n");
69
70
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
71
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
72
73
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
74
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
75
76
      fetiPreconditioner = FETI_LUMPED;
    } else {
77
      ERROR_EXIT("Preconditioner \"%s\" not available!\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
78
		 preconditionerName.c_str());
79
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
80

81
    preconditionerName = "";
82
    Parameters::get(name + "->right precon", preconditionerName);
83
84
    if (preconditionerName != "" && preconditionerName != "no") {
      ERROR_EXIT("FETI-DP does not support right preconditioning! (parameter \"%s->right precon\" has value \"%s\")\n",
85
		 name.c_str(), preconditionerName.c_str());
86
87
    }

88
    Parameters::get(name + "->feti->schur primal solver", schurPrimalSolver);
Thomas Witkowski's avatar
Thomas Witkowski committed
89
    TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
90
      ("Wrong solver \"%d\"for the Schur primal complement!\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
91
       schurPrimalSolver);
92

93
    Parameters::get(name + "->feti->stokes mode", stokesMode);
94
    if (stokesMode) {
95
      Parameters::get(name + "->feti->pressure component", pressureComponent);
96
97
98
      TEST_EXIT(pressureComponent >= 0)
	("FETI-DP in Stokes mode, no pressure component defined!\n");
    }
99

100
    Parameters::get(name + "->feti->augmented lagrange", augmentedLagrange);
101

102
    Parameters::get(name + "->feti->symmetric", isSymmetric);
103

104
105
106
107
    {
      MSG("WARNING: CHECK THIS HERE BEFORE GOING INTO RUNNING MULTILEVEL FETI-DP!\n");
      Parameters::get("parallel->level mode", levelMode);
    }
108

109
110
    {
      int tmp = 0;
111
      Parameters::get(name + "->feti->inexact", tmp);
112
113
114
115
116
      if (tmp == 1)
	fetiSolverType = INEXACT;
      if (tmp == 2)
	fetiSolverType = INEXACT_REDUCED;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
117
118

    Parameters::get("parallel->print timings", printTimings);
119
  }
120

121
122
123
124
125
126
127
128

  void PetscSolverFeti::init(vector<const FiniteElemSpace*> &fe0,
			     vector<const FiniteElemSpace*> &fe1,
			     bool createGlobalMapping)
  {
    FUNCNAME_DBG("PetscSolverFeti::init()");

    super::init(fe0, fe1, createGlobalMapping);
129

130
131
132
133
134
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
    int nLevels = levelData.getNumberOfLevels();
    TEST_EXIT_DBG(nLevels >= 1)("nLevels < 1! Should not happen!\n");

    if (createGlobalMapping) {
135
      if (meshLevel + 1 < nLevels &&
136
137
138
	  levelData.getMpiComm(meshLevel + 1) != MPI::COMM_SELF) {
	dofMapSubDomain.init(componentSpaces, feSpaces);
	dofMapSubDomain.setMpiComm(levelData.getMpiComm(meshLevel + 1));
139
	dofMapSubDomain.setDofComms(meshDistributor->getDofComms(), meshLevel + 1);
140
141
142
143
144
	dofMapSubDomain.clear();
	meshDistributor->registerDofMap(dofMapSubDomain);
      }
    }
  }
145
146


147
  void PetscSolverFeti::initialize()
148
  {
149
150
    FUNCNAME("PetscSolverFeti::initialize()");

151
#ifndef NDEBUG
152
153
    MSG("Init FETI-DP on mesh level %d\n", meshLevel);
#endif
154

155
    TEST_EXIT_DBG(meshLevel + 2 <=
156
157
		  meshDistributor->getMeshLevelData().getNumberOfLevels())
      ("Mesh hierarchy does not contain at least %d levels!\n", meshLevel + 2);
158

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

161
162
    subDomainIsLocal = (levelData.getMpiComm(meshLevel + 1) == MPI::COMM_SELF);

163
    if (subdomain == NULL) {
164
      string subSolverInitStr = name + "->subsolver";
165
166
      string solverType = "petsc";
      Parameters::get(subSolverInitStr, solverType);
167
168
      solverType = "p_" + solverType;
      LinearSolverCreator *solverCreator =
169
	dynamic_cast<LinearSolverCreator*>(CreatorMap<LinearSolverInterface>::getCreator(solverType, name));
170
      TEST_EXIT(solverCreator)
171
	("No valid solver type found in parameter \"%s\"\n",
172
	 name.c_str());
173
174
      solverCreator->setName(subSolverInitStr);
      subdomain = dynamic_cast<PetscSolver*>(solverCreator->create());
175
      subdomain->setSymmetric(isSymmetric);
176
177
      subdomain->setHandleDirichletRows(dirichletMode == 0);
      subdomain->setMeshDistributor(meshDistributor, meshLevel + 1);
178
      subdomain->init(componentSpaces, feSpaces);
179
180

      delete solverCreator;
181
    }
182

183
    primalDofMap.init(componentSpaces, feSpaces);
184
185
186
    dualDofMap.init(componentSpaces, feSpaces, false);
    localDofMap.init(componentSpaces, feSpaces, !subDomainIsLocal);
    lagrangeMap.init(componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188
    if (stokesMode)
189
      interfaceDofMap.init(componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
190

191
    if (fetiPreconditioner == FETI_DIRICHLET) {
192
      TEST_EXIT(levelMode == 1)
193
194
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

195
      interiorDofMap.init(componentSpaces, feSpaces, false);
196
    }
197
198
199
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
200
201
202
203
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
d    
Thomas Witkowski committed
204
    if (dirichletMode == 1) {
205
      int nComponents = mat.getNumRows();
Thomas Witkowski's avatar
d    
Thomas Witkowski committed
206
207
208
209
      for (int component = 0; component < nComponents; component++) {
	DOFMatrix* dofMat = mat[component][component];
	if (!dofMat)
	  continue;
210

Thomas Witkowski's avatar
d    
Thomas Witkowski committed
211
212
	dirichletRows[component] = dofMat->getDirichletRows();
      }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
213
214
215
216
    }
  }


217
218
219
220
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
221
222
    double timeCounter = MPI::Wtime();

223
224
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

225
226
227
228
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
229
    if (fetiPreconditioner == FETI_DIRICHLET)
230
231
      interiorDofMap.clear();

232
233
    primalDofMap.setDofComms(meshDistributor->getDofComms(), meshLevel);
    lagrangeMap.setDofComms(meshDistributor->getDofComms(), meshLevel);
234

235
236
237
238
    primalDofMap.setMpiComm(levelData.getMpiComm(meshLevel));
    dualDofMap.setMpiComm(levelData.getMpiComm(meshLevel));
    lagrangeMap.setMpiComm(levelData.getMpiComm(meshLevel));
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel + 1));
239
    if (fetiPreconditioner == FETI_DIRICHLET)
240
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel + 1));
241

242
    localDofMap.setDofComms(meshDistributor->getDofComms(), meshLevel + 1);
243

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    if (stokesMode) {
245
      interfaceDofMap.clear();
246
      interfaceDofMap.setDofComms(meshDistributor->getDofComms(), meshLevel);
247
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0));
248
249
    }

250
251
    int nComponents = componentSpaces.size();
    for (int component = 0; component < nComponents; component++) {
252
      createPrimals(component);
253
254
255
      createDuals(component);
      createInterfaceNodes(component);
      createIndexB(component);
256
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
257
258
259
260

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

262
    if (fetiPreconditioner == FETI_DIRICHLET)
263
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
264

Thomas Witkowski's avatar
Thomas Witkowski committed
265
    if (stokesMode)
266
267
      interfaceDofMap.update();

268
269
270
    for (int component = 0; component < nComponents; component++) {
      createLagrange(component);
      createAugmentedLagrange(component);
271
    }
272

273
274
275
    lagrangeMap.update();


276
277
    // === ===

278
279
280
    if (subDomainIsLocal) {
      MSG("WARNING: MAKE GENERAL!\n");

281
      rStartInterior = 0;
282
283
284
285
      int localDofs = localDofMap.getOverallDofs();

      mpi::getDofNumbering(domainComm, localDofs,
			   rStartInterior, nGlobalOverallInterior);
286
    } else {
287
288
      MSG("WARNING: MAKE GENERAL!\n");

289
290
291
292
293
294
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

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

295
      mpi::getDofNumbering(domainComm, groupRowsInterior,
296
297
298
299
300
301
			   rStartInterior, nGlobalOverallInterior);

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

302
      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1,
303
					MPI_INT, MPI_SUM);
304
305
    }

306

307
    for (int i = 0; i < static_cast<int>(componentSpaces.size()); i++) {
308
      const FiniteElemSpace *feSpace = componentSpaces[i];
309
310

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

312
      if (i == pressureComponent) {
Thomas Witkowski's avatar
Thomas Witkowski committed
313
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
314
	    interfaceDofMap[i].nRankDofs,
315
	    interfaceDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
316
      } else {
317
318
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n",
	    primalDofMap[i].nRankDofs,
319
320
	    primalDofMap[i].nLocalDofs,
	    primalDofMap[i].nOverallDofs);
321

Thomas Witkowski's avatar
Thomas Witkowski committed
322
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
323
	    dualDofMap[i].nRankDofs,
324
	    dualDofMap[i].nOverallDofs);
325

Thomas Witkowski's avatar
Thomas Witkowski committed
326
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
327
	    lagrangeMap[i].nRankDofs,
328
	    lagrangeMap[i].nOverallDofs);
329
330
331

	MSG("  nRankLocal = %d  nOverallLocal = %d\n",
	    localDofMap[i].nRankDofs,
332
	    localDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
333
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
334
    }
335

Thomas Witkowski's avatar
Thomas Witkowski committed
336
    subdomain->setDofMapping(&localDofMap);
337
    subdomain->setCoarseSpaceDofMapping(&primalDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
340
341

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
342
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
343
      timeCounter = MPI::Wtime() - timeCounter;
344
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
	  timeCounter);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
349
350
351
352


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


356
  void PetscSolverFeti::createPrimals(int component)
Thomas Witkowski's avatar
Thomas Witkowski committed
357
  {
358
    FUNCNAME("PetscSolverFeti::createPrimals()");
359

360
    if (component == pressureComponent)
361
362
      return;

363
364
    const FiniteElemSpace *feSpace = componentSpaces[component];

365
366
367
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
368
    // Set of DOF indices that are considered to be primal variables.
369
    DofContainerSet& vertices =
370
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
371
372

    DofIndexSet primals;
373
    for (DofContainerSet::iterator it = vertices.begin();
374
	 it != vertices.end(); ++it) {
375

376
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
377
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
378

379
380
381
      if (meshLevel == 1 && not (*interiorMap)[component].isSet(**it))
	continue;

382
      if (subDomainIsLocal) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
383
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
384
      } else {
385
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

389
390
391
392
	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
393
394
395
	    (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);
396
397
	} else {
	  MSG("OMMIT SOME PRIMAL!\n");
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
398
399
	}
      }
400
    }
401

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

405
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it) {
406
      if (dofMap[feSpace].isRankDof(*it)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
407
	primalDofMap[component].insertRankDof(*it);
408
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
409
  	primalDofMap[component].insertNonRankDof(*it);
410
411
      }
    }
412
413
414
  }


415
  void PetscSolverFeti::createDuals(int component)
416
417
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
418

419
    if (component == pressureComponent)
420
421
      return;

422
423
    const FiniteElemSpace *feSpace = componentSpaces[component];

424
425
426
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
427
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
428

429
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
430
	 it != allBoundaryDofs.end(); ++it) {
431
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
432
433
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
434
      if (isPrimal(component, **it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
435
	continue;
436
437
438

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

440
      if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
441
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
442
    }
443
444
  }

445

446
  void PetscSolverFeti::createInterfaceNodes(int component)
447
448
449
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

450
    if (component != pressureComponent)
451
452
      return;

453
    const FiniteElemSpace *feSpace = componentSpaces[component];
454
455
456
457
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
458
	 it != allBoundaryDofs.end(); ++it) {
459
      if (dirichletRows[component].count(**it))
460
461
	continue;

462
      if (dofMap[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
463
	interfaceDofMap[component].insertRankDof(**it);
464
      else
465
	interfaceDofMap[component].insertNonRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
466
    }
467
468
469
  }


470
  void PetscSolverFeti::createLagrange(int component)
471
472
473
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

474
    if (component == pressureComponent)
475
476
      return;

477
    const FiniteElemSpace *feSpace = componentSpaces[component];
478
    Mesh* mesh = feSpace->getMesh();
479
480
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
481
482
483
    // 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.
484
    map<int, std::set<DegreeOfFreedom> > subDomainRankDofs;
485

486
487
    if (not subDomainIsLocal) {
      StdMpi<vector<int> > stdMpi(domainComm);
488

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

492
493
	vector<int> dofs;
	dofs.reserve(it.getDofs().size());
494
495

	for (; !it.endDofIter(); it.nextDof()) {
496
	  if (dofMapSubDomain[feSpace].isRankDof(it.getDofIndex()))
497
	    dofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
498
	  else
499
	    dofs.push_back(0);
500
501
	}

502
	stdMpi.send(it.getRank(), dofs);
503
      }
504

505
      for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getSendDofs(), feSpace);
506
507
508
509
510
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

511
      for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getSendDofs(), feSpace);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
512
513
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
514
515
516
	  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
517
    }
518

Thomas Witkowski's avatar
Thomas Witkowski committed
519
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
520
521
522
523
524
525
      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)).         ===

526
    int mpiRank = domainComm.Get_rank();
527
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getSendDofs(), feSpace);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
528
529
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
530
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
531
532
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

533
534
	  // If the subdomain is local, always add the counterpart rank,
	  // otherwise check if the other rank is the owner of the DOF in
535
	  // its subdomain.
536
 	  if (subDomainIsLocal ||
537
	      subDomainRankDofs[it.getRank()].count(it.getDofIndex()))
538
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
539
540
541
542
	}
      }
    }

543
544
545
546

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

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

549
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getSendDofs(), feSpace);
550
551
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
552
	if (!isPrimal(component, it.getDofIndex()))
553
 	  if (subDomainIsLocal || subDomainRankDofs[it.getRank()].count(it.getDofIndex()))
554
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
555
556
557

    stdMpi.updateSendDataSize();

558
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getRecvDofs(), feSpace);
559
	 !it.end(); it.nextRank()) {
560
      bool recvFromRank = false;
561
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
562
	if (!isPrimal(component, it.getDofIndex())) {
563
 	  if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
564
565
566
	    recvFromRank = true;
	    break;
	  }
567
	}
568
      }
569
570

      if (recvFromRank)
571
	stdMpi.recv(it.getRank());
572
    }
573

574
575
    stdMpi.startCommunication();

576
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, meshLevel).getRecvDofs(), feSpace);
577
	 !it.end(); it.nextRank()) {
578
      int i = 0;
579
580
581
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(component, it.getDofIndex())) {
 	  if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(it.getDofIndex())) {
582
	    boundaryDofRanks[feSpace][it.getDofIndex()] =
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
583
	      stdMpi.getRecvData(it.getRank())[i++];
584
	  } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
585
	    lagrangeMap[component].insertNonRankDof(it.getDofIndex());
586
	  }
587
588
	}
      }
589
590
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
591

592
593
594
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

595
    int nRankLagrange = 0;
596
    DofMap& dualMap = dualDofMap[component].getMap();
597
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
598
      if (dofMap[feSpace].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
599
	lagrangeMap[component].insertRankDof(it->first, nRankLagrange);
600
	int degree = boundaryDofRanks[feSpace][it->first].size();
601
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
602
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
603
	lagrangeMap[component].insertNonRankDof(it->first);
604
605
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
606
    lagrangeMap[component].nRankDofs = nRankLagrange;
607
608
609
  }


610
  void PetscSolverFeti::createAugmentedLagrange(int component)
611
612
613
614
615
616
617
618
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


619
  void PetscSolverFeti::createIndexB(int component)
620
  {
621
    FUNCNAME("PetscSolverFeti::createIndexB()");
622

623
    const FiniteElemSpace *feSpace = componentSpaces[component];
624
    DOFAdmin* admin = feSpace->getAdmin();
625
626
627
628

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

630
    int nLocalInterior = 0;
631
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
632
      if (admin->isDofFree(i) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
633
634
635
	  isPrimal(component, i) ||
	  isDual(component, i) ||
	  isInterface(component, i) ||
636
	  dirichletRows[component].count(i))
637
	continue;
638

639
640
641
      if (meshLevel == 1 &&  not (*interiorMap)[component].isSet(i))
	continue;

642
      if (subDomainIsLocal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
643
	localDofMap[component].insertRankDof(i, nLocalInterior);
644

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
645
	if (fetiPreconditioner == FETI_DIRICHLET)
646
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);
647
648

	nLocalInterior++;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
649
      } else {
650
	if (dofMapSubDomain[feSpace].isRankDof(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
651
	  localDofMap[component].insertRankDof(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
652
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
653
	  localDofMap[component].insertNonRankDof(i);
654

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
655
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
656
	  ("Not yet implemnted!\n");
657
      }
658
    }
659

660
661
    // === And finally, add the global indicies of all dual nodes. ===

662
663
    for (DofMap::iterator it = dualDofMap[component].getMap().begin();
	 it != dualDofMap[component].getMap().end(); ++it) {
664
      if (subDomainIsLocal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
665
	localDofMap[component].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
666
      } else {
667
	if (dofMapSubDomain[feSpace].isRankDof(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
668
	  localDofMap[component].insertRankDof(it->first);
669
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
670
	  localDofMap[component].insertNonRankDof(it->first);
671
      }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
672
    }
673
674
675
  }


676
  void PetscSolverFeti::createMatLagrange()
677
678
679
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
680
    double wtime = MPI::Wtime();
681
    int mpiRank = domainComm.Get_rank();
Thomas Witkowski's avatar
Thomas Witkowski committed
682

683
684
    // === Create distributed matrix for Lagrange constraints. ===

685
    MatCreateAIJ(domainComm,
686
687
688
689
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
690
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
691

692

693
    Vec vec_scale_lagrange;
694
    createVec(lagrangeMap, vec_scale_lagrange);
695

696

697
698
699
700
701
702
703
    // === 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.                                                     ===

704
    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
705
      DofMap &dualMap = dualDofMap[component].getMap();
706

707
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
708
	TEST_EXIT_DBG(boundaryDofRanks[componentSpaces[component]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
709
	  ("Should not happen!\n");
710

Thomas Witkowski's avatar
Thomas Witkowski committed
711
	// Global index of the first Lagrange constriant for this node.
712
	int index = lagrangeMap.getMatIndex(component, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
713

Thomas Witkowski's avatar
Thomas Witkowski committed
714
	// Copy set of all ranks that contain this dual node.
715
	vector<int> W(boundaryDofRanks[componentSpaces[component]][it->first].begin(),
716
		      boundaryDofRanks[componentSpaces[component]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
717
718
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
719
720

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
721

722
	int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
723
724
725
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
726
727
	      MatSetValue(mat_lagrange,
			  index + counter,
728
			  localDofMap.getMatIndex(component, it->first) + rStartInterior,
729
730
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
731
	    }
732
733
734
735
736
737
	    counter++;
	  }
	}

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
738

739
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
740
741
742
743
744
745
	  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);
746
747
748
749
750
751
752
	  }
	}
      }
    }

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

754

755
#ifndef NDEBUG