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

namespace AMDiS {

  using namespace std;

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

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

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

    Parameters::get("parallel->multi level test", multiLevelTest);
75 76
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
77 78

    Parameters::get("parallel->print timings", printTimings);
79 80

    Parameters::get("parallel->feti->augmented lagrange", augmentedLagrange);
81 82

    Parameters::get("parallel->feti->symmetric", isSymmetric);
83 84 85
  }


86
  void PetscSolverFeti::initialize()
87
  {
88 89
    FUNCNAME("PetscSolverFeti::initialize()");

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
95 96 97 98 99 100
    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
101 102
    }
			   
Thomas Witkowski's avatar
Thomas Witkowski committed
103 104
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
105
      subdomain->setSymmetric(isSymmetric);
106
      subdomain->setHandleDirichletRows(false);
107

108
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
109
	subdomain->setMeshDistributor(meshDistributor, 
110
				      mpiCommGlobal, mpiCommLocal);
111
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
112
	subdomain->setMeshDistributor(meshDistributor, 
113 114
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
115
	subdomain->setLevel(meshLevel);
116 117
      }
    }
118

119 120 121 122
    primalDofMap.init(levelData, componentSpaces, feSpaces);   
    dualDofMap.init(levelData, componentSpaces, feSpaces, false);
    localDofMap.init(levelData, componentSpaces, feSpaces, meshLevel != 0);
    lagrangeMap.init(levelData, componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
123

Thomas Witkowski's avatar
Thomas Witkowski committed
124
    if (stokesMode)
125
      interfaceDofMap.init(levelData, componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
126

127
    if (fetiPreconditioner == FETI_DIRICHLET) {
128 129 130
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

131
      interiorDofMap.init(levelData, componentSpaces, feSpaces, false);
132
    }
133 134 135
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
136 137 138 139 140
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

    int nComponents = mat.getSize();
141 142
    for (int component = 0; component < nComponents; component++) {
      DOFMatrix* dofMat = mat[component][component];
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
143 144
      if (!dofMat)
	continue;
145 146

      dirichletRows[component] = dofMat->getDirichletRows();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
147 148 149 150
    }
  }


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

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

157 158
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

159 160 161 162
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
163
    if (fetiPreconditioner == FETI_DIRICHLET)
164 165
      interiorDofMap.clear();

166 167
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
168

169 170
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
171
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
172
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
173
    if (fetiPreconditioner == FETI_DIRICHLET)
174
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
175

176 177 178
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
179 180
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
181
    if (stokesMode) {
182 183 184 185 186
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

187 188 189 190 191 192
    int nComponents = componentSpaces.size();
    for (int component = 0; component < nComponents; component++) {
      createPrimals(component);  
      createDuals(component);
      createInterfaceNodes(component);
      createIndexB(component);
193
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
194 195 196 197

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
198
    if (fetiPreconditioner == FETI_DIRICHLET)
199
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
200

Thomas Witkowski's avatar
Thomas Witkowski committed
201
    if (stokesMode)
202 203
      interfaceDofMap.update();

204 205 206
    for (int component = 0; component < nComponents; component++) {
      createLagrange(component);
      createAugmentedLagrange(component);
207
    }
208

209 210 211
    lagrangeMap.update();


212 213 214 215 216 217 218 219 220 221 222 223
    // === ===

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

224
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
225 226 227 228 229 230 231 232 233
			   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);
    }

234
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
235 236
    for (unsigned int i = 0; i < componentSpaces.size(); i++) {
      const FiniteElemSpace *feSpace = componentSpaces[i];
237 238

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

240
      if (i == pressureComponent) {
Thomas Witkowski's avatar
Thomas Witkowski committed
241
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
242 243
	    interfaceDofMap[i].nRankDofs, 
	    interfaceDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
244
      } else {
245
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
246 247 248
	    primalDofMap[i].nRankDofs, 
	    primalDofMap[i].nLocalDofs,
	    primalDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
249 250
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
251 252
	    dualDofMap[i].nRankDofs, 
	    dualDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
253 254
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
255 256
	    lagrangeMap[i].nRankDofs, 
	    lagrangeMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
257
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
258
    }
259

Thomas Witkowski's avatar
Thomas Witkowski committed
260
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
261 262 263
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
264 265

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
266
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
267
      timeCounter = MPI::Wtime() - timeCounter;
268
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
269 270
	  timeCounter);
    }
271 272 273
  }


274
  void PetscSolverFeti::createPrimals(int component)
Thomas Witkowski's avatar
Thomas Witkowski committed
275
  {
276
    FUNCNAME("PetscSolverFeti::createPrimals()");  
277

278
    if (component == pressureComponent)
279 280
      return;

281 282
    const FiniteElemSpace *feSpace = componentSpaces[component];

283 284 285
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
286
    // Set of DOF indices that are considered to be primal variables.
287
    DofContainerSet& vertices = 
288
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
289 290

    DofIndexSet primals;
291
    for (DofContainerSet::iterator it = vertices.begin(); 
292
	 it != vertices.end(); ++it) {
293
      
294
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
295
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
296

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
297 298
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
299
      } else {
300
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
301 302 303
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
304 305 306 307 308 309 310
	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);
	}
      }
311
    }
312

313 314 315 316

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

317
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
318
      if (dofMap[feSpace].isRankDof(*it)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
319
	primalDofMap[component].insertRankDof(*it);
320
      } else
Thomas Witkowski's avatar
Thomas Witkowski committed
321
  	primalDofMap[component].insertNonRankDof(*it);
322 323 324
  }


325
  void PetscSolverFeti::createDuals(int component)
326 327
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
328

329
    if (component == pressureComponent)
330 331
      return;

332 333
    const FiniteElemSpace *feSpace = componentSpaces[component];

334 335 336
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
337
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
338
    
339
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
340
	 it != allBoundaryDofs.end(); ++it) {
341
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
342 343
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
344
      if (isPrimal(component, **it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
345 346 347
	continue;

      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
348
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
349
      } else {
350
	if (dofMapSd[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
351
	  dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
352 353
      }	  
    }
354 355 356
  }

  
357
  void PetscSolverFeti::createInterfaceNodes(int component)
358 359 360
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

361
    if (component != pressureComponent)
362 363
      return;

364
    const FiniteElemSpace *feSpace = componentSpaces[component];
365 366 367 368
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
369
	 it != allBoundaryDofs.end(); ++it) {
370
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
371 372
	continue;      
      
373
      if (dofMap[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
374
	interfaceDofMap[component].insertRankDof(**it);
375
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
376
	interfaceDofMap[component].insertNonRankDof(**it);      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
377
    }
378 379 380
  }


381
  void PetscSolverFeti::createLagrange(int component)
382 383 384
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

385
    if (component == pressureComponent)
386 387
      return;

388
    const FiniteElemSpace *feSpace = componentSpaces[component];
389 390
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
391 392 393
    // 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
394
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
395

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
396
    if (meshLevel > 0) {
397
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
398 399 400 401 402 403 404 405 406

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

Thomas Witkowski's avatar
Thomas Witkowski committed
432
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
433 434 435 436 437 438
      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)).         ===

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

447
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
448
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
449
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
450 451 452 453
	}
      }
    }

454 455 456 457

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

458
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
459

460
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
461 462
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
463
	if (!isPrimal(component, it.getDofIndex()))
464 465 466
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
467 468 469

    stdMpi.updateSendDataSize();

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

      if (recvFromRank)
485
	stdMpi.recv(it.getRank());
486
    }
487

488 489
    stdMpi.startCommunication();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
504

505 506 507
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


523
  void PetscSolverFeti::createAugmentedLagrange(int component)
524 525 526 527 528 529 530 531
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


532
  void PetscSolverFeti::createIndexB(int component)
533
  {
534
    FUNCNAME("PetscSolverFeti::createIndexB()");
535

536 537

    const FiniteElemSpace *feSpace = componentSpaces[component];
538
    DOFAdmin* admin = feSpace->getAdmin();
539 540 541 542

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

544
    int nLocalInterior = 0;    
545
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
546
      if (admin->isDofFree(i) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
547 548 549
	  isPrimal(component, i) ||
	  isDual(component, i) ||
	  isInterface(component, i) ||
550
	  dirichletRows[component].count(i))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
551
	continue;      
552

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
553
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
554
	localDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
555 556
	
	if (fetiPreconditioner == FETI_DIRICHLET)
557
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
558 559 560
	
	nLocalInterior++;	
      } else {
561
	if (dofMapSd[feSpace].isRankDof(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
562
	  localDofMap[component].insertRankDof(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
563
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
564
	  localDofMap[component].insertNonRankDof(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
565 566 567
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
568
      }
569
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
570
    
571 572
    // === And finally, add the global indicies of all dual nodes. ===

573 574
    for (DofMap::iterator it = dualDofMap[component].getMap().begin();
	 it != dualDofMap[component].getMap().end(); ++it) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
575
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
576
	localDofMap[component].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
577
      } else {
578
	if (dofMapSd[feSpace].isRankDof(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
579
	  localDofMap[component].insertRankDof(it->first);
580
	else 
Thomas Witkowski's avatar
Thomas Witkowski committed
581
	  localDofMap[component].insertNonRankDof(it->first);
582
      }
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
583
    }
584 585 586
  }


587
  void PetscSolverFeti::createMatLagrange()
588 589 590
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
591
    double wtime = MPI::Wtime();
592
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
593

594 595
    // === Create distributed matrix for Lagrange constraints. ===

596 597 598 599 600
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
601
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
602

603 604 605
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

606 607 608 609 610 611 612
    // === 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.                                                     ===

613
    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
614
      DofMap &dualMap = dualDofMap[component].getMap();
615

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

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

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
630
	
631
	int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
632 633 634
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
635 636
	      MatSetValue(mat_lagrange, 
			  index + counter, 
637
			  localDofMap.getMatIndex(component, it->first) + rStartInterior,
638 639
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
640
	    }
641 642 643 644 645 646 647
	    counter++;
	  }
	}

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
648
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
649 650 651 652 653 654
	  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);
655 656 657 658 659 660 661
	  }
	}
      }
    }

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

663

Thomas Witkowski's avatar
So on  
Thomas Witkowski committed
664 665 666 667 668 669 670 671 672 673 674 675
    int nZeroRows = PetscSolverFetiDebug::testZeroRows(mat_lagrange);
    int m,n;
    MatGetSize(mat_lagrange, &m ,&n);
    MSG("Lagrange matrix has %d zero rows and global size of %d %d!\n", nZeroRows, m, n);

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


676 677 678 679 680 681 682 683 684 685 686 687 688 689 690
    // === If required, create \ref mat_lagrange_scaled ===

    VecAssemblyBegin(vec_scale_lagrange);
    VecAssemblyEnd(vec_scale_lagrange);

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

    VecDestroy(&vec_scale_lagrange);


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

Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
691 692
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
693 694
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
695
    }
696 697
  }

698

Thomas Witkowski's avatar
Thomas Witkowski committed
699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714
  bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
  {
    Element *el = edge.el;
    int i0 = el->getVertexOfEdge(edge.ithObj, 0);
    int i1 = el->getVertexOfEdge(edge.ithObj, 1);
    DegreeOfFreedom d0 = el->getDof(i0, 0);
    DegreeOfFreedom d1 = el->getDof(i1, 0);
    WorldVector<double> c0, c1;
    el->getMesh()->getDofIndexCoords(d0, feSpace, c0);
    el->getMesh()->getDofIndexCoords(d1, feSpace, c1);
    bool xe = fabs(c0[0] - c1[0]) < 1e-8;
    bool ye = fabs(c0[1] - c1[1]) < 1e-8;
    bool ze = fabs(c0[2] - c1[2]) < 1e-8;
    int counter = static_cast<int>(xe) + static_cast<int>(ye) + static_cast<int>(ze);
    return (counter == 2);
  }
715

716

717
  void PetscSolverFeti::createMatAugmentedLagrange()
718 719 720 721 722 723 724 725 726 727
  {
    FUNCNAME("PetscSolverFeti::createMatAugmentedLagrange()");

    if (!augmentedLagrange)
      return;

    double wtime = MPI::Wtime();

    nOverallEdges = 0;
    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
Thomas Witkowski's avatar
Thomas Witkowski committed
728
    std::set<BoundaryObject> allEdges;
729
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
730 731 732
      if ((it->rankObj.subObj == FACE ||
	   (it->rankObj.subObj == EDGE && 
	    testWirebasketEdge(it->rankObj, feSpaces[0]))) && 
733 734 735 736 737
	  allEdges.count(it->rankObj) == 0) {
	bool dirichletOnlyEdge = true;

	DofContainer edgeDofs;
	it->rankObj.el->getAllDofs(feSpaces[0], it->rankObj, edgeDofs);
738
	
739 740
	for (DofContainer::iterator dit = edgeDofs.begin();
	     dit != edgeDofs.end(); ++dit) {
741
	  if (dirichletRows[0].count(**dit) == 0) {
742 743 744 745 746 747 748 749 750
	    dirichletOnlyEdge = false;
	    break;
	  }
	}

	if (!dirichletOnlyEdge)
	  allEdges.insert(it->rankObj);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
751 752

    nRankEdges = allEdges.size();    
753 754 755 756
    int rStartEdges = 0;
    mpi::getDofNumbering(mpiCommGlobal, nRankEdges, rStartEdges, nOverallEdges);

    MSG("nRankEdges = %d, nOverallEdges = %d\n", nRankEdges, nOverallEdges);
757
    
758 759 760
    nRankEdges *= componentSpaces.size();
    rStartEdges *= componentSpaces.size();
    nOverallEdges *= componentSpaces.size();
761 762

    MatCreateAIJ(mpiCommGlobal,
763 764
		 nRankEdges, lagrangeMap.getRankDofs(),
		 nOverallEdges, lagrangeMap.getOverallDofs(),
765
		 2, PETSC_NULL, 2, PETSC_NULL, 
766 767 768
		 &mat_augmented_lagrange);
    MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

769
    int rowCounter = rStartEdges;
Thomas Witkowski's avatar
Thomas Witkowski committed
770 771
    for (std::set<BoundaryObject>::iterator edgeIt = allEdges.begin(); 
	 edgeIt != allEdges.end(); ++edgeIt) {
772
      for (int component = 0; component < componentSpaces.size(); component++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
773
	DofContainer edgeDofs;
774
	edgeIt->el->getAllDofs(componentSpaces[component], *edgeIt, edgeDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
775 776 777
	
	for (DofContainer::iterator it = edgeDofs.begin();
	     it != edgeDofs.end(); ++it) {
Thomas Witkowski's avatar
Thomas Witkowski committed
778
	  TEST_EXIT_DBG(isPrimal(component, **it) == false)
Thomas Witkowski's avatar
Thomas Witkowski committed
779 780
	    ("Should not be primal!\n");
	  
781
	  if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
782 783
	    continue;