Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

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


Thomas Witkowski's avatar
Thomas Witkowski committed
13
#include "MatrixVector.h"
14
#include "parallel/PetscHelper.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
  PetscSolverFeti::PetscSolverFeti(string name)
31
    : PetscSolver(name),
32 33 34 35 36 37
      primalDofMap(COMPONENT_WISE),
      dualDofMap(COMPONENT_WISE),
      interfaceDofMap(COMPONENT_WISE),
      localDofMap(COMPONENT_WISE),
      lagrangeMap(COMPONENT_WISE),
      interiorDofMap(COMPONENT_WISE),
38
      schurPrimalSolver(0),
39
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
40
      subdomain(NULL),
41
      massMatrixSolver(NULL),
42 43
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
44
      nGlobalOverallInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
45
      printTimings(false),
Thomas Witkowski's avatar
d  
Thomas Witkowski committed
46
      dirichletMode(0),
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
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
54 55
    Parameters::get(initFileStr + "->left precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "no") {
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 73 74 75 76
    preconditionerName = "";
    Parameters::get(initFileStr + "->right precon", preconditionerName);
    if (preconditionerName != "" && preconditionerName != "no") {
      ERROR_EXIT("FETI-DP does not support right preconditioning! (parameter \"%s->right precon\" has value \"%s\")\n",
		 initFileStr.c_str(), preconditionerName.c_str());
    }

    Parameters::get(initFileStr + "->feti->schur primal solver", schurPrimalSolver);
Thomas Witkowski's avatar
Thomas Witkowski committed
77 78 79
    TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
      ("Wrong solver \"%d\"for the Schur primal complement!\n", 
       schurPrimalSolver);
80

81 82 83 84 85 86 87 88 89 90 91
    Parameters::get(initFileStr + "->feti->stokes mode", stokesMode);
    if (stokesMode) {
      Parameters::get(initFileStr + "->feti->pressure component", pressureComponent);
      TEST_EXIT(pressureComponent >= 0)
	("FETI-DP in Stokes mode, no pressure component defined!\n");
    }
			   
    Parameters::get(initFileStr + "->feti->augmented lagrange", augmentedLagrange);

    Parameters::get(initFileStr + "->feti->symmetric", isSymmetric);

92
    Parameters::get("parallel->multi level test", multiLevelTest);
93 94
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
95 96

    Parameters::get("parallel->print timings", printTimings);
97 98 99
  }


100
  void PetscSolverFeti::initialize()
101
  {
102 103
    FUNCNAME("PetscSolverFeti::initialize()");

104 105
    MSG("INIT WITH MESH-LEVEL: %d\n", meshLevel);

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
111
    if (subdomain == NULL) {
112 113 114 115 116 117 118 119 120 121 122
      string subSolverInitStr = initFileStr + "->subsolver";
      string solverType = "petsc";
      Parameters::get(subSolverInitStr, solverType);
      OEMSolverCreator *solverCreator =
	dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType, initFileStr));
      TEST_EXIT(solverCreator)
	("No valid solver type found in parameter \"%s\"\n", 
	 initFileStr.c_str());
      solverCreator->setName(subSolverInitStr);
      subdomain = dynamic_cast<PetscSolver*>(solverCreator->create());
      subdomain->initParameters();
123
      subdomain->setSymmetric(isSymmetric);
Thomas Witkowski's avatar
d  
Thomas Witkowski committed
124 125 126 127
      if (dirichletMode == 0)
	subdomain->setHandleDirichletRows(true);
      else
	subdomain->setHandleDirichletRows(false);
128

129
      if (meshLevel == 0) {
130
	subdomain->setMeshDistributor(meshDistributor);
131
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
132
	subdomain->setMeshDistributor(meshDistributor, 
133 134
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
135
	subdomain->setLevel(meshLevel);
136
      }
137 138

      delete solverCreator;
139
    }
140

141 142 143 144
    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
145

Thomas Witkowski's avatar
Thomas Witkowski committed
146
    if (stokesMode)
147
      interfaceDofMap.init(levelData, componentSpaces, feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
148

149
    if (fetiPreconditioner == FETI_DIRICHLET) {
150 151 152
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

153
      interiorDofMap.init(levelData, componentSpaces, feSpaces, false);
154
    }
155 156 157
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
158 159 160 161
  void PetscSolverFeti::createDirichletData(Matrix<DOFMatrix*> &mat)
  {
    FUNCNAME("PetscSolverFeti::createDirichletData()");

Thomas Witkowski's avatar
d  
Thomas Witkowski committed
162 163 164 165 166 167 168 169 170
    if (dirichletMode == 1) {
      int nComponents = mat.getSize();
      for (int component = 0; component < nComponents; component++) {
	DOFMatrix* dofMat = mat[component][component];
	if (!dofMat)
	  continue;
	
	dirichletRows[component] = dofMat->getDirichletRows();
      }
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
171 172 173 174
    }
  }


175 176 177 178
  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
179 180
    double timeCounter = MPI::Wtime();

181 182
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

183 184 185 186
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
187
    if (fetiPreconditioner == FETI_DIRICHLET)
188 189
      interiorDofMap.clear();

190 191
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
192

193 194
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
195
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
196
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
197
    if (fetiPreconditioner == FETI_DIRICHLET)
198
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
199

200 201 202
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
203 204
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

Thomas Witkowski's avatar
Thomas Witkowski committed
205
    if (stokesMode) {
206 207 208 209 210
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

211 212 213 214 215 216
    int nComponents = componentSpaces.size();
    for (int component = 0; component < nComponents; component++) {
      createPrimals(component);  
      createDuals(component);
      createInterfaceNodes(component);
      createIndexB(component);
217
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
218 219 220 221

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
222
    if (fetiPreconditioner == FETI_DIRICHLET)
223
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
224

Thomas Witkowski's avatar
Thomas Witkowski committed
225
    if (stokesMode)
226 227
      interfaceDofMap.update();

228 229 230
    for (int component = 0; component < nComponents; component++) {
      createLagrange(component);
      createAugmentedLagrange(component);
231
    }
232

233 234 235
    lagrangeMap.update();


236 237 238 239 240 241 242 243 244 245 246 247
    // === ===

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

248
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
249 250 251 252 253 254 255 256 257
			   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);
    }

258
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
259 260
    for (unsigned int i = 0; i < componentSpaces.size(); i++) {
      const FiniteElemSpace *feSpace = componentSpaces[i];
261 262

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

264
      if (i == pressureComponent) {
Thomas Witkowski's avatar
Thomas Witkowski committed
265
	MSG("  nRankInterface = %d  nOverallInterface = %d\n",
266 267
	    interfaceDofMap[i].nRankDofs, 
	    interfaceDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
268
      } else {
269
	MSG("  nRankPrimals = %d   nLocalPrimals = %d  nOverallPrimals = %d\n", 
270 271 272
	    primalDofMap[i].nRankDofs, 
	    primalDofMap[i].nLocalDofs,
	    primalDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
273 274
	
	MSG("  nRankDuals = %d  nOverallDuals = %d\n",
275 276
	    dualDofMap[i].nRankDofs, 
	    dualDofMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
277 278
	
	MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
279 280
	    lagrangeMap[i].nRankDofs, 
	    lagrangeMap[i].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
281
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
282
    }
283

Thomas Witkowski's avatar
Thomas Witkowski committed
284
    subdomain->setDofMapping(&localDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
285 286 287
    subdomain->setCoarseSpaceDofMapping(&primalDofMap); 
    if (stokesMode)
      subdomain->setCoarseSpaceDofMapping(&interfaceDofMap, pressureComponent);
Thomas Witkowski's avatar
Thomas Witkowski committed
288 289

    if (printTimings) {
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
290
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
291
      timeCounter = MPI::Wtime() - timeCounter;
292
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
293 294
	  timeCounter);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
295 296 297 298 299 300


    bool writePrimals = false;
    Parameters::get("parallel->debug->write primals", writePrimals);
    if (writePrimals)
      PetscSolverFetiDebug::writePrimalFiles(*this);
301 302 303
  }


304
  void PetscSolverFeti::createPrimals(int component)
Thomas Witkowski's avatar
Thomas Witkowski committed
305
  {
306
    FUNCNAME("PetscSolverFeti::createPrimals()");  
307

308
    if (component == pressureComponent)
309 310
      return;

311 312
    const FiniteElemSpace *feSpace = componentSpaces[component];

313 314 315
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
316
    // Set of DOF indices that are considered to be primal variables.
317
    DofContainerSet& vertices = 
318
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
319 320

    DofIndexSet primals;
321
    for (DofContainerSet::iterator it = vertices.begin(); 
322
	 it != vertices.end(); ++it) {
323
      
324
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
325
	continue;
Thomas Witkowski's avatar
Thomas Witkowski committed
326

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
327 328
      if (meshLevel == 0) {
	primals.insert(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
329
      } else {
330
	double e = 1e-8;
Thomas Witkowski's avatar
Thomas Witkowski committed
331 332 333
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

334 335 336 337
	if ((fabs(c[0]) < e && fabs(c[1] - 12.5) < e) ||
	    (fabs(c[0] - 25.0) < e && fabs(c[1] - 12.5) < e) ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1]) < e) ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 25.0) < e) ||
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
338 339 340 341 342
	    (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);
	}
      }
343
    }
344

345 346 347 348

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

349
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it) {
350
      if (dofMap[feSpace].isRankDof(*it)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
351
	primalDofMap[component].insertRankDof(*it);
352
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
353
  	primalDofMap[component].insertNonRankDof(*it);
354 355
      }
    }
356 357 358
  }


359
  void PetscSolverFeti::createDuals(int component)
360 361
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
362

363
    if (component == pressureComponent)
364 365
      return;

366 367
    const FiniteElemSpace *feSpace = componentSpaces[component];

368 369 370
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
371
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
372
    
373
    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
374
	 it != allBoundaryDofs.end(); ++it) {
375
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
376 377
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
378
      if (isPrimal(component, **it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
379 380 381
	continue;

      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
382
	dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
383
      } else {
384
	if (dofMapSd[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
385
	  dualDofMap[component].insertRankDof(**it);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
386 387
      }	  
    }
388 389 390
  }

  
391
  void PetscSolverFeti::createInterfaceNodes(int component)
392 393 394
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

395
    if (component != pressureComponent)
396 397
      return;

398
    const FiniteElemSpace *feSpace = componentSpaces[component];
399 400 401 402
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
403
	 it != allBoundaryDofs.end(); ++it) {
404
      if (dirichletRows[component].count(**it))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
405 406
	continue;      
      
407
      if (dofMap[feSpace].isRankDof(**it))
Thomas Witkowski's avatar
Thomas Witkowski committed
408
	interfaceDofMap[component].insertRankDof(**it);
409
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
410
	interfaceDofMap[component].insertNonRankDof(**it);      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
411
    }
412 413 414
  }


415
  void PetscSolverFeti::createLagrange(int component)
416 417 418
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

419
    if (component == pressureComponent)
420 421
      return;

422
    const FiniteElemSpace *feSpace = componentSpaces[component];
423 424
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
425 426 427
    // 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
428
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
429

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
430
    if (meshLevel > 0) {
431
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
432 433 434 435 436 437 438 439 440

      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()) {
441
	  if (dofMapSd[feSpace].isRankDof(it.getDofIndex()))
442
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
443
	  else
444 445 446 447 448 449 450 451 452 453 454 455 456 457 458
	    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
459 460
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
461
	  if (!isPrimal(component, it.getDofIndex()))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
462 463 464
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
465

Thomas Witkowski's avatar
Thomas Witkowski committed
466
    if (dualDofMap[component].nLocalDofs == 0)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
467 468 469 470 471 472
      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)).         ===

473
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
474 475 476 477
    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
478
	if (!isPrimal(component, it.getDofIndex())) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
479 480
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

481
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
482
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
483
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
484 485 486 487
	}
      }
    }

488 489 490 491

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

492
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
493

494
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
495 496
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
497
	if (!isPrimal(component, it.getDofIndex()))
498 499 500
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
501 502 503

    stdMpi.updateSendDataSize();

504
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
505
	 !it.end(); it.nextRank()) {
506
      bool recvFromRank = false;
507
      for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
508
	if (!isPrimal(component, it.getDofIndex())) {
509 510
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
511
	       dofMapSd[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
512 513 514
	    recvFromRank = true;
	    break;
	  }
515
	}
516
      }
517 518

      if (recvFromRank)
519
	stdMpi.recv(it.getRank());
520
    }
521

522 523
    stdMpi.startCommunication();

524
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
525
	 !it.end(); it.nextRank()) {
526
      int i = 0;
527
      for (; !it.endDofIter(); it.nextDof())
Thomas Witkowski's avatar
Thomas Witkowski committed
528
	if (!isPrimal(component, it.getDofIndex()))
529 530
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
531
	       dofMapSd[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
532 533
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
534
	  else
Thomas Witkowski's avatar
Thomas Witkowski committed
535
	    lagrangeMap[component].insertNonRankDof(it.getDofIndex());
536 537
    }

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
538

539 540 541
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

542
    int nRankLagrange = 0;
543
    DofMap& dualMap = dualDofMap[component].getMap();
544
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
545
      if (dofMap[feSpace].isRankDof(it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
546
	lagrangeMap[component].insertRankDof(it->first, nRankLagrange);
547
	int degree = boundaryDofRanks[feSpace][it->first].size();
548
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
549
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
550
	lagrangeMap[component].insertNonRankDof(it->first);
551 552
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
553
    lagrangeMap[component].nRankDofs = nRankLagrange;
554 555 556
  }


557
  void PetscSolverFeti::createAugmentedLagrange(int component)
558 559 560 561 562 563 564 565
  {
    FUNCNAME("PetscSolverFeti::createAugmentedLagrange()");

    if (!augmentedLagrange)
      return;
  }


566
  void PetscSolverFeti::createIndexB(int component)
567
  {
568
    FUNCNAME("PetscSolverFeti::createIndexB()");
569

570
    const FiniteElemSpace *feSpace = componentSpaces[component];
571
    DOFAdmin* admin = feSpace->getAdmin();
572 573 574 575

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

577
    int nLocalInterior = 0;    
578
    for (int i = 0; i < admin->getUsedSize(); i++) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
579
      if (admin->isDofFree(i) ||
Thomas Witkowski's avatar
Thomas Witkowski committed
580 581 582
	  isPrimal(component, i) ||
	  isDual(component, i) ||
	  isInterface(component, i) ||
583
	  dirichletRows[component].count(i))
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
584
	continue;      
585

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
586
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
587
	localDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
588 589
	
	if (fetiPreconditioner == FETI_DIRICHLET)
590
	  interiorDofMap[component].insertRankDof(i, nLocalInterior);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
591 592 593
	
	nLocalInterior++;	
      } else {
594
	if (dofMapSd[feSpace].isRankDof(i))
Thomas Witkowski's avatar
Thomas Witkowski committed
595
	  localDofMap[component].insertRankDof(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
596
	else
Thomas Witkowski's avatar
Thomas Witkowski committed
597
	  localDofMap[component].insertNonRankDof(i);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
598 599 600
	
	TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	  ("Not yet implemnted!\n");	
601
      }
602
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
603
    
604 605
    // === And finally, add the global indicies of all dual nodes. ===

606 607
    for (DofMap::iterator it = dualDofMap[component].getMap().begin();
	 it != dualDofMap[component].getMap().end(); ++it) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
608
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
609
	localDofMap[component].insertRankDof(it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
610
      } else {
611
	if (dofMapSd[feSpace].isRankDof(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
612
	  localDofMap[component].insertRankDof(it->first);
613
	else 
Thomas Witkowski's avatar
Thomas Witkowski committed
614
	  localDofMap[component].insertNonRankDof(it->first);
615
      }
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
616
    }
617 618 619
  }


620
  void PetscSolverFeti::createMatLagrange()
621 622 623
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
624
    double wtime = MPI::Wtime();
625
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
626

627 628
    // === Create distributed matrix for Lagrange constraints. ===

629 630 631 632 633
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
634
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
635

636 637 638
    Vec vec_scale_lagrange;
    lagrangeMap.createVec(vec_scale_lagrange);

639 640 641 642 643 644 645
    // === 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.                                                     ===

646
    for (unsigned int component = 0; component < componentSpaces.size(); component++) {
647
      DofMap &dualMap = dualDofMap[component].getMap();
648

649
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
650
	TEST_EXIT_DBG(boundaryDofRanks[componentSpaces[component]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
651 652 653
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
654
	int index = lagrangeMap.getMatIndex(component, it->first);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
655

Thomas Witkowski's avatar
Thomas Witkowski committed
656
	// Copy set of all ranks that contain this dual node.
657 658
	vector<int> W(boundaryDofRanks[componentSpaces[component]][it->first].begin(), 
		      boundaryDofRanks[componentSpaces[component]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
659 660
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
661 662

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
663
	
664
	int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
665 666 667
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
668 669
	      MatSetValue(mat_lagrange, 
			  index + counter, 
670
			  localDofMap.getMatIndex(component, it->first) + rStartInterior,
671 672
			  (W[i] == mpiRank ? 1.0 : -1.0),
			  INSERT_VALUES);
673
	    }
674 675 676 677 678 679 680
	    counter++;
	  }
	}

	// === Create scaling factors for scaling the lagrange matrix, which ===
	// === is required for FETI-DP preconditioners.                      ===
	
681
	if (dofMap[componentSpaces[component]].isRankDof(it->first)) {
682 683 684 685 686 687
	  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);
688 689 690 691 692 693 694
	  }
	}
      }
    }

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

696

Thomas Witkowski's avatar
So on  
Thomas Witkowski committed
697 698 699 700 701 702 703 704 705 706 707 708
    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);


709 710 711 712 713 714 715 716 717 718 719 720 721 722 723
    // === 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
724 725
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
726 727
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB  
Thomas Witkowski committed
728
    }
729 730
  }

731

Thomas Witkowski's avatar
Thomas Witkowski committed
732 733
  bool PetscSolverFeti::testWirebasketEdge(BoundaryObject &edge, const FiniteElemSpace *feSpace)
  {
Thomas Witkowski's avatar
a  
Thomas Witkowski committed
734 735
    FUNCNAME("PetscSolverFeti::testWirebasketEdge()");

Thomas Witkowski's avatar
da  
Thomas Witkowski committed
736 737
    return true;

Thomas Witkowski's avatar
a  
Thomas Witkowski committed
738 739 740
    if (meshDistributor->getMesh()->getDim() == 2)
      return true;

Thomas Witkowski's avatar
Thomas Witkowski committed
741 742 743
    if (meshDistributor->getIntBoundary().getDegreeOwn(edge) != 3)
      return false;

Thomas Witkowski's avatar
da  
Thomas Witkowski committed
744 745
    return true;

Thomas Witkowski's avatar
Thomas Witkowski committed
746 747 748 749 750 751 752 753 754 755 756 757 758 759
    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);
  }
760

761

762
  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseEdges()
763
  {
764
    FUNCNAME("PetscSolverFeti::getAugmentedLagrange()");
765 766

    InteriorBoundary &intBound = meshDistributor->getIntBoundary();
Thomas Witkowski's avatar
Thomas Witkowski committed
767
    std::set<BoundaryObject> allEdges;
768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787
    for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it)
      if (it->rankObj.subObj == EDGE && 
	  testWirebasketEdge(it->rankObj, feSpaces[0]) &&
	  allEdges.count(it->rankObj) == 0)
	allEdges.insert(it->rankObj);

    int nEmptyEdges = 0;
    vector<vector<BoundaryObject> > data;
    for (std::set<BoundaryObject>::iterator it = allEdges.begin();
	 it != allEdges.end(); ++it) {
      DofContainer edgeDofs;
      it->el->getAllDofs(feSpaces[0], *it, edgeDofs);
      if (edgeDofs.size() == 0) {
	nEmptyEdges++;
      } else {
	vector<BoundaryObject> oneBoundary;
	oneBoundary.push_back(*it);      
	data.push_back(oneBoundary);
      }
    }
788

789 790 791
    int nEdges = allEdges.size();
    mpi::globalAdd(nEdges);
    mpi::globalAdd(nEmptyEdges);
Thomas Witkowski's avatar
Thomas Witkowski committed
792

793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808
    MSG("Coarse space edges: %d (empty: %d)\n", nEdges, nEmptyEdges);
    
    return data;
  }


  vector<vector<BoundaryObject> > PetscSolverFeti::getCoarseFaces()
  {
    FUNCNAME