PetscSolverSchur.cc 15 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
20
21
22


#include "parallel/PetscSolverSchur.h"
23
#include "parallel/PetscHelper.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24
25
26
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

27
namespace AMDiS { namespace Parallel {
Thomas Witkowski's avatar
Thomas Witkowski committed
28
29
30
31
32
33
34
35

  using namespace std;

  void PetscSolverSchur::updateDofData(int nComponents)
  {
    FUNCNAME("PetscSolverSchur::updateDofData()");

    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
36
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
37

38
    const FiniteElemSpace *feSpace = componentSpaces[0];
39
    Mesh* mesh = feSpace->getMesh();
40
41
//     typedef map<int, DofContainer> RankToDofContainer;
#if 0
Thomas Witkowski's avatar
Thomas Witkowski committed
42
    typedef map<DegreeOfFreedom, bool> DofIndexToBool;
43
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
44
45
46

    boundaryDofs.clear();
    std::set<DegreeOfFreedom> boundaryLocalDofs;
47
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, 0).getSendDofs(), feSpace);
48
49
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
50
	boundaryLocalDofs.insert(it.getDofIndex());
51
	boundaryDofs.insert((*interiorMap)[feSpace][it.getDofIndex()].global);
Thomas Witkowski's avatar
Thomas Witkowski committed
52
      }
53

54

Thomas Witkowski's avatar
Thomas Witkowski committed
55
    nBoundaryDofs = boundaryDofs.size();
56
    mpi::getDofNumbering(domainComm, nBoundaryDofs,
Thomas Witkowski's avatar
Thomas Witkowski committed
57
58
			 rStartBoundaryDofs, nOverallBoundaryDofs);

59

60
    DofContainerSet& edgeDofs =
61
      meshDistributor->getBoundaryDofInfo(feSpace, 0).geoDofs[EDGE];
62
    DofContainerSet& vertexDofs =
63
      meshDistributor->getBoundaryDofInfo(feSpace, 0).geoDofs[VERTEX];
64
65
66
67
68
69
70
    int nEdgeDofs = edgeDofs.size();
    int nVertexDofs = vertexDofs.size();

    TEST_EXIT_DBG(nEdgeDofs + nVertexDofs == nBoundaryDofs)
      ("Should not happen!\n");

    int rStartEdgeDofs, nOverallEdgeDofs;
71
    mpi::getDofNumbering(domainComm, nEdgeDofs,
72
73
74
			 rStartEdgeDofs, nOverallEdgeDofs);

    int rStartVertexDofs, nOverallVertexDofs;
75
    mpi::getDofNumbering(domainComm, nVertexDofs,
76
77
78
79
			 rStartVertexDofs, nOverallVertexDofs);

    TEST_EXIT_DBG(nOverallEdgeDofs + nOverallVertexDofs == nOverallBoundaryDofs)
      ("Should not happen!\n");
80

81

Thomas Witkowski's avatar
Thomas Witkowski committed
82
    mapGlobalBoundaryDof.clear();
83
84
85
86
87
#if 1
    {
      int counter = rStartEdgeDofs;
      for (DofContainerSet::iterator it = edgeDofs.begin();
	   it != edgeDofs.end(); ++it)
88
	mapGlobalBoundaryDof[(*interiorMap)[feSpace][**it].global] =
89
90
91
92
93
94
	  counter++;
    }
    {
      int counter = nOverallEdgeDofs + rStartVertexDofs;
      for (DofContainerSet::iterator it = vertexDofs.begin();
	   it != vertexDofs.end(); ++it)
95
	mapGlobalBoundaryDof[(*interiorMap)[feSpace][**it].global] =
96
97
98
99
100
	  counter++;
    }
#else
 {
    int counter = rStartBoundaryDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
101
102
103
    for (std::set<DegreeOfFreedom>::iterator it = boundaryDofs.begin();
	 it != boundaryDofs.end(); ++it)
      mapGlobalBoundaryDof[*it] = counter++;
104
105
 }
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
106
107
108


    std::set<DegreeOfFreedom> otherBoundaryLocalDofs;
109
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, 0).getRecvDofs(), feSpace);
110
111
112
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
	otherBoundaryLocalDofs.insert(it.getDofIndex());
113

Thomas Witkowski's avatar
Thomas Witkowski committed
114
    interiorDofs.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
115
116
117

    ERROR_EXIT("Rewrite the following code block!\n");
#if 0
118
    DofIndexToBool& isRankDof = meshDistributor->getIsRankDof(feSpace);
119
    for (DofIndexToBool::iterator dofIt = isRankDof.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
120
	 dofIt != isRankDof.end(); ++dofIt) {
121
122
      if (dofIt->second &&
	  boundaryLocalDofs.count(dofIt->first) == 0 &&
Thomas Witkowski's avatar
Thomas Witkowski committed
123
	  otherBoundaryLocalDofs.count(dofIt->first) == 0)
124
	interiorDofs.insert(meshDistributor->mapDofToGlobal(feSpace, dofIt->first));
Thomas Witkowski's avatar
Thomas Witkowski committed
125
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
126
#endif
127

Thomas Witkowski's avatar
Thomas Witkowski committed
128
    nInteriorDofs = interiorDofs.size();
129
    mpi::getDofNumbering(domainComm, nInteriorDofs,
Thomas Witkowski's avatar
Thomas Witkowski committed
130
131
			 rStartInteriorDofs, nOverallInteriorDofs);

132
133
134
135
136
137
138
    {
      int counter = rStartInteriorDofs;
      mapGlobalInteriorDof.clear();
      for (std::set<DegreeOfFreedom>::iterator it = interiorDofs.begin();
	   it != interiorDofs.end(); ++it)
	mapGlobalInteriorDof[*it] = counter++;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
139
140
141
142
143


    TEST_EXIT_DBG(nInteriorDofs > 0)("Should not happen!\n");


144
    StdMpi<vector<DegreeOfFreedom> > stdMpi(domainComm);
145
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, 0).getSendDofs(), feSpace);
146
147
148
149
150
	 !it.end(); it.nextRank()) {
      stdMpi.getSendData(it.getRank()).resize(0);
      stdMpi.getSendData(it.getRank()).reserve(it.getDofs().size());

      for (; !it.endDofIter(); it.nextDof()) {
151
	int globalSendDof = (*interiorMap)[feSpace][it.getDofIndex()].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
152
153
154
155

	TEST_EXIT_DBG(mapGlobalBoundaryDof.count(globalSendDof))
	  ("No mapping for boundary DOF %d!\n", globalSendDof);

156
	stdMpi.getSendData(it.getRank()).push_back(mapGlobalBoundaryDof[globalSendDof]);
Thomas Witkowski's avatar
Thomas Witkowski committed
157
      }
158
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
159
160

    stdMpi.updateSendDataSize();
161

162
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, 0).getRecvDofs(), feSpace);
163
164
	 !it.end(); it.nextRank())
      stdMpi.recv(it.getRank());
Thomas Witkowski's avatar
Thomas Witkowski committed
165

166
    stdMpi.startCommunication();
167

168
    for (DofComm::Iterator it(meshDistributor->getDofComm(mesh, 0).getRecvDofs(), feSpace);
169
170
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
171
	int globalRecvDof = (*interiorMap)[feSpace][it.getDofIndex()].global;
172
	mapGlobalBoundaryDof[globalRecvDof] =
173
	  stdMpi.getRecvData(it.getRank())[it.getDofCounter()];
Thomas Witkowski's avatar
Thomas Witkowski committed
174
175
176
177
178
179
	boundaryDofs.insert(globalRecvDof);
      }


    // === Create PETSc IS structurs for interior and boundary DOFs. ===

180
    ISCreateStride(domainComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
181
		   nInteriorDofs * nComponents,
182
		   (rStartInteriorDofs + rStartBoundaryDofs) * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
		   1, &interiorIs);

185
    ISCreateStride(domainComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
186
		   nBoundaryDofs * nComponents,
187
		   (rStartInteriorDofs + rStartBoundaryDofs + nInteriorDofs) * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
188
189
190
191
		   1, &boundaryIs);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
192
  void PetscSolverSchur::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
193
  {
194
//     FUNCNAME("PetscSolverSchur::fillPetscMatrix()");
195

196
    createMatVec(*seqMat);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
197

198
    const FiniteElemSpace *feSpace = componentSpaces[0];
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
199
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
200
201
202
203
204
205
206
207
    updateDofData(nComponents);

    int nInteriorRows = nInteriorDofs * nComponents;
    int nOverallInteriorRows = nOverallInteriorDofs * nComponents;
    int nBoundaryRows = nBoundaryDofs * nComponents;
    int nOverallBoundaryRows = nOverallBoundaryDofs * nComponents;


208
209
    MatCreateAIJ(domainComm,
		 nInteriorRows, nInteriorRows,
210
211
212
		 nOverallInteriorRows, nOverallInteriorRows,
		 100, PETSC_NULL, 100, PETSC_NULL, &matA11);

213
214
    MatCreateAIJ(domainComm,
		 nBoundaryRows, nBoundaryRows,
215
216
217
		 nOverallBoundaryRows, nOverallBoundaryRows,
		 100, PETSC_NULL, 100, PETSC_NULL, &matA22);

218
219
    MatCreateAIJ(domainComm,
		 nInteriorRows, nBoundaryRows,
220
		 nOverallInteriorRows, nOverallBoundaryRows,
221
		 100, PETSC_NULL, 100, PETSC_NULL, &matA12);
222

223
224
    MatCreateAIJ(domainComm,
		 nBoundaryRows, nInteriorRows,
225
226
		 nOverallBoundaryRows, nOverallInteriorRows,
		 100, PETSC_NULL, 100, PETSC_NULL, &matA21);
Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
229
230


    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
231
232
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], nComponents, i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255

    MatAssemblyBegin(matA11, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(matA11, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(matA12, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(matA12, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(matA21, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(matA21, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(matA22, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(matA22, MAT_FINAL_ASSEMBLY);

    Mat tmpMat[2][2];
    tmpMat[0][0] = matA11;
    tmpMat[0][1] = matA12;
    tmpMat[1][0] = matA21;
    tmpMat[1][1] = matA22;

    IS tmpIS[2];
    tmpIS[0] = interiorIs;
    tmpIS[1] = boundaryIs;

256
    MatCreateNest(domainComm, 2, &tmpIS[0], 2, &tmpIS[0], &tmpMat[0][0],
257
258
		  &getMatInterior());
    MatNestSetVecType(getMatInterior(), VECNEST);
Thomas Witkowski's avatar
Thomas Witkowski committed
259

260
    matAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
261

262
263
    int nRankRows = (*interiorMap)[feSpace].nRankDofs * nComponents;
    int nOverallRows = (*interiorMap)[feSpace].nOverallDofs * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
264

265
    VecCreateMPI(domainComm, nRankRows, nOverallRows, &petscSolVec);
266
267
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
268

269
270
  void PetscSolverSchur::fillPetscRhs(SystemVector *vec)
  {
271
//     FUNCNAME("PetscSolverSchur::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
272

273
    const FiniteElemSpace *feSpace = componentSpaces[0];
274
    int nComponents = vec->getSize();
275
276
    int nRankRows = (*interiorMap)[feSpace].nRankDofs * nComponents;
    int nOverallRows = (*interiorMap)[feSpace].nOverallDofs * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
277

278
    VecCreateMPI(domainComm, nRankRows, nOverallRows, &(getVecRhsInterior()));
Thomas Witkowski's avatar
Thomas Witkowski committed
279
280

    for (int i = 0; i < nComponents; i++)
281
      setDofVector(getVecRhsInterior(), vec->getDOFVector(i), nComponents, i);
Thomas Witkowski's avatar
Thomas Witkowski committed
282

283
    vecRhsAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
284
285
286
  }


287
  void PetscSolverSchur::solvePetscMatrix(SystemVector &vec,
288
					  AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
289
  {
290
//     FUNCNAME("PetscSolverSchur::solvePetscMatrix()");
Thomas Witkowski's avatar
Thomas Witkowski committed
291

292
    const FiniteElemSpace *feSpace = componentSpaces[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
    int nComponents = vec.getSize();

295
    KSPCreate(domainComm, &kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
296

297
    petsc::ksp_set_operators(kspInterior, getMatInterior(), getMatInterior());
298
299
    KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetFromOptions(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
300

301
302
303
304
305
    KSPGetPC(kspInterior, &pcInterior);
    PCSetType(pcInterior, PCFIELDSPLIT);
    PCFieldSplitSetIS(pcInterior, "interior", interiorIs);
    PCFieldSplitSetIS(pcInterior, "boundary", boundaryIs);
    PCSetFromOptions(pcInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
306
307


308
    KSPSolve(kspInterior, getVecRhsInterior(), petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
309
310
311
312
313
314
315
316
317

    // === Transfere values from PETSc's solution vectors to AMDiS vectors. ===

    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>::Iterator dofIt(vec.getDOFVector(i), USED_DOFS);
      for (dofIt.reset(); !dofIt.end(); ++dofIt) {
318
	DegreeOfFreedom globalRowDof = (*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
319
320
321
322
323
	if (boundaryDofs.count(globalRowDof)) {
	  int index =
	    (mapGlobalBoundaryDof[globalRowDof] - rStartBoundaryDofs + nInteriorDofs) * (i + 1);
	  *dofIt = vecPointer[index];
	} else {
324
	  int index =
Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
	    (mapGlobalInteriorDof[globalRowDof] - rStartInteriorDofs) * (i + 1);
	  *dofIt = vecPointer[index];
	}
      }
    }

    VecRestoreArray(petscSolVec, &vecPointer);


    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===
    // === more than one partition.                                       ===
    meshDistributor->synchVector(vec);


    // === Destroy PETSC's variables. ===

341
342
343
344
    MatDestroy(&matA11);
    MatDestroy(&matA12);
    MatDestroy(&matA21);
    MatDestroy(&matA22);
345

346
347
    matDestroy();
    vecDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
348

349
    KSPDestroy(&kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
350
351
  }

352
353

  void PetscSolverSchur::setDofMatrix(DOFMatrix* seqMat, int dispMult,
Thomas Witkowski's avatar
Thomas Witkowski committed
354
355
356
357
				      int dispAddRow, int dispAddCol)
  {
    FUNCNAME("PetscSolverSchur::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
358
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
359

360
    const FiniteElemSpace* feSpace = componentSpaces[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
363
364
    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits= mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
365
366
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
367
368
369
370
371
372
373
374
375
376
377

    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

    vector<int> colsBoundary, colsInterior;
    vector<double> valuesBoundary, valuesInterior;
    colsBoundary.reserve(300);
    colsInterior.reserve(300);
    valuesBoundary.reserve(300);
    valuesInterior.reserve(300);

378
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()),
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
379
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
380
381

      // Global index of the current row DOF.
382
      int globalRowDof = (*interiorMap)[feSpace][cursor.value()].global;
383

Thomas Witkowski's avatar
Thomas Witkowski committed
384
385
386
387
      colsBoundary.clear();
      colsInterior.clear();
      valuesBoundary.clear();
      valuesInterior.clear();
388
389

      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
390
	   icursor != icend; ++icursor) {
391
	int globalColDof = (*interiorMap)[feSpace][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
394
395
396

	if (boundaryDofs.count(globalColDof)) {
	  TEST_EXIT_DBG(mapGlobalBoundaryDof.count(globalColDof))
	    ("Should not happen!\n");

397
	  int colIndex =
Thomas Witkowski's avatar
Thomas Witkowski committed
398
	    mapGlobalBoundaryDof[globalColDof] * dispMult + dispAddCol;
399

Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
402
403
	  colsBoundary.push_back(colIndex);
	  valuesBoundary.push_back(value(*icursor));
	} else {
	  TEST_EXIT_DBG(mapGlobalInteriorDof.count(globalColDof))
404
	    ("Cannot find global interior mapping for global column DOF %d!\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
	     globalColDof);

407
	  int colIndex =
Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
410
411
412
413
414
415
416
417
418
	    mapGlobalInteriorDof[globalColDof] * dispMult + dispAddCol;

	  colsInterior.push_back(colIndex);
	  valuesInterior.push_back(value(*icursor));
	}
      }

      if (boundaryDofs.count(globalRowDof)) {
	TEST_EXIT_DBG(mapGlobalBoundaryDof.count(globalRowDof))
	  ("Should not happen!\n");

419
	int rowIndex =
Thomas Witkowski's avatar
Thomas Witkowski committed
420
421
	  mapGlobalBoundaryDof[globalRowDof] * dispMult + dispAddRow;

422
 	MatSetValues(matA22, 1, &rowIndex, colsBoundary.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
423
 		     &(colsBoundary[0]), &(valuesBoundary[0]), ADD_VALUES);
424
 	MatSetValues(matA21, 1, &rowIndex, colsInterior.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
425
426
427
 		     &(colsInterior[0]), &(valuesInterior[0]), ADD_VALUES);
      } else {
	TEST_EXIT_DBG(mapGlobalInteriorDof.count(globalRowDof))
428
	  ("Cannot find global interior mapping for global row DOF %d!\n",
Thomas Witkowski's avatar
Thomas Witkowski committed
429
430
	   globalRowDof);

431
	int rowIndex =
Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
	  mapGlobalInteriorDof[globalRowDof] * dispMult + dispAddRow;

434
  	MatSetValues(matA11, 1, &rowIndex, colsInterior.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
435
  		     &(colsInterior[0]), &(valuesInterior[0]), ADD_VALUES);
436
  	MatSetValues(matA12, 1, &rowIndex, colsBoundary.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
437
438
439
440
441
442
443
444
445
  		     &(colsBoundary[0]), &(valuesBoundary[0]), ADD_VALUES);
      }
    }
  }


  void PetscSolverSchur::setDofVector(Vec& petscVec, DOFVector<double>* vec,
				      int dispMult, int dispAdd, bool rankOnly)
  {
446
    FUNCNAME_DBG("PetscSolverSchur::setDofVector()");
447

448
    const FiniteElemSpace *feSpace = componentSpaces[0];
449

Thomas Witkowski's avatar
Thomas Witkowski committed
450
451
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
452
      if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
453
454
	continue;

455
      // Calculate global row index of the DOF.
456
      DegreeOfFreedom globalRowDof =
457
	(*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
460
461
462
463
      double value = *dofIt;

      if (boundaryDofs.count(globalRowDof)) {
	TEST_EXIT_DBG(mapGlobalBoundaryDof.count(globalRowDof))
	  ("Should not happen!\n");

464
465
	int index =
	  (rStartInteriorDofs +
Thomas Witkowski's avatar
Thomas Witkowski committed
466
467
	   nInteriorDofs +
	   mapGlobalBoundaryDof[globalRowDof]) * dispMult + dispAdd;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
468
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
469
470
471
      } else {
	TEST_EXIT_DBG(mapGlobalInteriorDof.count(globalRowDof))
	  ("Should not happen!\n");
472
473
474

	int index =
	  (rStartBoundaryDofs +
Thomas Witkowski's avatar
Thomas Witkowski committed
475
	   mapGlobalInteriorDof[globalRowDof]) * dispMult + dispAdd;
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
476
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
477
478
479
480
      }
    }
  }

481
} } // end namespace Parallel, AMDiS