PetscSolverFeti.cc 68.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
#include "parallel/PetscSolverFetiStructs.h"
17 18
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
19 20
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
21
#include "parallel/PetscSolverGlobalMatrix.h"
22
#include "io/VtkWriter.h"
23 24 25 26 27

namespace AMDiS {

  using namespace std;

28 29
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
30
      schurPrimalSolver(0),
31
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
32
      subdomain(NULL),
33 34
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
35
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
36
      printTimings(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
37 38
      stokesMode(false),
      pressureComponent(-1)
39 40 41 42 43 44
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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

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

    Parameters::get("parallel->print timings", printTimings);
68 69 70
  }


71
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
72
  {
73 74
    FUNCNAME("PetscSolverFeti::initialize()");

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

78
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
79
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
80

Thomas Witkowski's avatar
Thomas Witkowski committed
81

Thomas Witkowski's avatar
Thomas Witkowski committed
82 83 84 85 86 87 88
    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");

Thomas Witkowski's avatar
Thomas Witkowski committed
89 90
      MSG("========= !!!! SUPER WARNING !!! ========\n");
      MSG("Use make use of stokes mode which is work in progress! We guarantee for nothing and even less than this!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
91
      pressureFeSpace = feSpaces[pressureComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
92 93
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
94 95
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
96

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

108 109 110 111
    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
112

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

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

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


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

Thomas Witkowski's avatar
Thomas Witkowski committed
129 130
    double timeCounter = MPI::Wtime();

131 132 133
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
134

135 136
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

137 138 139 140
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
141
    if (fetiPreconditioner == FETI_DIRICHLET)
142 143
      interiorDofMap.clear();

144 145
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
146

147 148
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
149
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
150
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
151
    if (fetiPreconditioner == FETI_DIRICHLET)
152
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
153

154 155 156
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
157 158
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
159
    if (stokesMode) {
160 161 162 163 164
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

165 166
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
167
    
168 169
      createPrimals(feSpace);  

170 171 172
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
173

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
174
      createIndexB(feSpace);     
175
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
176 177 178 179

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
180
    if (fetiPreconditioner == FETI_DIRICHLET)
181
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
182

Thomas Witkowski's avatar
Thomas Witkowski committed
183
    if (stokesMode)
184 185
      interfaceDofMap.update();

186 187 188 189
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
190

191 192 193
    lagrangeMap.update();


194 195 196 197 198 199 200 201 202 203 204 205
    // === ===

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

206
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
207 208 209 210 211 212 213 214 215
			   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);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
221
      if (feSpace == pressureFeSpace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241
	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);

	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
242
    }
243

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
245 246 247
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
248 249

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
250
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
251
      timeCounter = MPI::Wtime() - timeCounter;
252
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
253 254
	  timeCounter);
    }
255 256 257
  }


258
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
259
  {
260
    FUNCNAME("PetscSolverFeti::createPrimals()");  
261

Thomas Witkowski's avatar
Thomas Witkowski committed
262
    if (feSpace == pressureFeSpace)
263 264
      return;

265 266 267
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    // Set of DOF indices that are considered to be primal variables.
269
    DofContainerSet& vertices = 
270
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
271 272

    DofIndexSet primals;
273
    for (DofContainerSet::iterator it = vertices.begin(); 
274
	 it != vertices.end(); ++it) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
275
      if (meshLevel == 0) {
276
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
277
      } else {
278
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
279 280 281
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
282 283 284 285 286 287 288
	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);
	}
      }
289
    }
290

291 292 293 294

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

295
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
296 297
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
298
      else
299
  	primalDofMap[feSpace].insertNonRankDof(*it);
300 301 302
  }


303
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
304 305
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
306

Thomas Witkowski's avatar
Thomas Witkowski committed
307
    if (feSpace == pressureFeSpace)
308 309
      return;

310 311 312
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
313
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
314 315 316 317

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
318 319
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
320 321 322 323
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
324 325 326
  }

  
327 328 329 330
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

Thomas Witkowski's avatar
Thomas Witkowski committed
331
    if (feSpace != pressureFeSpace)
332 333 334 335 336 337 338
      return;

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
339
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
340
	interfaceDofMap[feSpace].insertRankDof(**it);
341
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
342
	interfaceDofMap[feSpace].insertNonRankDof(**it);
343 344 345
  }


346 347 348 349
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
350
    if (feSpace == pressureFeSpace)
351 352
      return;

353 354
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
355 356 357
    // 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
358
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
359

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
360
    if (meshLevel > 0) {
361
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
362 363 364 365 366 367 368 369 370

      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
371
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
372
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
373
	  else
374 375 376 377 378 379 380 381 382 383 384 385 386 387 388
	    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
389 390 391 392 393 394
	   !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());
    }
395

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
396 397 398 399 400 401 402
    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)).         ===

403
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
404 405 406 407 408 409 410
    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);

411
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
412
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
413
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
414 415 416 417
	}
      }
    }

418 419 420 421

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

422
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
423

424
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
425 426
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
427
	if (!isPrimal(feSpace, it.getDofIndex()))
428 429 430
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
431 432 433

    stdMpi.updateSendDataSize();

434
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
435
	 !it.end(); it.nextRank()) {
436
      bool recvFromRank = false;
437
      for (; !it.endDofIter(); it.nextDof()) {
438
	if (!isPrimal(feSpace, it.getDofIndex())) {
439 440 441
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
442 443 444
	    recvFromRank = true;
	    break;
	  }
445
	}
446
      }
447 448

      if (recvFromRank)
449
	stdMpi.recv(it.getRank());
450
    }
451

452 453
    stdMpi.startCommunication();

454
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
455
	 !it.end(); it.nextRank()) {
456
      int i = 0;
457
      for (; !it.endDofIter(); it.nextDof())
458
	if (!isPrimal(feSpace, it.getDofIndex()))
459 460 461
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
462 463
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
464 465
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
466 467
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
468

469 470 471
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

472
    int nRankLagrange = 0;
473
    DofMap& dualMap = dualDofMap[feSpace].getMap();
474
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
475 476
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
477
	int degree = boundaryDofRanks[feSpace][it->first].size();
478
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
479
      } else {
480
	lagrangeMap[feSpace].insertNonRankDof(it->first);
481 482
      }
    }
483
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
484 485 486
  }


487
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
488
  {
489
    FUNCNAME("PetscSolverFeti::createIndexB()");
490

491
    DOFAdmin* admin = feSpace->getAdmin();
492 493 494 495

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

497
    int nLocalInterior = 0;    
498
    for (int i = 0; i < admin->getUsedSize(); i++) {
499
      if (admin->isDofFree(i) == false && 
500
	  isPrimal(feSpace, i) == false &&
501 502
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
503 504
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
505

506
	  if (fetiPreconditioner == FETI_DIRICHLET)
507 508 509 510 511 512 513 514
	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	  
	  nLocalInterior++;	  
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	    localDofMap[feSpace].insertRankDof(i);
	  else
	    localDofMap[feSpace].insertNonRankDof(i);
515 516 517

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
518
	}
519
      }
520
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
521
    
522 523
    // === And finally, add the global indicies of all dual nodes. ===

524 525 526 527 528 529 530 531 532 533
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
      if (meshLevel == 0)
	localDofMap[feSpace].insertRankDof(it->first);
      else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
534 535 536
  }


537
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
538 539 540
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
541
    double wtime = MPI::Wtime();
542
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
543

544 545
    // === Create distributed matrix for Lagrange constraints. ===

546 547 548 549 550
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
551
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
552

553 554 555 556 557 558 559
    // === 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
560
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
561
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
562

563
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
564
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
565 566 567
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
568
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
569

Thomas Witkowski's avatar
Thomas Witkowski committed
570
	// Copy set of all ranks that contain this dual node.
571 572
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
573 574
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
575 576

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
577 578 579 580
	
	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
581 582
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
583

584
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
585
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
586
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
587
	    index++;	      
588 589 590 591 592 593 594
	  }
	}
      }
    }

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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
596 597
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
598 599
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
600
    }
601 602 603
  }


604
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
605
  {
606
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
607

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

Thomas Witkowski's avatar
Thomas Witkowski committed
611
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
612
      
613 614
      localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
      primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
615

616
      MatCreateShell(mpiCommGlobal,
617 618 619 620
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
621 622 623 624 625
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
626
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
627
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
628 629
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
630 631
      KSPSetFromOptions(ksp_schur_primal);
    } else {
632 633
      MSG("Create direct schur primal solver!\n");

634 635
      double wtime = MPI::Wtime();

636 637 638
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

639 640 641
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
642

Thomas Witkowski's avatar
Thomas Witkowski committed
643
      Mat matBPi;
644 645 646
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
647 648 649
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

650

651 652 653
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
654

655 656
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

657
      for (int i = 0; i < nRowsRankB; i++) {
658
	PetscInt row = localDofMap.getStartDofs() + i;
659
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
660 661 662 663 664

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

665
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
666 667
      }

668
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
669
		primalDofMap.getLocalDofs())
670
	("Should not happen!\n");
671

672 673 674
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
675
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
676 677 678 679 680 681 682

 	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);
683
	
Thomas Witkowski's avatar
Thomas Witkowski committed
684
	subdomain->solve(tmpVec, tmpVec);
685

Thomas Witkowski's avatar
Thomas Witkowski committed
686
	PetscScalar *tmpValues;
687
	VecGetArray(tmpVec, &tmpValues);
688
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
689
	  MatSetValue(matBPi, 
690
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
691 692 693
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
694 695 696 697 698
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
699 700
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
701

702 703
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
704 705

      Mat matPrimal;
706 707
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
708
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
709 710

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
711
      MatDestroy(&matBPi);
712

713
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
714
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
715
		      SAME_NONZERO_PATTERN);
716 717 718 719 720 721
      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);
722
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
723

Thomas Witkowski's avatar
Thomas Witkowski committed
724
      if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
725
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
726 727 728 729 730 731 732 733 734 735
	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
736
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
737 738 739
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
740
    }
741 742 743 744 745 746 747
  }


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

748
    if (schurPrimalSolver == 0) {
749
      schurPrimalData.subSolver = NULL;
750

751 752 753
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
    }
754 755 756
    
    MatDestroy(&mat_schur_primal);
    KSPDestroy(&ksp_schur_primal);
757 758 759
  }


760
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
761 762 763
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

764 765
    // === Create FETI-DP solver object. ===

766 767 768 769
    fetiData.mat_lagrange = &mat_lagrange;
    fetiData.subSolver = subdomain;
    fetiData.ksp_schur_primal = &ksp_schur_primal;

Thomas Witkowski's avatar
Thomas Witkowski committed
770
    localDofMap.createVec(fetiData.tmp_vec_b0, nGlobalOverallInterior);
771
    lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
772
    primalDofMap.createVec(fetiData.tmp_vec_primal0);
773

Thomas Witkowski's avatar
Thomas Witkowski committed
774
    if (stokesMode == false) {      
775 776 777 778 779 780 781 782
      MatCreateShell(mpiCommGlobal,
		     lagrangeMap.getRankDofs(), 
		     lagrangeMap.getRankDofs(),
		     lagrangeMap.getOverallDofs(), 
		     lagrangeMap.getOverallDofs(),
		     &fetiData, &mat_feti);
      MatShellSetOperation(mat_feti, MATOP_MULT, 
			   (void(*)(void))petscMultMatFeti);
783
    }  else {
Thomas Witkowski's avatar
Thomas Witkowski committed
784 785 786 787
      localDofMap.createVec(fetiData.tmp_vec_b1, nGlobalOverallInterior);
      primalDofMap.createVec(fetiData.tmp_vec_primal1);
      interfaceDofMap.createVec(fetiData.tmp_vec_interface);

788 789 790 791 792
      MatCreateShell(mpiCommGlobal,
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(), 
		     interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(), 
		     interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
793
		     &fetiData, &mat_feti);
794
      MatShellSetOperation(mat_feti, MATOP_MULT, 
Thomas Witkowski's avatar
Thomas Witkowski committed
795
			   (void(*)(void))petscMultMatFetiInterface);      
796
    }
797

798