PetscSolverFeti.cc 87.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * 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.
 * 
 ******************************************************************************/
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

35
namespace AMDiS { namespace Parallel {
36
37
38

  using namespace std;

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

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

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

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

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

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

103
104
105
106
    {
      MSG("WARNING: CHECK THIS HERE BEFORE GOING INTO RUNNING MULTILEVEL FETI-DP!\n");
      Parameters::get("parallel->level mode", levelMode);
    }
107
108
109
    
    {
      int tmp = 0;
110
      Parameters::get(name + "->feti->inexact", tmp);
111
112
113
114
115
      if (tmp == 1)
	fetiSolverType = INEXACT;
      if (tmp == 2)
	fetiSolverType = INEXACT_REDUCED;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
116
117

    Parameters::get("parallel->print timings", printTimings);
118
  }
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
  

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

    super::init(fe0, fe1, createGlobalMapping);
    
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
    int nLevels = levelData.getNumberOfLevels();
    TEST_EXIT_DBG(nLevels >= 1)("nLevels < 1! Should not happen!\n");

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


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

150
151
152
#if (DEBUG != 0)
    MSG("Init FETI-DP on mesh level %d\n", meshLevel);
#endif
153

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

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

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

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

      delete solverCreator;
180
    }
181

182
183
184
185
    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
186

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

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

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


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

Thomas Witkowski's avatar
d  
Thomas Witkowski committed
203
    if (dirichletMode == 1) {
204
      int nComponents = mat.getNumRows();
Thomas Witkowski's avatar
d  
Thomas Witkowski committed
205
206
207
208
209
210
211
      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
212
213
214
215
    }
  }


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

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

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

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

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

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

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

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

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

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

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

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

272
273
274
    lagrangeMap.update();


275
276
    // === ===

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

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

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

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

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

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

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

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

305

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

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

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

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

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

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

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


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


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

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

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

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

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

    DofIndexSet primals;
372
    for (DofContainerSet::iterator it = vertices.begin(); 
373
	 it != vertices.end(); ++it) {
374
      
375
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
376
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
377

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

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

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

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

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


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

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

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

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

    DofContainer allBoundaryDofs;
426
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
427
    
428
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
429
	 it != allBoundaryDofs.end(); ++it) {
430
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
431
432
	continue;

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

      if (meshLevel == 1 && not (*interiorMap)[component].isSet(**it))
	continue;
438
      
439
      if (subDomainIsLocal || dofMapSubDomain[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
440
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
441
    }
442
443
444
  }

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

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

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
457
	 it != allBoundaryDofs.end(); ++it) {
458
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
459
460
	continue;      
      
461
      if (dofMap[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
462
	interfaceDofMap[component].insertRankDof(**it);
463
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
464
	interfaceDofMap[component].insertNonRankDof(**it);      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
465
    }
466
467
468
  }


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

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

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

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

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

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

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

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

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

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

      stdMpi.startCommunication();

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

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

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

532
533
534
535
536
	  // 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()))
537
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
538
539
540
541
	}
      }
    }

542
543
544
545

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

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

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

    stdMpi.updateSendDataSize();

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

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

573
574
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
590

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

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


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

    if (!augmentedLagrange)
      return;
  }


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

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

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

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

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

641
      if (subDomainIsLocal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
642
	localDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
643
644
	
	if (fetiPreconditioner == FETI_DIRICHLET)
645
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
646
647
648
	
	nLocalInterior++;	
      } else {
649
	if (dofMapSubDomain[feSpace].isRankDof(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
650
	  localDofMap[component].insertRankDof(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
651
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
652
	  localDofMap[component].insertNonRankDof(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
653
654
655
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
656
      }
657
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
658
    
659
660
    // === And finally, add the global indicies of all dual nodes. ===

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


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

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

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

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

691

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

695

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

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

706
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
707
	TEST_EXIT_DBG(boundaryDofRanks[componentSpaces[component]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
708
709
710
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
711
	int index = lagrangeMap.getMatIndex(component, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
712

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

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

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
738
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
739
740
741
742
743
744
	  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);
745
746
747
748
749
750
751
	  }
	}
      }
    }

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

753

754
#if (DEBUG != 0)
755
756
757
758
    {
      int nZeroRows = PetscSolverFetiDebug::testZeroRows(mat_lagrange);
      int m,n;
      MatGetSize(mat_lagrange, &m ,&n);
759
      mpi::globalAdd(domainComm, nZeroRows);
760
      MSG("Lagrange matrix has %d zero rows and global size of %d %d!\n", nZeroRows, m, n);
761
762

      TEST_EXIT(nZeroRows == 0)("Lagrange matrix has zero rows!\n");
763
    }
764
#endif
Thomas Witkowski's avatar
So on  
Thomas Witkowski committed
765
766


767
768
769
770
771
    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);