PetscSolverSchur.cc 14.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * 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.
 * 
 ******************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
20
21
22
23
24
25


#include "parallel/PetscSolverSchur.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

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

  using namespace std;

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

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

37
    const FiniteElemSpace *feSpace = componentSpaces[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
38
39
40
41
42
    typedef map<int, DofContainer> RankToDofContainer;
    typedef map<DegreeOfFreedom, bool> DofIndexToBool;

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

50
      
Thomas Witkowski's avatar
Thomas Witkowski committed
51
    nBoundaryDofs = boundaryDofs.size();
52
    mpi::getDofNumbering(domainComm, nBoundaryDofs, 
Thomas Witkowski's avatar
Thomas Witkowski committed
53
54
			 rStartBoundaryDofs, nOverallBoundaryDofs);

55

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

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

    int rStartEdgeDofs, nOverallEdgeDofs;
67
    mpi::getDofNumbering(domainComm, nEdgeDofs, 
68
69
70
			 rStartEdgeDofs, nOverallEdgeDofs);

    int rStartVertexDofs, nOverallVertexDofs;
71
    mpi::getDofNumbering(domainComm, nVertexDofs, 
72
73
74
75
76
77
			 rStartVertexDofs, nOverallVertexDofs);

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

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


    std::set<DegreeOfFreedom> otherBoundaryLocalDofs;
105
    for (DofComm::Iterator it(meshDistributor->getDofComm(0).getRecvDofs(), feSpace);
106
107
108
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
	otherBoundaryLocalDofs.insert(it.getDofIndex());
Thomas Witkowski's avatar
Thomas Witkowski committed
109
110
      
    interiorDofs.clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
111
112
113

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

128
129
130
131
132
133
134
    {
      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
135
136
137
138
139


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


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

      for (; !it.endDofIter(); it.nextDof()) {
147
	int globalSendDof = (*interiorMap)[feSpace][it.getDofIndex()].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
148
149
150
151

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

152
	stdMpi.getSendData(it.getRank()).push_back(mapGlobalBoundaryDof[globalSendDof]);
Thomas Witkowski's avatar
Thomas Witkowski committed
153
      }
154
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
155
156

    stdMpi.updateSendDataSize();
157

158
    for (DofComm::Iterator it(meshDistributor->getDofComm(0).getRecvDofs(), feSpace);
159
160
	 !it.end(); it.nextRank())
      stdMpi.recv(it.getRank());
Thomas Witkowski's avatar
Thomas Witkowski committed
161

162
    stdMpi.startCommunication();
163

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


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

176
    ISCreateStride(domainComm, 
Thomas Witkowski's avatar
Thomas Witkowski committed
177
178
179
180
		   nInteriorDofs * nComponents,
		   (rStartInteriorDofs + rStartBoundaryDofs) * nComponents, 
		   1, &interiorIs);

181
    ISCreateStride(domainComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
182
183
184
185
186
187
		   nBoundaryDofs * nComponents,
		   (rStartInteriorDofs + rStartBoundaryDofs + nInteriorDofs) * nComponents, 
		   1, &boundaryIs);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
188
  void PetscSolverSchur::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
189
190
  {
    FUNCNAME("PetscSolverSchur::fillPetscMatrix()");
191

192
    createMatVec(*seqMat);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
193

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

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


204
    MatCreateAIJ(domainComm, 
205
206
207
208
		 nInteriorRows, nInteriorRows, 
		 nOverallInteriorRows, nOverallInteriorRows,
		 100, PETSC_NULL, 100, PETSC_NULL, &matA11);

209
    MatCreateAIJ(domainComm, 
210
211
212
213
		 nBoundaryRows, nBoundaryRows, 
		 nOverallBoundaryRows, nOverallBoundaryRows,
		 100, PETSC_NULL, 100, PETSC_NULL, &matA22);

214
    MatCreateAIJ(domainComm, 
215
216
217
218
		 nInteriorRows, nBoundaryRows, 
		 nOverallInteriorRows, nOverallBoundaryRows,
		 100, PETSC_NULL, 100, PETSC_NULL, &matA12);    

219
    MatCreateAIJ(domainComm, 
220
221
222
		 nBoundaryRows, nInteriorRows, 
		 nOverallBoundaryRows, nOverallInteriorRows,
		 100, PETSC_NULL, 100, PETSC_NULL, &matA21);
Thomas Witkowski's avatar
Thomas Witkowski committed
223
224
225
226


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

    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;

252
    MatCreateNest(domainComm, 2, &tmpIS[0], 2, &tmpIS[0], &tmpMat[0][0], 
253
254
		  &getMatInterior());
    MatNestSetVecType(getMatInterior(), VECNEST);
Thomas Witkowski's avatar
Thomas Witkowski committed
255

256
    matAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
257

258
259
    int nRankRows = (*interiorMap)[feSpace].nRankDofs * nComponents;
    int nOverallRows = (*interiorMap)[feSpace].nOverallDofs * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
260

261
    VecCreateMPI(domainComm, nRankRows, nOverallRows, &petscSolVec);
262
263
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
264

265
266
267
  void PetscSolverSchur::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverSchur::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
268

269
    const FiniteElemSpace *feSpace = componentSpaces[0];
270
    int nComponents = vec->getSize();
271
272
    int nRankRows = (*interiorMap)[feSpace].nRankDofs * nComponents;
    int nOverallRows = (*interiorMap)[feSpace].nOverallDofs * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
273

274
    VecCreateMPI(domainComm, nRankRows, nOverallRows, &(getVecRhsInterior()));
Thomas Witkowski's avatar
Thomas Witkowski committed
275
276

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

279
    vecRhsAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
282
  }


283
284
  void PetscSolverSchur::solvePetscMatrix(SystemVector &vec, 
					  AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
287
  {
    FUNCNAME("PetscSolverSchur::solvePetscMatrix()");

288
    const FiniteElemSpace *feSpace = componentSpaces[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
289
290
    int nComponents = vec.getSize();

291
    KSPCreate(domainComm, &kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
292

293
294
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(), 
		    SAME_NONZERO_PATTERN); 
295
296
    KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetFromOptions(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
297

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


305
    KSPSolve(kspInterior, getVecRhsInterior(), petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
306
307
308
309
310
311
312
313
314

    // === 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) {
315
	DegreeOfFreedom globalRowDof = (*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
	if (boundaryDofs.count(globalRowDof)) {
	  int index =
	    (mapGlobalBoundaryDof[globalRowDof] - rStartBoundaryDofs + nInteriorDofs) * (i + 1);
	  *dofIt = vecPointer[index];
	} else {
	  int index = 
	    (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. ===

338
339
340
341
    MatDestroy(&matA11);
    MatDestroy(&matA12);
    MatDestroy(&matA21);
    MatDestroy(&matA22);
342

343
344
    matDestroy();
    vecDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
345

346
    KSPDestroy(&kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
349
  }
  

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
350
  void PetscSolverSchur::setDofMatrix(DOFMatrix* seqMat, int dispMult, 
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353
354
				      int dispAddRow, int dispAddCol)
  {
    FUNCNAME("PetscSolverSchur::setDofMatrix()");

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

357
    const FiniteElemSpace* feSpace = componentSpaces[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
358
359
360
361
    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
362
363
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
366
367
368
369
370
371
372
373
374

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

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

      // Global index of the current row DOF.
379
      int globalRowDof = (*interiorMap)[feSpace][cursor.value()].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
380
381
382
383
384
385
386
387
     
      colsBoundary.clear();
      colsInterior.clear();
      valuesBoundary.clear();
      valuesInterior.clear();
      
      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	   icursor != icend; ++icursor) {
388
	int globalColDof = (*interiorMap)[feSpace][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442

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

	  int colIndex = 
	    mapGlobalBoundaryDof[globalColDof] * dispMult + dispAddCol;
	  
	  colsBoundary.push_back(colIndex);
	  valuesBoundary.push_back(value(*icursor));
	} else {
	  TEST_EXIT_DBG(mapGlobalInteriorDof.count(globalColDof))
	    ("Cannot find global interior mapping for global column DOF %d!\n", 
	     globalColDof);

	  int colIndex = 
	    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");

	int rowIndex = 
	  mapGlobalBoundaryDof[globalRowDof] * dispMult + dispAddRow;

 	MatSetValues(matA22, 1, &rowIndex, colsBoundary.size(), 
 		     &(colsBoundary[0]), &(valuesBoundary[0]), ADD_VALUES);
 	MatSetValues(matA21, 1, &rowIndex, colsInterior.size(), 
 		     &(colsInterior[0]), &(valuesInterior[0]), ADD_VALUES);
      } else {
	TEST_EXIT_DBG(mapGlobalInteriorDof.count(globalRowDof))
	  ("Cannot find global interior mapping for global row DOF %d!\n", 
	   globalRowDof);

	int rowIndex = 
	  mapGlobalInteriorDof[globalRowDof] * dispMult + dispAddRow;

  	MatSetValues(matA11, 1, &rowIndex, colsInterior.size(), 
  		     &(colsInterior[0]), &(valuesInterior[0]), ADD_VALUES);
  	MatSetValues(matA12, 1, &rowIndex, colsBoundary.size(), 
  		     &(colsBoundary[0]), &(valuesBoundary[0]), ADD_VALUES);
      }
    }
  }


  void PetscSolverSchur::setDofVector(Vec& petscVec, DOFVector<double>* vec,
				      int dispMult, int dispAdd, bool rankOnly)
  {
443
    const FiniteElemSpace *feSpace = componentSpaces[0];
444

Thomas Witkowski's avatar
Thomas Witkowski committed
445
446
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
447
      if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
448
449
	continue;

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

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

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

476
} } // end namespace Parallel, AMDiS