From 39e4e6880d785bd3c3fdb6342d1ae66511a0115d Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 5 Mar 2012 14:48:10 +0000
Subject: [PATCH] SHIT, SHIT, SHIT! JUST FOR JADE....
---
AMDiS/src/parallel/PetscSolverFeti.cc | 335 +++++---------------------
AMDiS/src/parallel/PetscSolverFeti.h | 19 +-
2 files changed, 75 insertions(+), 279 deletions(-)
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index 2c2fb91e..4e945523 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -221,10 +221,9 @@ namespace AMDiS {
createPrimals(feSpace);
createDuals(feSpace);
+ createLagrange(feSpace);
createIndexB(feSpace);
}
-
- createLagrange(meshDistributor->getFeSpaces());
}
@@ -324,8 +323,7 @@ namespace AMDiS {
for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
!it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof()) {
- // If DOF is not primal, i.e., its a dual node
- if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
+ if (!isPrimal(feSpace, it.getDofIndex())) {
boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
}
@@ -340,7 +338,7 @@ namespace AMDiS {
for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
!it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof())
- if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
+ if (!isPrimal(feSpace, it.getDofIndex()))
stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
stdMpi.updateSendDataSize();
@@ -349,7 +347,7 @@ namespace AMDiS {
!it.end(); it.nextRank()) {
bool recvFromRank = false;
for (; !it.endDofIter(); it.nextDof()) {
- if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
+ if (!isPrimal(feSpace, it.getDofIndex())) {
recvFromRank = true;
break;
}
@@ -365,7 +363,7 @@ namespace AMDiS {
!it.end(); it.nextRank()) {
int i = 0;
for (; !it.endDofIter(); it.nextDof())
- if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
+ if (!isPrimal(feSpace, it.getDofIndex()))
boundaryDofRanks[it.getDofIndex()] =
stdMpi.getRecvData(it.getRank())[i++];
}
@@ -377,25 +375,24 @@ namespace AMDiS {
for (DofContainer::iterator it = allBoundaryDofs.begin();
it != allBoundaryDofs.end(); ++it)
- if (primalDofMap[feSpace].isSet(**it) == false)
+ if (!isPrimal(feSpace, **it))
dualDofMap[feSpace].insertRankDof(**it);
dualDofMap[feSpace].update(false);
MSG("nRankDuals = %d nOverallDuals = %d\n",
- dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs);
+ dualDofMap[feSpace].nRankDofs,
+ dualDofMap[feSpace].nOverallDofs);
}
- void PetscSolverFeti::createLagrange(vector<const FiniteElemSpace*> &feSpaces)
+ void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
{
FUNCNAME("PetscSolverFeti::createLagrange()");
// === Reserve for each dual node, on the rank that owns this node, the ===
// === appropriate number of Lagrange constraints. ===
-
-#if 0
int nRankLagrange = 0;
DofMapping& dualMap = dualDofMap[feSpace].getMap();
for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
@@ -406,10 +403,11 @@ namespace AMDiS {
}
}
lagrangeMap[feSpace].nRankDofs = nRankLagrange;
- lagrangeMap[feSpace].update(false);
+ lagrangeMap[feSpace].update();
MSG("nRankLagrange = %d nOverallLagrange = %d\n",
- lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs);
+ lagrangeMap[feSpace].nRankDofs,
+ lagrangeMap[feSpace].nOverallDofs);
// === Communicate lagrangeMap to all other ranks. ===
@@ -419,7 +417,7 @@ namespace AMDiS {
for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
!it.end(); it.nextRank()) {
for (; !it.endDofIter(); it.nextDof())
- if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
+ if (!isPrimal(feSpace, it.getDofIndex())) {
DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
stdMpi.getSendData(it.getRank()).push_back(d);
}
@@ -431,7 +429,7 @@ namespace AMDiS {
!it.end(); it.nextRank()) {
bool recvData = false;
for (; !it.endDofIter(); it.nextDof())
- if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
+ if (!isPrimal(feSpace, it.getDofIndex())) {
recvData = true;
break;
}
@@ -446,13 +444,12 @@ namespace AMDiS {
!it.end(); it.nextRank()) {
int counter = 0;
for (; !it.endDofIter(); it.nextDof()) {
- if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
+ if (!isPrimal(feSpace, it.getDofIndex())) {
DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
lagrangeMap[feSpace].insert(it.getDofIndex(), d);
}
}
}
-#endif
}
@@ -469,7 +466,7 @@ namespace AMDiS {
nLocalInterior = 0;
for (int i = 0; i < admin->getUsedSize(); i++) {
if (admin->isDofFree(i) == false &&
- primalDofMap[feSpace].isSet(i) == false &&
+ isPrimal(feSpace, i) == false &&
dualDofMap[feSpace].isSet(i) == false) {
localDofMap[feSpace].insertRankDof(i);
nLocalInterior++;
@@ -492,10 +489,8 @@ namespace AMDiS {
localDofMap[feSpace].update(false);
-#if 0
dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs +
nLocalInterior);
-#endif
}
@@ -503,13 +498,13 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::createMatLagrange()");
+ const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
+
// === Create distributed matrix for Lagrange constraints. ===
MatCreateMPIAIJ(PETSC_COMM_WORLD,
- lagrangeMap.getRankDofs(feSpaces),
- localDofMap.getRankDofs(),
- lagrangeMap.getOverallDofs(feSpaces),
- localDofMap.getOverallDofs(),
+ lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
+ lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),
2, PETSC_NULL, 2, PETSC_NULL,
&mat_lagrange);
@@ -520,38 +515,33 @@ namespace AMDiS {
// === m == r, than the rank sets -1.0 for the corresponding ===
// === constraint. ===
- for (unsigned int feIdx = 0; feIdx < feSpaces.size(); feIdx++) {
- const FiniteElemSpace *feSpace = feSpaces[feIdx];
-
- DofMapping &dualMap = dualDofMap[feSpace].getMap();
- for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
- TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
- ("Should not happen!\n");
-
- // Global index of the first Lagrange constriant for this node.
- int constraintMatIndex = lagrangeMap.mapGlobal(it->first, feIdx);
- // Copy set of all ranks that contain this dual node.
- vector<int> W(boundaryDofRanks[it->first].begin(),
- boundaryDofRanks[it->first].end());
- // Number of ranks that contain this dual node.
- int degree = W.size();
-
- for (int i = 0; i < degree; i++) {
- for (int j = i + 1; j < degree; j++) {
- if (W[i] == mpiRank || W[j] == mpiRank) {
-#if 1
- // Get global index of the dual variable
- int colIndex = localDofMap.mapGlobal(it->first, feIdx);
-#else
+ DofMapping &dualMap = dualDofMap[feSpace].getMap();
+ for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
+ TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
+ ("Should not happen!\n");
+
+ // Global index of the first Lagrange constriant for this node.
+ int index = lagrangeMap[feSpace][it->first];
+ // Copy set of all ranks that contain this dual node.
+ vector<int> W(boundaryDofRanks[it->first].begin(),
+ boundaryDofRanks[it->first].end());
+ // Number of ranks that contain this dual node.
+ int degree = W.size();
+
+ for (int i = 0; i < degree; i++) {
+ for (int j = i + 1; j < degree; j++) {
+ if (W[i] == mpiRank || W[j] == mpiRank) {
+ // Set the constraint for all components of the system.
+ for (int k = 0; k < nComponents; k++) {
+ int rowIndex = index * nComponents + k;
int colIndex = it->second * nComponents + k;
-#endif
double value = (W[i] == mpiRank ? 1.0 : -1.0);
- MatSetValue(mat_lagrange, constraintMatIndex, colIndex, value,
- INSERT_VALUES);
+ MatSetValue(mat_lagrange, rowIndex, colIndex, value,
+ INSERT_VALUES);
}
-
- constraintMatIndex++;
}
+
+ index++;
}
}
}
@@ -579,12 +569,12 @@ namespace AMDiS {
localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
&(schurPrimalData.tmp_vec_b));
VecCreateMPI(PETSC_COMM_WORLD,
- primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
+ primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
&(schurPrimalData.tmp_vec_primal));
MatCreateShell(PETSC_COMM_WORLD,
- primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nRankDofs * nComponents,
- primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
+ primalDofMap.getRankDofs(), primalDofMap.getRankDofs(),
+ primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
&schurPrimalData,
&mat_schur_primal);
MatShellSetOperation(mat_schur_primal, MATOP_MULT,
@@ -601,8 +591,8 @@ namespace AMDiS {
double wtime = MPI::Wtime();
- int nRowsRankPrimal = primalDofMap[feSpace].nRankDofs * nComponents;
- int nRowsOverallPrimal = primalDofMap[feSpace].nOverallDofs * nComponents;
+ int nRowsRankPrimal = primalDofMap.getRankDofs();
+ int nRowsOverallPrimal = primalDofMap.getOverallDofs();
int nRowsRankB = localDofMap.getRankDofs();
int nRowsOverallB = localDofMap.getOverallDofs();
@@ -634,7 +624,8 @@ namespace AMDiS {
mpi::globalMax(maxLocalPrimal);
TEST_EXIT(mat_b_primal_cols.size() ==
- primalDofMap[feSpace].size() * nComponents)("Should not happen!\n");
+ primalDofMap[feSpace].size() * nComponents)
+ ("Should not happen!\n");
for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
it != mat_b_primal_cols.end(); ++it) {
Vec tmpVec;
@@ -711,8 +702,6 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::createFetiKsp()");
- const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
-
// === Create FETI-DP solver object. ===
fetiData.mat_primal_b = &mat_primal_b;
@@ -725,19 +714,15 @@ namespace AMDiS {
localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
&(fetiData.tmp_vec_b));
VecCreateMPI(PETSC_COMM_WORLD,
- lagrangeMap[feSpace].nRankDofs * nComponents,
- lagrangeMap[feSpace].nOverallDofs * nComponents,
+ lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
&(fetiData.tmp_vec_lagrange));
VecCreateMPI(PETSC_COMM_WORLD,
- primalDofMap[feSpace].nRankDofs * nComponents,
- primalDofMap[feSpace].nOverallDofs * nComponents,
+ primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
&(fetiData.tmp_vec_primal));
MatCreateShell(PETSC_COMM_WORLD,
- lagrangeMap[feSpace].nRankDofs * nComponents,
- lagrangeMap[feSpace].nRankDofs * nComponents,
- lagrangeMap[feSpace].nOverallDofs * nComponents,
- lagrangeMap[feSpace].nOverallDofs * nComponents,
+ lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
+ lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
&fetiData, &mat_feti);
MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);
@@ -990,9 +975,9 @@ namespace AMDiS {
int nRowsRankB = localDofMap.getRankDofs();
int nRowsOverallB = localDofMap.getOverallDofs();
- int nRowsRankPrimal = primalDofMap.getRankDofs(feSpaces);
- int nRowsOverallPrimal = primalDofMap.getOverallDofs(feSpaces);
- int nRowsDual = dualDofMap.getRankDofs(feSpaces);
+ int nRowsRankPrimal = primalDofMap.getRankDofs();
+ int nRowsOverallPrimal = primalDofMap.getOverallDofs();
+ int nRowsDual = dualDofMap.getRankDofs();
int nRowsInterior = nRowsRankB - nRowsDual;
MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
@@ -1341,14 +1326,14 @@ namespace AMDiS {
VecCreateMPI(PETSC_COMM_WORLD,
localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &f_b);
VecCreateMPI(PETSC_COMM_WORLD,
- primalDofMap[feSpace].nRankDofs * nComponents,
- primalDofMap[feSpace].nOverallDofs * nComponents, &f_primal);
+ primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
+ &f_primal);
for (int i = 0; i < nComponents; i++) {
DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
int index = dofIt.getDOFIndex();
- if (primalDofMap[feSpace].isSet(index)) {
+ if (isPrimal(feSpace, index)) {
index = primalDofMap[feSpace][index] * nComponents + i;
double value = *dofIt;
VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
@@ -1408,203 +1393,9 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
-#if 0
- const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
-
- int nOverallNest = localDofMap.getOverallDofs() +
- (primalDofMap[feSpace].nOverallDofs +
- lagrangeMap[feSpace].nOverallDofs) * nComponents;
-
- Mat mat_lagrange_transpose;
- MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose);
-
- Mat A_global;
- MatCreateMPIAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,
- nOverallNest, nOverallNest, 30, PETSC_NULL, 30, PETSC_NULL,
- &A_global);
- Vec A_global_rhs, A_global_sol;
- MatGetVecs(A_global, &A_global_rhs, &A_global_sol);
-
- PetscInt nCols;
- const PetscInt *cols;
- const PetscScalar *vals;
-
- // === mat_b_b ===
- for (int i = 0; i < localDofMap[feSpace].nRankDofs * nComponents; i++) {
- int rowIndex = localDofMap[feSpace].rStartDofs * nComponents + i;
- MatGetRow(mat_b_b, i, &nCols, &cols, &vals);
- MatSetValues(A_global, 1, &rowIndex, nCols, cols, vals, ADD_VALUES);
- MatRestoreRow(mat_b_b, rowIndex, &nCols, &cols, &vals);
- }
-
- PetscScalar *g_f_b;
- VecGetArray(f_b, &g_f_b);
- for (int i = 0; i < localDofMap[feSpace].nRankDofs * nComponents; i++)
- VecSetValue(A_global_rhs, localDofMap[feSpace].rStartDofs * nComponents + i, g_f_b[i], INSERT_VALUES);
- VecRestoreArray(f_b, &g_f_b);
-
-
-
- // === mat_primal_primal ===
-
- for (int i = 0; i < primalDofMap[feSpace].nRankDofs * nComponents; i++) {
- int rowIndex = primalDofMap[feSpace].rStartDofs * nComponents + i;
- MatGetRow(mat_primal_primal, rowIndex, &nCols, &cols, &vals);
-
- int rowIndexA = localDofMap[feSpace].nOverallDofs * nComponents + rowIndex;
- vector<int> colsA(nCols);
- for (int j = 0; j < nCols; j++)
- colsA[j] = localDofMap[feSpace].nOverallDofs * nComponents + cols[j];
-
- MatSetValues(A_global, 1, &rowIndexA, nCols, &(colsA[0]), vals, ADD_VALUES);
- MatRestoreRow(mat_primal_primal, rowIndex, &nCols, &cols, &vals);
- }
-
- PetscScalar *g_f_primal;
- VecGetArray(f_primal, &g_f_primal);
- for (int i = 0; i < primalDofMap[feSpace].nRankDofs * nComponents; i++)
- VecSetValue(A_global_rhs,
- (localDofMap[feSpace].nOverallDofs + primalDofMap[feSpace].rStartDofs) * nComponents + i, g_f_primal[i],
- INSERT_VALUES);
- VecRestoreArray(f_primal, &g_f_primal);
-
-
- // === mat_b_primal ===
-
- for (int i = 0; i < localDofMap[feSpace].nRankDofs * nComponents; i++) {
- int rowIndex = localDofMap[feSpace].rStartDofs * nComponents + i;
- MatGetRow(mat_b_primal, rowIndex, &nCols, &cols, &vals);
-
- vector<int> colsA(nCols);
- for (int j = 0; j < nCols; j++)
- colsA[j] = localDofMap[feSpace].nOverallDofs * nComponents + cols[j];
-
- MatSetValues(A_global, 1, &rowIndex, nCols, &(colsA[0]), vals, ADD_VALUES);
- MatRestoreRow(mat_b_primal, rowIndex, &nCols, &cols, &vals);
- }
-
- // === mat_primal_b ===
-
- for (int i = 0; i < primalDofMap[feSpace].nRankDofs * nComponents; i++) {
- int rowIndex = primalDofMap[feSpace].rStartDofs * nComponents + i;
- MatGetRow(mat_primal_b, rowIndex, &nCols, &cols, &vals);
-
- int rowIndexA = localDofMap[feSpace].nOverallDofs * nComponents + rowIndex;
-
- MatSetValues(A_global, 1, &rowIndexA, nCols, cols, vals, ADD_VALUES);
- MatRestoreRow(mat_primal_b, rowIndex, &nCols, &cols, &vals);
- }
-
- // === mat_lagrange ===
-
- for (int i = 0; i < lagrangeMap[feSpace].nRankDofs * nComponents; i++) {
- int rowIndex = lagrangeMap[feSpace].rStartDofs * nComponents + i;
- MatGetRow(mat_lagrange, rowIndex, &nCols, &cols, &vals);
-
- int rowIndexA = (localDofMap[feSpace].nOverallDofs + primalDofMap[feSpace].nOverallDofs) * nComponents + rowIndex;
-
- MatSetValues(A_global, 1, &rowIndexA, nCols, cols, vals, ADD_VALUES);
- MatRestoreRow(mat_lagrange, rowIndex, &nCols, &cols, &vals);
- }
-
- // === mat_lagrange_transpose ===
-
- for (int i = 0; i < localDofMap[feSpace].nRankDofs * nComponents; i++) {
- int rowIndex = localDofMap[feSpace].rStartDofs * nComponents + i;
- MatGetRow(mat_lagrange_transpose, rowIndex, &nCols, &cols, &vals);
-
- vector<int> colsA(nCols);
- for (int j = 0; j < nCols; j++)
- colsA[j] = (localDofMap[feSpace].nOverallDofs + primalDofMap[feSpace].nOverallDofs) * nComponents + cols[j];
-
- MatSetValues(A_global, 1, &rowIndex, nCols, &(colsA[0]), vals, ADD_VALUES);
- MatRestoreRow(mat_lagrange_transpose, rowIndex, &nCols, &cols, &vals);
- }
-
- MatAssemblyBegin(A_global, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(A_global, MAT_FINAL_ASSEMBLY);
-
- VecAssemblyBegin(A_global_rhs);
- VecAssemblyEnd(A_global_rhs);
-
- // === Create solver and solve the overall FETI system. ===
-
- KSP ksp;
- KSPCreate(PETSC_COMM_WORLD, &ksp);
- KSPSetOperators(ksp, A_global, A_global, SAME_NONZERO_PATTERN);
- KSPSetFromOptions(ksp);
- KSPSolve(ksp, A_global_rhs, A_global_sol);
-
-
- Vec u_b, u_primal;
- VecDuplicate(f_b, &u_b);
- VecDuplicate(f_primal, &u_primal);
-
- vector<PetscInt> localBIndex, globalBIndex;
- localBIndex.reserve(localDofMap[feSpace].nRankDofs * nComponents);
- globalBIndex.reserve(localDofMap[feSpace].nRankDofs * nComponents);
-
- for (int i = 0; i < localDofMap[feSpace].nRankDofs * nComponents; i++) {
- localBIndex.push_back(localDofMap[feSpace].rStartDofs * nComponents + i);
- globalBIndex.push_back(localDofMap[feSpace].rStartDofs * nComponents + i);
- }
- IS localBIs, globalBIs;
- ISCreateGeneral(PETSC_COMM_SELF,
- localBIndex.size(),
- &(localBIndex[0]),
- PETSC_USE_POINTER,
- &localBIs);
- ISCreateGeneral(PETSC_COMM_SELF,
- globalBIndex.size(),
- &(globalBIndex[0]),
- PETSC_USE_POINTER,
- &globalBIs);
-
- VecScatter vecscat;
- VecScatterCreate(A_global_sol, globalBIs, u_b, localBIs, &vecscat);
- VecScatterBegin(vecscat, A_global_sol, u_b, INSERT_VALUES, SCATTER_FORWARD);
- VecScatterEnd(vecscat, A_global_sol, u_b, INSERT_VALUES, SCATTER_FORWARD);
-
- ISDestroy(&localBIs);
- ISDestroy(&globalBIs);
- VecScatterDestroy(&vecscat);
-
-
-
-
- localBIndex.resize(0);
- globalBIndex.resize(0);
- localBIndex.reserve(primalDofMap[feSpace].nRankDofs * nComponents);
- globalBIndex.reserve(primalDofMap[feSpace].nRankDofs * nComponents);
-
- for (int i = 0; i < primalDofMap[feSpace].nRankDofs * nComponents; i++) {
- localBIndex.push_back(primalDofMap[feSpace].rStartDofs * nComponents + i);
- globalBIndex.push_back(localDofMap[feSpace].nOverallDofs * nComponents + primalDofMap[feSpace].rStartDofs * nComponents + i);
- }
- ISCreateGeneral(PETSC_COMM_SELF,
- localBIndex.size(),
- &(localBIndex[0]),
- PETSC_USE_POINTER,
- &localBIs);
- ISCreateGeneral(PETSC_COMM_SELF,
- globalBIndex.size(),
- &(globalBIndex[0]),
- PETSC_USE_POINTER,
- &globalBIs);
-
- VecScatterCreate(A_global_sol, globalBIs, u_primal, localBIs, &vecscat);
- VecScatterBegin(vecscat, A_global_sol, u_primal, INSERT_VALUES, SCATTER_FORWARD);
- VecScatterEnd(vecscat, A_global_sol, u_primal, INSERT_VALUES, SCATTER_FORWARD);
-
- ISDestroy(&localBIs);
- ISDestroy(&globalBIs);
- VecScatterDestroy(&vecscat);
-
- recoverSolution(u_b, u_primal, vec);
-
- KSPDestroy(&ksp);
-
-#endif
+ // The function was removed when FETI-DP was reformulated for mixed finite
+ // elements. To reconstruct this function, check for revision number 2024.
+ ERROR_EXIT("This function was removed!\n");
}
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index 7bd61b59..5926d28e 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -83,7 +83,7 @@ namespace AMDiS {
/// Create Lagrange multiplier variables corresponding to the dual
/// variables.
- void createLagrange(vector<const FiniteElemSpace*> &feSpaces);
+ void createLagrange(const FiniteElemSpace *feSpace);
/// Creates a global index of the B variables.
void createIndexB(const FiniteElemSpace *feSpace);
@@ -144,6 +144,12 @@ namespace AMDiS {
void printStatistics();
+ /// Checks whether a given DOF in a gifen FE space is a primal DOF.
+ inline bool isPrimal(const FiniteElemSpace *feSpace, DegreeOfFreedom dof)
+ {
+ return primalDofMap[feSpace].isSet(dof);
+ }
+
protected:
/// Number of components in the PDE to be solved.
int nComponents;
@@ -151,18 +157,17 @@ namespace AMDiS {
/// Mapping from primal DOF indices to a global index of primals.
FeSpaceData<GlobalDofMap> primalDofMap;
- /// Index for each non primal DOF to the global index of
- /// B variables.
- FeSpaceData<GlobalDofMap> localDofMap;
-
- /// Mapping from dual DOF indices to a global index within the local (B)
- /// variables.
+ /// Mapping from dual DOF indices to a global index of duals.
FeSpaceData<GlobalDofMap> dualDofMap;
/// Stores to each dual DOF index the index of the first Lagrange
/// constraint that is assigned to this DOF.
FeSpaceData<GlobalDofMap> lagrangeMap;
+ /// Index for each non primal DOF to the global index of
+ /// B variables.
+ FeSpaceData<GlobalDofMap> localDofMap;
+
/// Stores to each dual boundary DOF the set of ranks in which the DOF
/// is contained in.
DofIndexToPartitions boundaryDofRanks;
--
GitLab