PetscSolverFeti.cc 63.7 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.


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.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 31
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
32
      schurPrimalSolver(0),
33
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
34
      subdomain(NULL),
35
      massMatrixSolver(NULL),
36 37
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
38
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
39
      printTimings(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
40 41
      stokesMode(false),
      pressureComponent(-1)
42 43 44 45 46 47
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

    Parameters::get("parallel->feti->schur primal solver", schurPrimalSolver);
    TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
      ("Wrong solver \"%d\"for the Schur primal complement!\n", 
       schurPrimalSolver);
65 66

    Parameters::get("parallel->multi level test", multiLevelTest);
67 68
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
69 70

    Parameters::get("parallel->print timings", printTimings);
71 72 73
  }


74
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
75
  {
76 77
    FUNCNAME("PetscSolverFeti::initialize()");

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

81
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
82
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
83

Thomas Witkowski's avatar
Thomas Witkowski committed
84

Thomas Witkowski's avatar
Thomas Witkowski committed
85 86 87 88 89 90 91 92
    stokesMode = false;
    Parameters::get("parallel->feti->stokes mode", stokesMode);
    if (stokesMode) {
      Parameters::get("parallel->feti->pressure component", pressureComponent);
      TEST_EXIT(pressureComponent >= 0)
	("FETI-DP in Stokes mode, no pressure component defined!\n");

      pressureFeSpace = feSpaces[pressureComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
93 94
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
95 96
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
97

98
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
99
	subdomain->setMeshDistributor(meshDistributor, 
100
				      mpiCommGlobal, mpiCommLocal);
101
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
102
	subdomain->setMeshDistributor(meshDistributor, 
103 104
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
105
	subdomain->setLevel(meshLevel);
106 107
      }
    }
108

109 110 111 112
    primalDofMap.init(levelData, feSpaces, uniqueFe);   
    dualDofMap.init(levelData, feSpaces, uniqueFe, false);
    localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0);
    lagrangeMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
113

Thomas Witkowski's avatar
Thomas Witkowski committed
114
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
115
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
116

117
    if (fetiPreconditioner == FETI_DIRICHLET) {
118 119 120
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

121
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
122
    }
123 124 125
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
126 127 128 129
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
130 131 132 133 134
    bool removeDirichletRows = false;
    Parameters::get("parallel->feti->remove dirichlet", removeDirichletRows);
    if (!removeDirichletRows)
      return;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151
    int nComponents = mat.getSize();
    for (int i = 0; i < nComponents; i++) {
      DOFMatrix* dofMat = mat[i][i];
      if (!dofMat)
	continue;
      
      const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
      std::set<DegreeOfFreedom>& dRows = dofMat->getDirichletRows();
      if (dirichletRows.count(feSpace)) {
	WARNING("Implement test that for all components dirichlet rows are equal!\n");
      } else {
	dirichletRows[feSpace] = dRows;
      }
    }
  }


152 153 154 155
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
156 157
    double timeCounter = MPI::Wtime();

158 159 160
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
161

162 163
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

164 165 166 167
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
168
    if (fetiPreconditioner == FETI_DIRICHLET)
169 170
      interiorDofMap.clear();

171 172
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
173

174 175
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
176
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
177
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
178
    if (fetiPreconditioner == FETI_DIRICHLET)
179
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
180

181 182 183
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
184 185
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
186
    if (stokesMode) {
187 188 189 190 191
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

192 193
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
194
    
195 196
      createPrimals(feSpace);  

197 198 199
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
200

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
201
      createIndexB(feSpace);     
202
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
203 204 205 206

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
207
    if (fetiPreconditioner == FETI_DIRICHLET)
208
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
209

Thomas Witkowski's avatar
Thomas Witkowski committed
210
    if (stokesMode)
211 212
      interfaceDofMap.update();

213 214 215 216
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
217

218 219 220
    lagrangeMap.update();


221 222 223 224 225 226 227 228 229 230 231 232
    // === ===

    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();

233
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
234 235 236 237 238 239 240 241 242
			   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);
    }

243
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
Thomas Witkowski's avatar
Thomas Witkowski committed
244 245
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
246
      MSG("FETI-DP data for %d-ith FE space: %p\n", i, feSpace);
Thomas Witkowski's avatar
Thomas Witkowski committed
247

Thomas Witkowski's avatar
Thomas Witkowski committed
248
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
	MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
	    primalDofMap[feSpace].nRankDofs, 
	    primalDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
	    dualDofMap[feSpace].nRankDofs, 
	    dualDofMap[feSpace].nOverallDofs);
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
	    lagrangeMap[feSpace].nRankDofs, 
	    lagrangeMap[feSpace].nOverallDofs);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
265 266 267
// 	TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
// 		      static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
// 	  ("Should not happen!\n");	
Thomas Witkowski's avatar
Thomas Witkowski committed
268
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    }
270

Thomas Witkowski's avatar
Thomas Witkowski committed
271
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
272 273 274
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
275 276

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
277
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
278
      timeCounter = MPI::Wtime() - timeCounter;
279
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
280 281
	  timeCounter);
    }
282 283 284
  }


285
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
286
  {
287
    FUNCNAME("PetscSolverFeti::createPrimals()");  
288

Thomas Witkowski's avatar
Thomas Witkowski committed
289
    if (feSpace == pressureFeSpace)
290 291
      return;

292 293 294
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
295
    // Set of DOF indices that are considered to be primal variables.
296
    DofContainerSet& vertices = 
297
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
298 299

    DofIndexSet primals;
300
    for (DofContainerSet::iterator it = vertices.begin(); 
301
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
302 303
      if (dirichletRows[feSpace].count(**it))
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
304

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
305 306
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
307
      } else {
308
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
309 310 311
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
312 313 314 315 316 317 318
	if (fabs(c[0]) < e || fabs(c[1]) < e ||
	    fabs(c[0] - 25.0) < e || fabs(c[1] - 25.0) < e ||
	    (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);
	}
      }
319
    }
320

321 322 323 324

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

325
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
326 327
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
328
      else
329
  	primalDofMap[feSpace].insertNonRankDof(*it);
330 331 332
  }


333
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
334 335
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
336

Thomas Witkowski's avatar
Thomas Witkowski committed
337
    if (feSpace == pressureFeSpace)
338 339
      return;

340 341 342
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
343
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
344
    
345
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
346 347 348 349 350 351 352 353 354 355 356
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;

      if (isPrimal(feSpace, **it))
	continue;

      if (meshLevel == 0) {
	dualDofMap[feSpace].insertRankDof(**it);
      } else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
357
	  dualDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
358 359
      }	  
    }
360 361 362
  }

  
363 364 365 366
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
367
    if (feSpace != pressureFeSpace)
368 369 370 371 372 373
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
374 375 376 377
	 it != allBoundaryDofs.end(); ++it) {
      if (dirichletRows[feSpace].count(**it))
	continue;      
      
378
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
379
	interfaceDofMap[feSpace].insertRankDof(**it);
380
      else
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
381 382
	interfaceDofMap[feSpace].insertNonRankDof(**it);      
    }
383 384 385
  }


386 387 388 389
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
390
    if (feSpace == pressureFeSpace)
391 392
      return;

393 394
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
395 396 397
    // 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
398
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
399

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
400
    if (meshLevel > 0) {
401
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
402 403 404 405 406 407 408 409 410

      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()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
411
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
412
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
413
	  else
414 415 416 417 418 419 420 421 422 423 424 425 426 427 428
	    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
429 430 431 432 433 434
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex()))
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
435

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
436 437 438 439 440 441 442
    if (dualDofMap[feSpace].nLocalDofs == 0)
      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)).         ===

443
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
444 445 446 447 448 449 450
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(feSpace, it.getDofIndex())) {
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

451
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
452
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
453
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
454 455 456 457
	}
      }
    }

458 459 460 461

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

462
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
463

464
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
465 466
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
467
	if (!isPrimal(feSpace, it.getDofIndex()))
468 469 470
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
471 472 473

    stdMpi.updateSendDataSize();

474
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
475
	 !it.end(); it.nextRank()) {
476
      bool recvFromRank = false;
477
      for (; !it.endDofIter(); it.nextDof()) {
478
	if (!isPrimal(feSpace, it.getDofIndex())) {
479 480 481
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
482 483 484
	    recvFromRank = true;
	    break;
	  }
485
	}
486
      }
487 488

      if (recvFromRank)
489
	stdMpi.recv(it.getRank());
490
    }
491

492 493
    stdMpi.startCommunication();

494
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
495
	 !it.end(); it.nextRank()) {
496
      int i = 0;
497
      for (; !it.endDofIter(); it.nextDof())
498
	if (!isPrimal(feSpace, it.getDofIndex()))
499 500 501
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
502 503
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
504 505
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
506 507
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
508

509 510 511
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

512
    int nRankLagrange = 0;
513
    DofMap& dualMap = dualDofMap[feSpace].getMap();
514
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
515 516
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
517
	int degree = boundaryDofRanks[feSpace][it->first].size();
518
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
519
      } else {
520
	lagrangeMap[feSpace].insertNonRankDof(it->first);
521 522
      }
    }
523
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
524 525 526
  }


527
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
528
  {
529
    FUNCNAME("PetscSolverFeti::createIndexB()");
530

531
    DOFAdmin* admin = feSpace->getAdmin();
532 533 534 535

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

537
    int nLocalInterior = 0;    
538
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
539 540 541 542 543 544
      if (admin->isDofFree(i) ||
	  isPrimal(feSpace, i) ||
	  isDual(feSpace, i) ||
	  isInterface(feSpace, i) ||
	  dirichletRows[feSpace].count(i))
	continue;      
545

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
546 547 548 549 550 551 552 553 554 555 556 557 558 559 560
      if (meshLevel == 0) {
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);
	
	if (fetiPreconditioner == FETI_DIRICHLET)
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	
	nLocalInterior++;	
      } else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	  localDofMap[feSpace].insertRankDof(i);
	else
	  localDofMap[feSpace].insertNonRankDof(i);
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
561
      }
562
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
563
    
564 565
    // === And finally, add the global indicies of all dual nodes. ===

566
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
567 568
	 it != dualDofMap[feSpace].getMap().end(); ++it) {
      if (meshLevel == 0) {
569
	localDofMap[feSpace].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
570
      } else {
571 572 573 574 575
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
576
    }
577 578 579
  }


580
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
581 582 583
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
584
    double wtime = MPI::Wtime();
585
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
586

587 588
    // === Create distributed matrix for Lagrange constraints. ===

589 590 591 592 593
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
594
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
595

596 597 598 599 600 601 602
    // === 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.                                                     ===

Thomas Witkowski's avatar
Thomas Witkowski committed
603
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
604
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
605

606
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
607
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
608 609 610
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
611
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
612

Thomas Witkowski's avatar
Thomas Witkowski committed
613
	// Copy set of all ranks that contain this dual node.
614 615
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
616 617
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
618 619

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
620 621 622 623
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
624 625
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
626

627
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
628
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
629
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
630
	    index++;	      
631 632 633 634 635 636 637
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
639 640
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
641 642
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
643
    }
644 645 646
  }


647
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
648
  {
649
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
650

Thomas Witkowski's avatar
Thomas Witkowski committed
651
    if (schurPrimalSolver == 0) {
652 653
      MSG("Create iterative schur primal solver!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
654
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
655
      
656 657
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
658

659
      MatCreateShell(mpiCommGlobal,
660 661 662 663
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
664 665 666 667 668
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
669
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
670
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
671 672
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
673 674
      KSPSetFromOptions(ksp_schur_primal);
    } else {
675 676
      MSG("Create direct schur primal solver!\n");

677 678
      double wtime = MPI::Wtime();

679 680 681
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

682 683 684
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
685

Thomas Witkowski's avatar
Thomas Witkowski committed
686
      Mat matBPi;
687 688 689
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
690 691 692
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

693

694 695 696
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
697

698 699
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

700
      for (int i = 0; i < nRowsRankB; i++) {
701
	PetscInt row = localDofMap.getStartDofs() + i;
702
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
703 704 705 706 707

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

708
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
709 710
      }

711
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
712
		primalDofMap.getLocalDofs())
713
	("Should not happen!\n");
714

715 716 717
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
718
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
719 720 721 722 723 724 725

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
726

Thomas Witkowski's avatar
Thomas Witkowski committed
727
	subdomain->solve(tmpVec, tmpVec);
728

Thomas Witkowski's avatar
Thomas Witkowski committed
729
	PetscScalar *tmpValues;
730
	VecGetArray(tmpVec, &tmpValues);
731
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
732
	  MatSetValue(matBPi, 
733
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
734 735 736
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
737 738 739 740 741
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
742 743
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
744

745 746
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
747 748

      Mat matPrimal;
749 750
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
751
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
752 753

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
754
      MatDestroy(&matBPi);
755

756
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
757
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
758
		      SAME_NONZERO_PATTERN);
759 760 761 762 763 764
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPPREONLY);
      PC pc_schur_primal;      
      KSPGetPC(ksp_schur_primal, &pc_schur_primal);
      PCSetType(pc_schur_primal, PCLU);
      PCFactorSetMatSolverPackage(pc_schur_primal, MATSOLVERMUMPS);
765
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
766

Thomas Witkowski's avatar
Thomas Witkowski committed
767
      if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
768
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
769 770 771 772 773 774 775 776 777 778
	MatInfo minfo;
	MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
	MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);
	
	MSG("FETI-DP timing 06: %.5f seconds (creation of schur primal matrix)\n",
	    MPI::Wtime() - wtime);

	wtime = MPI::Wtime();
	KSPSetUp(ksp_schur_primal);
	KSPSetUpOnBlocks(ksp_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
779
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
780 781 782
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
783
    }
784 785 786 787 788 789 790
  }


  void PetscSolverFeti::destroySchurPrimalKsp()
  {
    FUNCNAME("PetscSolverFeti::destroySchurPrimal()");

791
    if (schurPrimalSolver == 0) {
792
      schurPrimalData.subSolver = NULL;
793

794 795 796
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
797 798 799
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
800 801 802
  }


803
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
804 805 806
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

807 808
    // === Create FETI-DP solver object. ===

809 810 811 812
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
813
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
814
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
815
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
816

Thomas Witkowski's avatar
Thomas Witkowski committed
817
    if (stokesMode == false) {      
818 819 820 821 822 823 824 825
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
826
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
827 828 829 830
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

831 832 833 834 835
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
836
		     &fetiData, &mat_feti);
837
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
838
			   (void(*)(void))petscMultMatFetiInterface);      
839
    }
840

841
    KSPCreate(mpiCommGlobal, &ksp_feti);
842
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
843 844 845
    KSPSetOptionsPrefix(ks