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