From aa82de2b6042bc130f99335af0a7d6e1e124f50d Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Wed, 24 Oct 2012 10:53:15 +0000
Subject: [PATCH] Fixed problem for Dirichlet boundary conditions on non
diagonal blocks in parallel.
---
AMDiS/src/DirichletBC.h | 6 +--
AMDiS/src/parallel/PetscSolverGlobalMatrix.cc | 54 +++++++++++++++----
2 files changed, 45 insertions(+), 15 deletions(-)
diff --git a/AMDiS/src/DirichletBC.h b/AMDiS/src/DirichletBC.h
index 2babebca..1773e229 100644
--- a/AMDiS/src/DirichletBC.h
+++ b/AMDiS/src/DirichletBC.h
@@ -108,10 +108,8 @@ namespace AMDiS {
/// DOFVector containing the boundary values
DOFVectorBase<double> *dofVec;
- /** \brief
- * Defines, if the boundary condition must be applied to the matrix. See
- * comment of \ref BoundaryCondition::applyBoundaryCondition.
- */
+ /// Defines, if the boundary condition must be applied to the matrix. See
+ /// comment of \ref BoundaryCondition::applyBoundaryCondition.
bool applyBC;
};
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index d4cdceea..d413a2e6 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -476,33 +476,65 @@ namespace AMDiS {
return;
int nComponents = seqMat->getSize();
- vector<PetscInt> zeroRows;
+ vector<int> dirichletRows;
+ vector<int> dirichletCols;
+ vector<double> dirichletValues;
for (int i = 0; i < nComponents; i++) {
- DOFMatrix *dofMat = (*seqMat)[i][i];
- if (!dofMat)
+ bool dirichletRow = false;
+ int dirichletMainCol = -1;
+
+ for (int j = 0; j < nComponents; j++) {
+ DOFMatrix *dofMat = (*seqMat)[i][j];
+ if (!dofMat)
+ continue;
+
+ BoundaryIndexMap& bounds =
+ const_cast<BoundaryIndexMap&>(dofMat->getBoundaryManager()->getBoundaryConditionMap());
+
+ for (BoundaryIndexMap::iterator bIt = bounds.begin(); bIt != bounds.end(); ++bIt) {
+ if (bIt->second && bIt->second->isDirichlet()) {
+ dirichletRow = true;
+ if ((dynamic_cast<DirichletBC*>(bIt->second))->applyBoundaryCondition()) {
+ TEST_EXIT_DBG(dirichletMainCol == -1)("Should not happen!\n");
+ dirichletMainCol = j;
+ }
+ }
+ }
+ }
+
+ if (!dirichletRow)
continue;
+ DOFMatrix *dofMat = (*seqMat)[i][dirichletMainCol];
const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
+ TEST_EXIT(dofMat->getRowFeSpace() == dofMat->getColFeSpace())
+ ("I have to think about this scenario! Really possible?\n");
std::set<DegreeOfFreedom> &dRows = dofMat->getDirichletRows();
for (std::set<DegreeOfFreedom>::iterator dofIt = dRows.begin();
dofIt != dRows.end(); ++dofIt) {
MultiIndex multiIndex;
if (dofMap[feSpace].find(*dofIt, multiIndex) &&
- dofMap[feSpace].isRankDof(*dofIt))
- zeroRows.push_back(dofMap.getMatIndex(i, multiIndex.global));
+ dofMap[feSpace].isRankDof(*dofIt)) {
+ int rowIndex = dofMap.getMatIndex(i, multiIndex.global);
+ int colIndex = dofMap.getMatIndex(dirichletMainCol, multiIndex.global);
+ dirichletRows.push_back(rowIndex);
+ dirichletCols.push_back(colIndex);
+ dirichletValues.push_back(1.0);
+ }
}
}
- MatZeroRows(mpiMat, zeroRows.size(), &(zeroRows[0]), 1.0,
+ MatZeroRows(mpiMat, dirichletRows.size(), &(dirichletRows[0]), 1.0,
PETSC_NULL, PETSC_NULL);
- PetscViewer view;
- PetscViewerBinaryOpen(PETSC_COMM_WORLD, "interior.mat",
- FILE_MODE_WRITE, &view);
- MatView(mpiMat, view);
- PetscViewerDestroy(&view);
+ for (int i = 0; i < static_cast<int>(dirichletRows.size()); i++) {
+ MatSetValue(mpiMat, dirichletRows[i], dirichletCols[i], dirichletValues[i], INSERT_VALUES);
+ }
+
+ MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY);
}
--
GitLab