diff --git a/AMDiS/src/parallel/ParallelCoarseSpaceSolver.cc b/AMDiS/src/parallel/ParallelCoarseSpaceSolver.cc
index a7247f50e6c13329780477a4ae785c5b37a867a5..558dc46a5e08466b41fb4108a1b09be102f69bc1 100644
--- a/AMDiS/src/parallel/ParallelCoarseSpaceSolver.cc
+++ b/AMDiS/src/parallel/ParallelCoarseSpaceSolver.cc
@@ -154,14 +154,13 @@ namespace AMDiS {
MatCreateSeqAIJ(mpiCommLocal, nRankInteriorRows, nRankInteriorRows,
0, nnz[0][0].dnnz,
&mat[0][0]);
- MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
} else {
MatCreateAIJ(mpiCommGlobal, nRankInteriorRows, nRankInteriorRows,
nOverallInteriorRows, nOverallInteriorRows,
0, nnz[0][0].dnnz, 0, nnz[0][0].onnz,
- &mat[0][0]);
- MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
+ &mat[0][0]);
}
+ MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecSol[0]);
VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecRhs[0]);
@@ -178,6 +177,7 @@ namespace AMDiS {
nOverallCoarseRows, nOverallCoarseRows,
0, nnz[i + 1][i + 1].dnnz, 0, nnz[i + 1][i + 1].onnz,
&mat[i + 1][i + 1]);
+ MatSetOption(mat[i + 1][i + 1], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
cMap->createVec(vecSol[i + 1]);
cMap->createVec(vecRhs[i + 1]);
}
diff --git a/AMDiS/src/parallel/ParallelCoarseSpaceSolver.h b/AMDiS/src/parallel/ParallelCoarseSpaceSolver.h
index 452bb0c0eb6f10d5f8a3741a7d9c26b16e2bea2b..f4c21fd61313e9efc326dbe7aca8f94a88b50e49 100644
--- a/AMDiS/src/parallel/ParallelCoarseSpaceSolver.h
+++ b/AMDiS/src/parallel/ParallelCoarseSpaceSolver.h
@@ -234,6 +234,12 @@ namespace AMDiS {
}
+ /// Returns whether the solver has a coarse grid.
+ inline bool hasCoarseSpace()
+ {
+ return (!coarseSpaceMap.empty());
+ }
+
protected:
/// Prepare internal data structures. First, it create \ref uniqueCoarseMap
/// and \ref componentIthCoarseMap . Both are used to create the correct
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 6597dc2c57e437a568280c0d33d73de7e86c4a3f..5cdc1ce415ace8847bc8f527e1e66763ca548810 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -108,7 +108,8 @@ namespace AMDiS {
if (subdomain == NULL) {
subdomain = new PetscSolverGlobalMatrix("");
subdomain->setSymmetric(isSymmetric);
- subdomain->setHandleDirichletRows(false);
+ // subdomain->setHandleDirichletRows(false);
+ subdomain->setHandleDirichletRows(true);
if (meshLevel == 0) {
subdomain->setMeshDistributor(meshDistributor,
@@ -142,6 +143,8 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::createDirichletData()");
+ return;
+
int nComponents = mat.getSize();
for (int component = 0; component < nComponents; component++) {
DOFMatrix* dofMat = mat[component][component];
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 18e75fa99b630f5bf47064a046da8c90b1825bfe..a7b916a4c92d22cbf37a2978800610d1db43019d 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -17,6 +17,17 @@
namespace AMDiS {
+
+ PetscSolverGlobalMatrix::PetscSolverGlobalMatrix(string name)
+ : PetscSolver(name),
+ zeroStartVector(false),
+ printMatInfo(false)
+ {
+ Parameters::get("parallel->use zero start vector", zeroStartVector);
+ Parameters::get("parallel->print matrix info", printMatInfo);
+ }
+
+
void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
{
FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");
@@ -54,7 +65,7 @@ namespace AMDiS {
matAssembly();
- removeDirichletRows(seqMat, dofMap, getMatInterior());
+ removeDirichletRows(seqMat);
if (printMatInfo) {
MatInfo matInfo;
@@ -225,6 +236,8 @@ namespace AMDiS {
matAssembly();
+ removeDirichletRows(seqMat);
+
// === Create solver for the non primal (thus local) variables. ===
KSPCreate(mpiCommLocal, &kspInterior);
@@ -266,7 +279,7 @@ namespace AMDiS {
vecRhsAssembly();
- removeDirichletRows(vec, dofMap, getVecRhsInterior());
+ removeDirichletRows(vec);
// === For debugging allow to write the rhs vector to a file. ===
@@ -456,26 +469,25 @@ namespace AMDiS {
}
- void PetscSolverGlobalMatrix::removeDirichletRows(Matrix<DOFMatrix*> *seqMat,
- ParallelDofMapping &dofMap,
- Mat mpiMat)
+ void PetscSolverGlobalMatrix::removeDirichletRows(Matrix<DOFMatrix*> *seqMat)
{
FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()");
if (!handleDirichletRows)
return;
+ vector<int> dRowsInterior, dRowsCoarse;
+ vector<int> dColsInterior, dColsCoarse;
+ vector<double> dValuesInterior, dValuesCoarse;
+
int nComponents = seqMat->getSize();
- vector<int> dirichletRows;
- vector<int> dirichletCols;
- vector<double> dirichletValues;
- for (int i = 0; i < nComponents; i++) {
+ for (int rowComp = 0; rowComp < nComponents; rowComp++) {
bool dirichletRow = false;
int dirichletMainCol = -1;
- for (int j = 0; j < nComponents; j++) {
- DOFMatrix *dofMat = (*seqMat)[i][j];
+ for (int colComp = 0; colComp < nComponents; colComp++) {
+ DOFMatrix *dofMat = (*seqMat)[rowComp][colComp];
if (!dofMat)
continue;
@@ -486,7 +498,7 @@ namespace AMDiS {
if (bIt->second && bIt->second->isDirichlet()) {
dirichletRow = true;
if ((dynamic_cast<DirichletBC*>(bIt->second))->applyBoundaryCondition()) {
- dirichletMainCol = j;
+ dirichletMainCol = colComp;
break;
}
}
@@ -496,65 +508,152 @@ namespace AMDiS {
if (!dirichletRow)
continue;
- DOFMatrix *dofMat = (*seqMat)[i][dirichletMainCol];
- const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
+ DOFMatrix *dofMat = (*seqMat)[rowComp][dirichletMainCol];
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)) {
- 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);
+ if (hasCoarseSpace()) {
+ bool isRowCoarse = isCoarseSpace(rowComp, *dofIt);
+ bool isColCoarse = isCoarseSpace(dirichletMainCol, *dofIt);
+ TEST_EXIT(isRowCoarse == isColCoarse)
+ ("Really possible? So reimplement AMDiS from the scratch!\n");
+
+ if (isRowCoarse) {
+ ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComp];
+ ParallelDofMapping *colCoarseSpace = coarseSpaceMap[dirichletMainCol];
+
+ MultiIndex multiIndex;
+ if ((*rowCoarseSpace)[rowComp].find(*dofIt, multiIndex) &&
+ (*rowCoarseSpace)[rowComp].isRankDof(*dofIt)) {
+ int rowIndex = rowCoarseSpace->getLocalMatIndex(rowComp, *dofIt);
+ int colIndex = colCoarseSpace->getLocalMatIndex(dirichletMainCol, *dofIt);
+ dRowsCoarse.push_back(rowIndex);
+ dColsCoarse.push_back(colIndex);
+ dValuesCoarse.push_back(1.0);
+ }
+ } else {
+
+ MultiIndex multiIndex;
+ if ((*interiorMap)[rowComp].find(*dofIt, multiIndex) &&
+ (*interiorMap)[rowComp].isRankDof(*dofIt)) {
+ int rowIndex = interiorMap->getLocalMatIndex(rowComp, *dofIt);
+ int colIndex = interiorMap->getLocalMatIndex(dirichletMainCol, *dofIt);
+ dRowsInterior.push_back(rowIndex);
+ dColsInterior.push_back(colIndex);
+ dValuesInterior.push_back(1.0);
+ }
+ }
+ } else {
+ MultiIndex multiIndex;
+ if ((*interiorMap)[rowComp].find(*dofIt, multiIndex) &&
+ (*interiorMap)[rowComp].isRankDof(*dofIt)) {
+ int rowIndex = interiorMap->getMatIndex(rowComp, multiIndex.global);
+ int colIndex = interiorMap->getMatIndex(dirichletMainCol, multiIndex.global);
+ dRowsInterior.push_back(rowIndex);
+ dColsInterior.push_back(colIndex);
+ dValuesInterior.push_back(1.0);
+ }
}
}
}
- MatZeroRows(mpiMat, dirichletRows.size(), &(dirichletRows[0]), 0.0,
- PETSC_NULL, PETSC_NULL);
+ {
+ Mat mpiMat = getMatInterior();
+ MatZeroRows(mpiMat, dRowsInterior.size(), &(dRowsInterior[0]), 0.0,
+ PETSC_NULL, PETSC_NULL);
+
+ for (int i = 0; i < static_cast<int>(dRowsInterior.size()); i++)
+ MatSetValue(mpiMat, dRowsInterior[i], dColsInterior[i],
+ dValuesInterior[i], INSERT_VALUES);
+
+ MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY);
+ }
- 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);
+ if (hasCoarseSpace()) {
+ Mat mpiMat = getMatCoarse();
+ MatZeroRows(mpiMat, dRowsCoarse.size(), &(dRowsCoarse[0]), 0.0,
+ PETSC_NULL, PETSC_NULL);
+
+ for (int i = 0; i < static_cast<int>(dRowsCoarse.size()); i++)
+ MatSetValue(mpiMat, dRowsCoarse[i], dColsCoarse[i],
+ dValuesCoarse[i], INSERT_VALUES);
+
+ MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY);
+ }
}
- void PetscSolverGlobalMatrix::removeDirichletRows(SystemVector *seqVec,
- ParallelDofMapping &dofMap,
- Vec mpiVec)
+ void PetscSolverGlobalMatrix::removeDirichletRows(SystemVector *seqVec)
{
FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()");
if (!handleDirichletRows)
return;
- int nComponents = seqVec->getSize();
-
- for (int i = 0; i < nComponents; i++) {
- DOFVector<double> *dofVec = seqVec->getDOFVector(i);
- const FiniteElemSpace *feSpace = dofVec->getFeSpace();
+ int cInterior = 0;
+ int cCoarse = 0;
+ int nComponents = seqVec->getSize();
+ for (int component = 0; component < nComponents; component++) {
+ DOFVector<double> *dofVec = seqVec->getDOFVector(component);
map<DegreeOfFreedom, double>& dValues = dofVec->getDirichletValues();
+
for (map<DegreeOfFreedom, double>::iterator dofIt = dValues.begin();
dofIt != dValues.end(); ++dofIt) {
- MultiIndex multiIndex;
- if (dofMap[feSpace].find(dofIt->first, multiIndex) &&
- dofMap[feSpace].isRankDof(dofIt->first))
- VecSetValue(mpiVec, dofMap.getMatIndex(i, multiIndex.global),
- dofIt->second, INSERT_VALUES);
+ if (hasCoarseSpace()) {
+ if (isCoarseSpace(component, dofIt->first)) {
+ ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[component];
+ MultiIndex multiIndex;
+ if ((*rowCoarseSpace)[component].find(dofIt->first, multiIndex) &&
+ (*rowCoarseSpace)[component].isRankDof(dofIt->first)) {
+ if (dofIt->second > 0.5)
+ cCoarse++;
+ VecSetValue(getVecRhsCoarse(),
+ rowCoarseSpace->getMatIndex(component, multiIndex.global),
+ dofIt->second, INSERT_VALUES);
+ }
+ } else {
+ MultiIndex multiIndex;
+ if ((*interiorMap)[component].find(dofIt->first, multiIndex) &&
+ (*interiorMap)[component].isRankDof(dofIt->first)) {
+ if (dofIt->second > 0.5)
+ cInterior++;
+ VecSetValue(getVecRhsInterior(),
+ interiorMap->getLocalMatIndex(component, dofIt->first),
+ dofIt->second, INSERT_VALUES);
+ }
+ }
+ } else {
+ MultiIndex multiIndex;
+ if ((*interiorMap)[component].find(dofIt->first, multiIndex) &&
+ (*interiorMap)[component].isRankDof(dofIt->first)) {
+ if (dofIt->second > 0.5)
+ cInterior++;
+ VecSetValue(getVecRhsInterior(),
+ interiorMap->getMatIndex(component, multiIndex.global),
+ dofIt->second, INSERT_VALUES);
+ }
+ }
}
}
- VecAssemblyBegin(mpiVec);
- VecAssemblyEnd(mpiVec);
+ VecAssemblyBegin(getVecRhsInterior());
+ VecAssemblyEnd(getVecRhsInterior());
+
+ if (hasCoarseSpace()) {
+ VecAssemblyBegin(getVecRhsCoarse());
+ VecAssemblyEnd(getVecRhsCoarse());
+ }
+
+ mpi::globalAdd(cInterior);
+ mpi::globalAdd(cCoarse);
+ MSG("WORKED ON DIRICHLET DOFS: %d %d\n", cInterior, cCoarse);
}
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
index b3ff065d8b4a7812c570a05739c0cfdf226ae9ec..6ae8a3343eef3101ffd465cdceaca25a7c69782f 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
@@ -39,14 +39,7 @@ namespace AMDiS {
class PetscSolverGlobalMatrix : public PetscSolver
{
public:
- PetscSolverGlobalMatrix(string name)
- : PetscSolver(name),
- zeroStartVector(false),
- printMatInfo(false)
- {
- Parameters::get("parallel->use zero start vector", zeroStartVector);
- Parameters::get("parallel->print matrix info", printMatInfo);
- }
+ PetscSolverGlobalMatrix(string name);
void fillPetscMatrix(Matrix<DOFMatrix*> *mat);
@@ -63,13 +56,9 @@ namespace AMDiS {
void destroyVectorData();
protected:
- void removeDirichletRows(Matrix<DOFMatrix*> *seqMat,
- ParallelDofMapping &dofMap,
- Mat mpiMat);
+ void removeDirichletRows(Matrix<DOFMatrix*> *seqMat);
- void removeDirichletRows(SystemVector *seqVec,
- ParallelDofMapping &dofMap,
- Vec mpiVec);
+ void removeDirichletRows(SystemVector *seqVec);
/// Reads field split information and creats a splitting based on
/// component numbers.