PetscSolverFeti.cc 80.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include "MatrixVector.h"
14
#include "parallel/PetscHelper.h"
15
#include "parallel/PetscSolverFeti.h"
16
17
#include "parallel/PetscSolverFetiDebug.h"
#include "parallel/PetscSolverFetiMonitor.h"
18
#include "parallel/PetscSolverFetiStructs.h"
19
20
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
21
22
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
23
#include "parallel/PetscSolverGlobalMatrix.h"
24
#include "io/VtkWriter.h"
25
26
27
28
29

namespace AMDiS {

  using namespace std;

30
  PetscSolverFeti::PetscSolverFeti(string name)
31
    : PetscSolver(name),
32
33
34
35
36
37
      primalDofMap(COMPONENT_WISE),
      dualDofMap(COMPONENT_WISE),
      interfaceDofMap(COMPONENT_WISE),
      localDofMap(COMPONENT_WISE),
      lagrangeMap(COMPONENT_WISE),
      interiorDofMap(COMPONENT_WISE),
38
      schurPrimalSolver(0),
39
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
40
      subdomain(NULL),
41
      massMatrixSolver(NULL),
42
43
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
44
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
45
      printTimings(false),
Thomas Witkowski's avatar
d    
Thomas Witkowski committed
46
      dirichletMode(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
47
      stokesMode(false),
48
      augmentedLagrange(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
49
      pressureComponent(-1)
50
51
52
53
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

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

92
    Parameters::get("parallel->multi level test", multiLevelTest);
93
94
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
95
96

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


100
  void PetscSolverFeti::initialize()
101
  {
102
103
    FUNCNAME("PetscSolverFeti::initialize()");

104
105
    MSG("INIT WITH MESH-LEVEL: %d\n", meshLevel);

106
107
108
    TEST_EXIT_DBG(meshLevel + 1 == meshDistributor->getMeshLevelData().getLevelNumber())
      ("Mesh hierarchy does not contain %d levels!\n", meshLevel + 1);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
111
    if (subdomain == NULL) {
112
113
114
115
116
117
118
119
120
121
122
      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();
123
      subdomain->setSymmetric(isSymmetric);
Thomas Witkowski's avatar
d    
Thomas Witkowski committed
124
125
126
127
      if (dirichletMode == 0)
	subdomain->setHandleDirichletRows(true);
      else
	subdomain->setHandleDirichletRows(false);
128

129
      if (meshLevel == 0) {
130
	subdomain->setMeshDistributor(meshDistributor);
131
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
132
	subdomain->setMeshDistributor(meshDistributor, 
133
134
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
135
	subdomain->setLevel(meshLevel);
136
      }
137
138

      delete solverCreator;
139
    }
140

141
142
143
144
    primalDofMap.init(levelData, componentSpaces, feSpaces);   
    dualDofMap.init(levelData, componentSpaces, feSpaces, false);
    localDofMap.init(levelData, componentSpaces, feSpaces, meshLevel != 0);
    lagrangeMap.init(levelData, componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
145

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

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

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


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

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


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

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

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

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

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

193
194
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
195
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
196
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
197
    if (fetiPreconditioner == FETI_DIRICHLET)
198
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
199

200
201
202
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
203
204
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
205
    if (stokesMode) {
206
207
208
209
210
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

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

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
222
    if (fetiPreconditioner == FETI_DIRICHLET)
223
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
224

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

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

233
234
235
    lagrangeMap.update();


236
237
238
239
240
241
242
243
244
245
246
247
    // === ===

    if (meshLevel == 0) {
      rStartInterior = 0;
      nGlobalOverallInterior = localDofMap.getOverallDofs();
    } else {
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

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

248
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
249
250
251
252
253
254
255
256
257
			   rStartInterior, nGlobalOverallInterior);

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

      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }

258
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
259
260
    for (unsigned int i = 0; i < componentSpaces.size(); i++) {
      const FiniteElemSpace *feSpace = componentSpaces[i];
261
262

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

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

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

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


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


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

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

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

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

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

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

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

334
335
336
337
	if ((fabs(c[0]) < e && fabs(c[1] - 12.5) < e) ||
	    (fabs(c[0] - 25.0) < e && fabs(c[1] - 12.5) < e) ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1]) < e) ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 25.0) < e) ||
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
338
339
340
341
342
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 12.5) < e)) {
	  MSG("PRIMAL COORD %f %f\n", c[0], c[1]);
	  primals.insert(**it);
	}
      }
343
    }
344

345
346
347
348

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

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


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

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

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

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

    DofContainer allBoundaryDofs;
371
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
372
    
373
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
374
	 it != allBoundaryDofs.end(); ++it) {
375
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
376
377
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
378
      if (isPrimal(component, **it))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
379
380
381
	continue;

      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
382
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
383
      } else {
384
	if (dofMapSd[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
385
	  dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
386
387
      }	  
    }
388
389
390
  }

  
391
  void PetscSolverFeti::createInterfaceNodes(int component)
392
393
394
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

395
    if (component != pressureComponent)
396
397
      return;

398
    const FiniteElemSpace *feSpace = componentSpaces[component];
399
400
401
402
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

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


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

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

422
    const FiniteElemSpace *feSpace = componentSpaces[component];
423
424
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
427
    // 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.
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
428
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
429

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
430
    if (meshLevel > 0) {
431
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
432
433
434
435
436
437
438
439
440

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

	vector<int> subdomainRankDofs;
	subdomainRankDofs.reserve(it.getDofs().size());

	for (; !it.endDofIter(); it.nextDof()) {
441
	  if (dofMapSd[feSpace].isRankDof(it.getDofIndex()))
442
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
443
	  else
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
	    subdomainRankDofs.push_back(0);
	}

	stdMpi.send(it.getRank(), subdomainRankDofs);
      }	     

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace); 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
459
460
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
461
	  if (!isPrimal(component, it.getDofIndex()))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
462
463
464
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
465

Thomas Witkowski's avatar
Thomas Witkowski committed
466
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
467
468
469
470
471
472
      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)).         ===

473
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
474
475
476
477
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
478
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
479
480
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

481
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
482
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
483
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
484
485
486
487
	}
      }
    }

488
489
490
491

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

492
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
493

494
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
495
496
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
497
	if (!isPrimal(component, it.getDofIndex()))
498
499
500
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
501
502
503

    stdMpi.updateSendDataSize();

504
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
505
	 !it.end(); it.nextRank()) {
506
      bool recvFromRank = false;
507
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
508
	if (!isPrimal(component, it.getDofIndex())) {
509
510
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
511
	       dofMapSd[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
512
513
514
	    recvFromRank = true;
	    break;
	  }
515
	}
516
      }
517
518

      if (recvFromRank)
519
	stdMpi.recv(it.getRank());
520
    }
521

522
523
    stdMpi.startCommunication();

524
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
525
	 !it.end(); it.nextRank()) {
526
      int i = 0;
527
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
528
	if (!isPrimal(component, it.getDofIndex()))
529
530
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
531
	       dofMapSd[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
532
533
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
534
	  else
Thomas Witkowski's avatar
Thomas Witkowski committed
535
	    lagrangeMap[component].insertNonRankDof(it.getDofIndex());
536
537
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
538

539
540
541
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


557
  void PetscSolverFeti::createAugmentedLagrange(int component)
558
559
560
561
562
563
564
565
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


566
  void PetscSolverFeti::createIndexB(int component)
567
  {
568
    FUNCNAME("PetscSolverFeti::createIndexB()");
569

570
    const FiniteElemSpace *feSpace = componentSpaces[component];
571
    DOFAdmin* admin = feSpace->getAdmin();
572
573
574
575

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

577
    int nLocalInterior = 0;    
578
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
579
      if (admin->isDofFree(i) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
580
581
582
	  isPrimal(component, i) ||
	  isDual(component, i) ||
	  isInterface(component, i) ||
583
	  dirichletRows[component].count(i))
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
584
	continue;      
585

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
586
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
587
	localDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
588
589
	
	if (fetiPreconditioner == FETI_DIRICHLET)
590
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
591
592
593
	
	nLocalInterior++;	
      } else {
594
	if (dofMapSd[feSpace].isRankDof(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
595
	  localDofMap[component].insertRankDof(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
596
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
597
	  localDofMap[component].insertNonRankDof(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
598
599
600
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
601
      }
602
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
603
    
604
605
    // === And finally, add the global indicies of all dual nodes. ===

606
607
    for (DofMap::iterator it = dualDofMap[component].getMap().begin();
	 it != dualDofMap[component].getMap().end(); ++it) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
608
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
609
	localDofMap[component].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
610
      } else {
611
	if (dofMapSd[feSpace].isRankDof(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
612
	  localDofMap[component].insertRankDof(it->first);
613
	else 
Thomas Witkowski's avatar
Thomas Witkowski committed
614
	  localDofMap[component].insertNonRankDof(it->first);
615
      }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
616
    }
617
618
619
  }


620
  void PetscSolverFeti::createMatLagrange()
621
622
623
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
624
    double wtime = MPI::Wtime();
625
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
626

627
628
    // === Create distributed matrix for Lagrange constraints. ===

629
630
631
632
633
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
634
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
635

636
637
638
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

639
640
641
642
643
644
645
    // === 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.                                                     ===

646
    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
647
      DofMap &dualMap = dualDofMap[component].getMap();
648

649
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
650
	TEST_EXIT_DBG(boundaryDofRanks[componentSpaces[component]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
651
652
653
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
654
	int index = lagrangeMap.getMatIndex(component, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
655

Thomas Witkowski's avatar
Thomas Witkowski committed
656
	// Copy set of all ranks that contain this dual node.
657
658
	vector<int> W(boundaryDofRanks[componentSpaces[component]][it->first].begin(), 
		      boundaryDofRanks[componentSpaces[component]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
659
660
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
661
662

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

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
681
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
682
683
684
685
686
687
	  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);
688
689
690
691
692
693
694
	  }
	}
      }
    }

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

696

Thomas Witkowski's avatar
So on    
Thomas Witkowski committed
697
698
699
700
701
702
703
704
705
706
707
708
    int nZeroRows = PetscSolverFetiDebug::testZeroRows(mat_lagrange);
    int m,n;
    MatGetSize(mat_lagrange, &m ,&n);
    MSG("Lagrange matrix has %d zero rows and global size of %d %d!\n", nZeroRows, m, n);

    PetscViewer petscView;
    PetscViewerBinaryOpen(PETSC_COMM_WORLD, "lagrange.mat", 
			  FILE_MODE_WRITE, &petscView);
    MatView(mat_lagrange, petscView);
    PetscViewerDestroy(&petscView);


709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

    if (fetiPreconditioner != FETI_NONE || stokesMode) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatDiagonalScale(mat_lagrange_scaled, vec_scale_lagrange, PETSC_NULL);
    }

    VecDestroy(&vec_scale_lagrange);


    // === Print final timings. ===

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
724
725
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
726
727
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
728
    }
729
730
  }

731

Thomas Witkowski's avatar
Thomas Witkowski committed
732
733
  bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
  {
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
734
735
    FUNCNAME("PetscSolverFeti::testWirebasketEdge()");

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
736
737
    return true;

Thomas Witkowski's avatar
a    
Thomas Witkowski committed
738
739
740
    if (meshDistributor->getMesh()->getDim() == 2)
      return true;

Thomas Witkowski's avatar
Thomas Witkowski committed
741
742
743
    if (meshDistributor->getIntBoundary().getDegreeOwn(edge) != 3)
      return false;

Thomas Witkowski's avatar
da    
Thomas Witkowski committed
744
745
    return true;

Thomas Witkowski's avatar
Thomas Witkowski committed
746
747
748
749
750
751
752
753
754
755
756
757
758
759
    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);
  }
760

761

762
  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseEdges()
763
  {
764
    FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");
765
766

    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
Thomas Witkowski's avatar
Thomas Witkowski committed
767
    std::set<BoundaryObject> allEdges;
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
    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);
      }
    }
788

789
790
791
    int nEdges = allEdges.size();
    mpi::globalAdd(nEdges);
    mpi::globalAdd(nEmptyEdges);
Thomas Witkowski's avatar
Thomas Witkowski committed
792

793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
    MSG("Coarse space edges: %d (empty: %d)\n", nEdges, nEmptyEdges);
    
    return data;
  }


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

    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
    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
809
#if 0
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
    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
826
827
828
829
830
831
    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
832
833
	      intBound.getDegreeOwn(bs[i].rankObj) == 1 &&
	      allMyEdges.count(bs[i].rankObj)) {
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
834
835
	    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
836
	    it->second.insert(bs[i].rankObj);
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
837
838
839
840
	  }	      
	}	  
      }
    }
Thomas Witkowski's avatar
da    
Thomas Witkowski committed
841
#endif
Thomas Witkowski's avatar
a    
Thomas Witkowski committed
842
843


844
845
846
847
848
849
850
851
852
853
854
    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);
855
      }
856
857
858
859
      if (oneBoundary.size()) 
	data.push_back(oneBoundary);
      else 
	nEmptyFaces++;
860
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
861

862
    int nFaces = allFaces.size();
863
    mpi::globalAdd(nFaces);
864
    mpi::globalAdd(nEmptyFaces);
865

866
867
868
    MSG("Coarse space faces: %d (empty: %d)\n", nFaces, nEmptyFaces);
   
    return data;
869
870
871
872
873
874
875
876
877
878
879
880
  }


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

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

881
882
883
884
885
    vector<vector<BoundaryObject> > allEdges = getCoarseEdges();
    vector<vector<BoundaryObject> > allFaces = getCoarseFaces(); 
    allEdges.insert(allEdges.end(), allFaces.begin(), allFaces.end());

    nRankEdges = allEdges.size();
886
887
888
889
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

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

891
892
893
    nRankEdges *= componentSpaces.size();
    rStartEdges *= componentSpaces.size();
    nOverallEdges *= componentSpaces.size();
894
895

    MatCreateAIJ(mpiCommGlobal,
896
897
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
898
		 2, PETSC_NULL, 2, PETSC_NULL, 
899
900
901
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

902
    int rowCounter = rStartEdges;
903
904
    for (vector<vector<BoundaryObject> >::iterator it = allEdges.begin(); 
	 it != allEdges.end(); ++it) {
905
      for (int component = 0; component < componentSpaces.size(); component++) {
906
907
908
909
910
	for (vector<BoundaryObject>::iterator edgeIt = it->begin();
	     edgeIt != it->end(); ++edgeIt) {
	  
	  DofContainer edgeDofs;
	  edgeIt->el->getAllDofs(componentSpaces[component], *edgeIt, edgeDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
911
	  
912
	  TEST_EXIT(edgeDofs.size())("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
913
	  
914
915
916
917
918
919
920
921
922
	  for (DofContainer::iterator it = edgeDofs.begin();
	       it != edgeDofs.end(); ++it) {
	    TEST_EXIT(isPrimal(component, **it) == false)
	      ("Should not be primal!\n");
	    
	    int col = lagrangeMap.getMatIndex(component, **it);
	    double value = 1.0;
	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
	  }
923
	}