PetscSolverFeti.cc 69.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
      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 257 258 259 260 261 262 263 264 265 266 267 268
	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
269 270 271
// 	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
272
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
273
    }
274

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

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


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

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

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

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

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

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

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

325 326 327 328

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

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


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

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

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

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

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

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

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

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


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

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

397 398
    boundaryDofRanks[feSpace].clear();

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

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

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

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

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

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

462 463 464 465

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

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

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

    stdMpi.updateSendDataSize();

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

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

496 497
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
512

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

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


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

    if (!augmentedLagrange)
      return;
  }


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

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

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

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

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

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


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

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

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

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

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

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
633 634 635 636
	
	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
637 638
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
639

640
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
641
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
642
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
643
	    index++;	      
644 645 646 647 648 649 650
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
652 653
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
654 655
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
656
    }
657 658 659
  }


660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685
  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,
686 687
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
688 689 690 691
		 1, PETSC_NULL, 1, PETSC_NULL, 
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

692
    int rowCounter = rStartEdges;
693 694 695 696 697 698 699 700 701 702 703 704 705
    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");

706
	    int col = lagrangeMap.getMatIndex(i, **it);
707
	    double value = 1.0;
708
	    MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
709 710
	  }

711
	  rowCounter++;
712 713 714 715 716 717 718 719 720 721 722 723 724 725 726
	}
      }
    }

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


727
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
728
  {
729
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
730

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

734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769
      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
770

771
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
772
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
773 774
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
775 776
      KSPSetFromOptions(ksp_schur_primal);
    } else {
777 778
      MSG("Create direct schur primal solver!\n");

779 780
      TEST_EXIT(!augmentedLagrange)("Not yet supported!\n");

781 782
      double wtime = MPI::Wtime();

783 784 785
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

786 787 788
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
789

Thomas Witkowski's avatar
Thomas Witkowski committed
790
      Mat matBPi;
791 792 793
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
794 795 796
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

797

798 799 800
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
801

802 803
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

804
      for (int i = 0; i < nRowsRankB; i++) {
805
	PetscInt row = localDofMap.getStartDofs() + i;
806
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
807 808 809 810 811

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

812
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
813 814
      }

815
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
816
		primalDofMap.getLocalDofs())
817
	("Should not happen!\n");
818

819 820 821
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
822
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
823 824 825 826 827 828 829

 	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
830

Thomas Witkowski's avatar
Thomas Witkowski committed
831
	subdomain->solve(tmpVec, tmpVec);
832

Thomas Witkowski's avatar
Thomas Witkowski committed
833
	PetscScalar *tmpValues;
834
	VecGetArray(tmpVec, &tmpValues);
835
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
836
	  MatSetValue(matBPi, 
837
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
838 839 840
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
841 842 843 844 845
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
846 847
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
848

849 850
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
851 852

      Mat matPrimal;
853 854
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
855
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
856 857

      MatDestroy(&matPrimal);