PetscSolverFeti.cc 70.5 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
      stokesMode(false),
41
      augmentedLagrange(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
42
      pressureComponent(-1)
43 44 45 46 47 48
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

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

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

    Parameters::get("parallel->feti->augmented lagrange", augmentedLagrange);
74 75 76
  }


77
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
78
  {
79 80
    FUNCNAME("PetscSolverFeti::initialize()");

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

84
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
85
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
86

Thomas Witkowski's avatar
Thomas Witkowski committed
87

Thomas Witkowski's avatar
Thomas Witkowski committed
88 89 90 91 92 93 94 95
    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
96 97
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
98 99
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
100

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

112 113 114 115
    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
116

Thomas Witkowski's avatar
Thomas Witkowski committed
117
    if (stokesMode)
Thomas Witkowski's avatar
Thomas Witkowski committed
118
      interfaceDofMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
119

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

124
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
125
    }
126 127 128
  }


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

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
    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;
      }
    }
  }


155 156 157 158
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
159 160
    double timeCounter = MPI::Wtime();

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

165 166
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

167 168 169 170
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
171
    if (fetiPreconditioner == FETI_DIRICHLET)
172 173
      interiorDofMap.clear();

174 175
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
176

177 178
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
179
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
180
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
181
    if (fetiPreconditioner == FETI_DIRICHLET)
182
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
183

184 185 186
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
187 188
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

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

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

200 201 202
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
203

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
204
      createIndexB(feSpace);     
205
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
206 207 208 209

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
210
    if (fetiPreconditioner == FETI_DIRICHLET)
211
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
212

Thomas Witkowski's avatar
Thomas Witkowski committed
213
    if (stokesMode)
214 215
      interfaceDofMap.update();

216 217 218
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
219
      createAugmentedLagrange(feSpace);
220
    }
221

222 223 224
    lagrangeMap.update();


225 226 227 228 229 230 231 232 233 234 235 236
    // === ===

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

237
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
238 239 240 241 242 243 244 245 246
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
252
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
253 254 255 256
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
	    interfaceDofMap[feSpace].nRankDofs, 
	    interfaceDofMap[feSpace].nOverallDofs);
      } else {
257
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
258
	    primalDofMap[feSpace].nRankDofs, 
259
	    primalDofMap[feSpace].nLocalDofs,
Thomas Witkowski's avatar
Thomas Witkowski committed
260 261 262 263 264 265 266 267 268 269
	    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
270 271 272
// 	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
273
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
274
    }
275

Thomas Witkowski's avatar
Thomas Witkowski committed
276
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
277 278 279
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
280 281

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


290
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
291
  {
292
    FUNCNAME("PetscSolverFeti::createPrimals()");  
293

Thomas Witkowski's avatar
Thomas Witkowski committed
294
    if (feSpace == pressureFeSpace)
295 296
      return;

297 298 299
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
300
    // Set of DOF indices that are considered to be primal variables.
301
    DofContainerSet& vertices = 
302
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
303 304

    DofIndexSet primals;
305
    for (DofContainerSet::iterator it = vertices.begin(); 
306
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
307 308
      if (dirichletRows[feSpace].count(**it))
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
309

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
310 311
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
312
      } else {
313
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
314 315 316
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
317 318 319 320 321 322 323
	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);
	}
      }
324
    }
325

326 327 328 329

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

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


338
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
339 340
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
341

Thomas Witkowski's avatar
Thomas Witkowski committed
342
    if (feSpace == pressureFeSpace)
343 344
      return;

345 346 347
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
348
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
349
    
350
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
351 352 353 354 355 356 357 358 359 360 361
	 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))
362
	  dualDofMap[feSpace].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
363 364
      }	  
    }
365 366 367
  }

  
368 369 370 371
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
372
    if (feSpace != pressureFeSpace)
373 374 375 376 377 378
      return;

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

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


391 392 393 394
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
395
    if (feSpace == pressureFeSpace)
396 397
      return;

398 399
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
400 401 402
    // 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
403
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
404

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
405
    if (meshLevel > 0) {
406
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
407 408 409 410 411 412 413 414 415

      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
416
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
417
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
418
	  else
419 420 421 422 423 424 425 426 427 428 429 430 431 432 433
	    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
434 435 436 437 438 439
	   !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());
    }
440

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
441 442 443 444 445 446 447
    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)).         ===

448
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
449 450 451 452 453 454 455
    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);

456
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
457
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
458
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
459 460 461 462
	}
      }
    }

463 464 465 466

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

467
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
468

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

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
494
	stdMpi.recv(it.getRank());
495
    }
496

497 498
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
513

514 515 516
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


532 533 534 535 536 537 538 539 540
  void PetscSolverFeti::createAugmentedLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


541
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
542
  {
543
    FUNCNAME("PetscSolverFeti::createIndexB()");
544

545
    DOFAdmin* admin = feSpace->getAdmin();
546 547 548 549

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

551
    int nLocalInterior = 0;    
552
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
553 554 555 556 557 558
      if (admin->isDofFree(i) ||
	  isPrimal(feSpace, i) ||
	  isDual(feSpace, i) ||
	  isInterface(feSpace, i) ||
	  dirichletRows[feSpace].count(i))
	continue;      
559

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
560 561 562 563 564 565 566 567 568 569 570 571 572 573 574
      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");	
575
      }
576
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
577
    
578 579
    // === And finally, add the global indicies of all dual nodes. ===

580
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
581 582
	 it != dualDofMap[feSpace].getMap().end(); ++it) {
      if (meshLevel == 0) {
583
	localDofMap[feSpace].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
584
      } else {
585 586 587 588 589
	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
590
    }
591 592 593
  }


594
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
595 596 597
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
598
    double wtime = MPI::Wtime();
599
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
600

601 602
    // === Create distributed matrix for Lagrange constraints. ===

603 604 605 606 607
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
608
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
609

610 611 612
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

613 614 615 616 617 618 619
    // === 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
620
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
621
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
622

623
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
624
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
625 626 627
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
628
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
629

Thomas Witkowski's avatar
Thomas Witkowski committed
630
	// Copy set of all ranks that contain this dual node.
631 632
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
633 634
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
635 636

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
637
	
638
	int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
639 640 641
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
642 643 644 645 646
	      MatSetValue(mat_lagrange, 
			  index + counter, 
			  localDofMap.getMatIndex(k, it->first) + rStartInterior,
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
647
	    }
648 649 650 651 652 653 654 655 656 657 658 659 660 661
	    counter++;
	  }
	}

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
	if (meshDistributor->getDofMap()[feSpaces[k]].isRankDof(it->first)) {
	  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);
662 663 664 665 666 667 668
	  }
	}
      }
    }

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

670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687

    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

    VecView(vec_scale_lagrange, PETSC_VIEWER_STDOUT_WORLD);

    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
688 689
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
690 691
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
692
    }
693 694 695
  }


696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721
  void PetscSolverFeti::createMatAugmentedLagrange(vector<const FiniteElemSpace*> &feSpaces)
  {
    FUNCNAME("PetscSolverFeti::createMatAugmentedLagrange()");

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

    nRankEdges = 0;
    nOverallEdges = 0;
    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
      if (it->rankObj.subObj == EDGE)
	nRankEdges++;
    
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

    MSG("nRankEdges = %d, nOverallEdges = %d\n", nRankEdges, nOverallEdges);

    nRankEdges *= feSpaces.size();
    rStartEdges *= feSpaces.size();
    nOverallEdges *= feSpaces.size();

    MatCreateAIJ(mpiCommGlobal,
722 723
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
724 725 726 727
		 1, PETSC_NULL, 1, PETSC_NULL, 
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

728
    int rowCounter = rStartEdges;
729 730 731 732 733 734 735 736 737 738 739 740 741
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
      if (it->rankObj.subObj == EDGE) {
	for (int i = 0; i < feSpaces.size(); i++) {
	  DofContainer edgeDofs;
	  it->rankObj.el->getAllDofs(feSpaces[i], it->rankObj, edgeDofs);

	  MSG("SIZE = %d\n", edgeDofs.size());

	  for (DofContainer::iterator it = edgeDofs.begin();
	       it != edgeDofs.end(); it++) {
	    TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false)
	      ("Should not be primal!\n");

742
	    int col = lagrangeMap.getMatIndex(i, **it);
743
	    double value = 1.0;
744
	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
745 746
	  }

747
	  rowCounter++;
748 749 750 751 752 753 754 755 756 757 758 759 760 761 762
	}
      }
    }

    MatAssemblyBegin(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_augmented_lagrange, MAT_FINAL_ASSEMBLY);

    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
      MSG("FETI-DP timing 05a: %.5f seconds (creation of augmented lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
    }
  }


763
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
764
  {
765
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
766

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

770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805
      if (augmentedLagrange == false) {
	schurPrimalData.subSolver = subdomain;
	
	localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
	
	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getRankDofs(), 
		       primalDofMap.getOverallDofs(), 
		       primalDofMap.getOverallDofs(),
		       &schurPrimalData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimal);	
      } else {
	schurPrimalAugmentedData.subSolver = subdomain;

	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b0, nGlobalOverallInterior);
	localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b1, nGlobalOverallInterior);
	primalDofMap.createVec(schurPrimalAugmentedData.tmp_vec_primal);
	lagrangeMap.createVec(schurPrimalAugmentedData.tmp_vec_lagrange);

	schurPrimalAugmentedData.mat_lagrange = &mat_lagrange;
	schurPrimalAugmentedData.mat_augmented_lagrange = &mat_augmented_lagrange;

	MatCreateShell(mpiCommGlobal,
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getRankDofs() + nRankEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges, 
		       primalDofMap.getOverallDofs() + nOverallEdges,
		       &schurPrimalAugmentedData, 
		       &mat_schur_primal);
	MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			     (void(*)(void))petscMultMatSchurPrimalAugmented);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
806

807
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
808
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
809 810
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
811 812
      KSPSetFromOptions(ksp_schur_primal);
    } else {
813 814
      MSG("Create direct schur primal solver!\n");

815 816
      TEST_EXIT(!augmentedLagrange)("Not yet supported!\n");

817 818
      double wtime = MPI::Wtime();

819 820 821
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

822 823 824
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
825

Thomas Witkowski's avatar
Thomas Witkowski committed
826
      Mat matBPi;
827 828 829
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
830 831 832
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

833

834 835 836
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
837

838 839
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

840
      for (int i = 0; i < nRowsRankB; i++) {
841
	PetscInt row = localDofMap.getStartDofs() + i;
842
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
843 844 845 846 847

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

848
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
849 850
      }

851
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
852
		primalDofMap.getLocalDofs())
853
	("Should not happen!\n");
854

855 856 857
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
858
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
859 860 861 862 863 864 865

 	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
866

Thomas Witkowski's avatar
Thomas Witkowski committed
867
	subdomain->solve(tmpVec, tmpVec);
868

Thomas Witkowski's avatar
Thomas Witkowski committed
869
	PetscScalar *tmpValues;
870
	VecGetArray(tmpVec, &tmpValues);
871
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
872
	  MatSetValue(matBPi, 
873
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
874 875 876
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
877 878 879 880 881
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
882 883
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
884

885 886
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
887 888

      Mat matPrimal;
889 890
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
891
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
892 893

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
894
      MatDestroy(&matBPi);
895

896
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
897
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
898
		      SAME_NONZERO_PATTERN);
899 900 901 902 903 904
      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);