Commit 39e4e688 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

SHIT, SHIT, SHIT! JUST FOR JADE....

parent 66c4eef8
...@@ -221,10 +221,9 @@ namespace AMDiS { ...@@ -221,10 +221,9 @@ namespace AMDiS {
createPrimals(feSpace); createPrimals(feSpace);
createDuals(feSpace); createDuals(feSpace);
createLagrange(feSpace);
createIndexB(feSpace); createIndexB(feSpace);
} }
createLagrange(meshDistributor->getFeSpaces());
} }
...@@ -324,8 +323,7 @@ namespace AMDiS { ...@@ -324,8 +323,7 @@ namespace AMDiS {
for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
!it.end(); it.nextRank()) !it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof()) { for (; !it.endDofIter(); it.nextDof()) {
// If DOF is not primal, i.e., its a dual node if (!isPrimal(feSpace, it.getDofIndex())) {
if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
boundaryDofRanks[it.getDofIndex()].insert(mpiRank); boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
boundaryDofRanks[it.getDofIndex()].insert(it.getRank()); boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
} }
...@@ -340,7 +338,7 @@ namespace AMDiS { ...@@ -340,7 +338,7 @@ namespace AMDiS {
for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
!it.end(); it.nextRank()) !it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof()) 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.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
stdMpi.updateSendDataSize(); stdMpi.updateSendDataSize();
...@@ -349,7 +347,7 @@ namespace AMDiS { ...@@ -349,7 +347,7 @@ namespace AMDiS {
!it.end(); it.nextRank()) { !it.end(); it.nextRank()) {
bool recvFromRank = false; bool recvFromRank = false;
for (; !it.endDofIter(); it.nextDof()) { for (; !it.endDofIter(); it.nextDof()) {
if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) { if (!isPrimal(feSpace, it.getDofIndex())) {
recvFromRank = true; recvFromRank = true;
break; break;
} }
...@@ -365,7 +363,7 @@ namespace AMDiS { ...@@ -365,7 +363,7 @@ namespace AMDiS {
!it.end(); it.nextRank()) { !it.end(); it.nextRank()) {
int i = 0; int i = 0;
for (; !it.endDofIter(); it.nextDof()) for (; !it.endDofIter(); it.nextDof())
if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) if (!isPrimal(feSpace, it.getDofIndex()))
boundaryDofRanks[it.getDofIndex()] = boundaryDofRanks[it.getDofIndex()] =
stdMpi.getRecvData(it.getRank())[i++]; stdMpi.getRecvData(it.getRank())[i++];
} }
...@@ -377,25 +375,24 @@ namespace AMDiS { ...@@ -377,25 +375,24 @@ namespace AMDiS {
for (DofContainer::iterator it = allBoundaryDofs.begin(); for (DofContainer::iterator it = allBoundaryDofs.begin();
it != allBoundaryDofs.end(); ++it) it != allBoundaryDofs.end(); ++it)
if (primalDofMap[feSpace].isSet(**it) == false) if (!isPrimal(feSpace, **it))
dualDofMap[feSpace].insertRankDof(**it); dualDofMap[feSpace].insertRankDof(**it);
dualDofMap[feSpace].update(false); dualDofMap[feSpace].update(false);
MSG("nRankDuals = %d nOverallDuals = %d\n", 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()"); FUNCNAME("PetscSolverFeti::createLagrange()");
// === Reserve for each dual node, on the rank that owns this node, the === // === Reserve for each dual node, on the rank that owns this node, the ===
// === appropriate number of Lagrange constraints. === // === appropriate number of Lagrange constraints. ===
#if 0
int nRankLagrange = 0; int nRankLagrange = 0;
DofMapping& dualMap = dualDofMap[feSpace].getMap(); DofMapping& dualMap = dualDofMap[feSpace].getMap();
for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) { for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
...@@ -406,10 +403,11 @@ namespace AMDiS { ...@@ -406,10 +403,11 @@ namespace AMDiS {
} }
} }
lagrangeMap[feSpace].nRankDofs = nRankLagrange; lagrangeMap[feSpace].nRankDofs = nRankLagrange;
lagrangeMap[feSpace].update(false); lagrangeMap[feSpace].update();
MSG("nRankLagrange = %d nOverallLagrange = %d\n", MSG("nRankLagrange = %d nOverallLagrange = %d\n",
lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs); lagrangeMap[feSpace].nRankDofs,
lagrangeMap[feSpace].nOverallDofs);
// === Communicate lagrangeMap to all other ranks. === // === Communicate lagrangeMap to all other ranks. ===
...@@ -419,7 +417,7 @@ namespace AMDiS { ...@@ -419,7 +417,7 @@ namespace AMDiS {
for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
!it.end(); it.nextRank()) { !it.end(); it.nextRank()) {
for (; !it.endDofIter(); it.nextDof()) for (; !it.endDofIter(); it.nextDof())
if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) { if (!isPrimal(feSpace, it.getDofIndex())) {
DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()]; DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
stdMpi.getSendData(it.getRank()).push_back(d); stdMpi.getSendData(it.getRank()).push_back(d);
} }
...@@ -431,7 +429,7 @@ namespace AMDiS { ...@@ -431,7 +429,7 @@ namespace AMDiS {
!it.end(); it.nextRank()) { !it.end(); it.nextRank()) {
bool recvData = false; bool recvData = false;
for (; !it.endDofIter(); it.nextDof()) for (; !it.endDofIter(); it.nextDof())
if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) { if (!isPrimal(feSpace, it.getDofIndex())) {
recvData = true; recvData = true;
break; break;
} }
...@@ -446,13 +444,12 @@ namespace AMDiS { ...@@ -446,13 +444,12 @@ namespace AMDiS {
!it.end(); it.nextRank()) { !it.end(); it.nextRank()) {
int counter = 0; int counter = 0;
for (; !it.endDofIter(); it.nextDof()) { for (; !it.endDofIter(); it.nextDof()) {
if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) { if (!isPrimal(feSpace, it.getDofIndex())) {
DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++]; DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
lagrangeMap[feSpace].insert(it.getDofIndex(), d); lagrangeMap[feSpace].insert(it.getDofIndex(), d);
} }
} }
} }
#endif
} }
...@@ -469,7 +466,7 @@ namespace AMDiS { ...@@ -469,7 +466,7 @@ namespace AMDiS {
nLocalInterior = 0; nLocalInterior = 0;
for (int i = 0; i < admin->getUsedSize(); i++) { for (int i = 0; i < admin->getUsedSize(); i++) {
if (admin->isDofFree(i) == false && if (admin->isDofFree(i) == false &&
primalDofMap[feSpace].isSet(i) == false && isPrimal(feSpace, i) == false &&
dualDofMap[feSpace].isSet(i) == false) { dualDofMap[feSpace].isSet(i) == false) {
localDofMap[feSpace].insertRankDof(i); localDofMap[feSpace].insertRankDof(i);
nLocalInterior++; nLocalInterior++;
...@@ -492,10 +489,8 @@ namespace AMDiS { ...@@ -492,10 +489,8 @@ namespace AMDiS {
localDofMap[feSpace].update(false); localDofMap[feSpace].update(false);
#if 0
dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs + dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs +
nLocalInterior); nLocalInterior);
#endif
} }
...@@ -503,13 +498,13 @@ namespace AMDiS { ...@@ -503,13 +498,13 @@ namespace AMDiS {
{ {
FUNCNAME("PetscSolverFeti::createMatLagrange()"); FUNCNAME("PetscSolverFeti::createMatLagrange()");
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
// === Create distributed matrix for Lagrange constraints. === // === Create distributed matrix for Lagrange constraints. ===
MatCreateMPIAIJ(PETSC_COMM_WORLD, MatCreateMPIAIJ(PETSC_COMM_WORLD,
lagrangeMap.getRankDofs(feSpaces), lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
localDofMap.getRankDofs(), lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),
lagrangeMap.getOverallDofs(feSpaces),
localDofMap.getOverallDofs(),
2, PETSC_NULL, 2, PETSC_NULL, 2, PETSC_NULL, 2, PETSC_NULL,
&mat_lagrange); &mat_lagrange);
...@@ -520,38 +515,33 @@ namespace AMDiS { ...@@ -520,38 +515,33 @@ namespace AMDiS {
// === m == r, than the rank sets -1.0 for the corresponding === // === m == r, than the rank sets -1.0 for the corresponding ===
// === constraint. === // === constraint. ===
for (unsigned int feIdx = 0; feIdx < feSpaces.size(); feIdx++) { DofMapping &dualMap = dualDofMap[feSpace].getMap();
const FiniteElemSpace *feSpace = feSpaces[feIdx]; for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
DofMapping &dualMap = dualDofMap[feSpace].getMap(); ("Should not happen!\n");
for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
TEST_EXIT_DBG(boundaryDofRanks.count(it->first)) // Global index of the first Lagrange constriant for this node.
("Should not happen!\n"); int index = lagrangeMap[feSpace][it->first];
// Copy set of all ranks that contain this dual node.
// Global index of the first Lagrange constriant for this node. vector<int> W(boundaryDofRanks[it->first].begin(),
int constraintMatIndex = lagrangeMap.mapGlobal(it->first, feIdx); boundaryDofRanks[it->first].end());
// Copy set of all ranks that contain this dual node. // Number of ranks that contain this dual node.
vector<int> W(boundaryDofRanks[it->first].begin(), int degree = W.size();
boundaryDofRanks[it->first].end());
// Number of ranks that contain this dual node. for (int i = 0; i < degree; i++) {
int degree = W.size(); for (int j = i + 1; j < degree; j++) {
if (W[i] == mpiRank || W[j] == mpiRank) {
for (int i = 0; i < degree; i++) { // Set the constraint for all components of the system.
for (int j = i + 1; j < degree; j++) { for (int k = 0; k < nComponents; k++) {
if (W[i] == mpiRank || W[j] == mpiRank) { int rowIndex = index * nComponents + k;
#if 1
// Get global index of the dual variable
int colIndex = localDofMap.mapGlobal(it->first, feIdx);
#else
int colIndex = it->second * nComponents + k; int colIndex = it->second * nComponents + k;
#endif
double value = (W[i] == mpiRank ? 1.0 : -1.0); double value = (W[i] == mpiRank ? 1.0 : -1.0);
MatSetValue(mat_lagrange, constraintMatIndex, colIndex, value, MatSetValue(mat_lagrange, rowIndex, colIndex, value,
INSERT_VALUES); INSERT_VALUES);
} }
constraintMatIndex++;
} }
index++;
} }
} }
} }
...@@ -579,12 +569,12 @@ namespace AMDiS { ...@@ -579,12 +569,12 @@ namespace AMDiS {
localDofMap.getRankDofs(), localDofMap.getOverallDofs(), localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
&(schurPrimalData.tmp_vec_b)); &(schurPrimalData.tmp_vec_b));
VecCreateMPI(PETSC_COMM_WORLD, VecCreateMPI(PETSC_COMM_WORLD,
primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
&(schurPrimalData.tmp_vec_primal)); &(schurPrimalData.tmp_vec_primal));
MatCreateShell(PETSC_COMM_WORLD, MatCreateShell(PETSC_COMM_WORLD,
primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap.getRankDofs(), primalDofMap.getRankDofs(),
primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
&schurPrimalData, &schurPrimalData,
&mat_schur_primal); &mat_schur_primal);
MatShellSetOperation(mat_schur_primal, MATOP_MULT, MatShellSetOperation(mat_schur_primal, MATOP_MULT,
...@@ -601,8 +591,8 @@ namespace AMDiS { ...@@ -601,8 +591,8 @@ namespace AMDiS {
double wtime = MPI::Wtime(); double wtime = MPI::Wtime();
int nRowsRankPrimal = primalDofMap[feSpace].nRankDofs * nComponents; int nRowsRankPrimal = primalDofMap.getRankDofs();
int nRowsOverallPrimal = primalDofMap[feSpace].nOverallDofs * nComponents; int nRowsOverallPrimal = primalDofMap.getOverallDofs();
int nRowsRankB = localDofMap.getRankDofs(); int nRowsRankB = localDofMap.getRankDofs();
int nRowsOverallB = localDofMap.getOverallDofs(); int nRowsOverallB = localDofMap.getOverallDofs();
...@@ -634,7 +624,8 @@ namespace AMDiS { ...@@ -634,7 +624,8 @@ namespace AMDiS {
mpi::globalMax(maxLocalPrimal); mpi::globalMax(maxLocalPrimal);
TEST_EXIT(mat_b_primal_cols.size() == 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(); for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
it != mat_b_primal_cols.end(); ++it) { it != mat_b_primal_cols.end(); ++it) {
Vec tmpVec; Vec tmpVec;
...@@ -711,8 +702,6 @@ namespace AMDiS { ...@@ -711,8 +702,6 @@ namespace AMDiS {
{ {
FUNCNAME("PetscSolverFeti::createFetiKsp()"); FUNCNAME("PetscSolverFeti::createFetiKsp()");
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
// === Create FETI-DP solver object. === // === Create FETI-DP solver object. ===
fetiData.mat_primal_b = &mat_primal_b; fetiData.mat_primal_b = &mat_primal_b;
...@@ -725,19 +714,15 @@ namespace AMDiS { ...@@ -725,19 +714,15 @@ namespace AMDiS {
localDofMap.getRankDofs(), localDofMap.getOverallDofs(), localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
&(fetiData.tmp_vec_b)); &(fetiData.tmp_vec_b));
VecCreateMPI(PETSC_COMM_WORLD, VecCreateMPI(PETSC_COMM_WORLD,
lagrangeMap[feSpace].nRankDofs * nComponents, lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
lagrangeMap[feSpace].nOverallDofs * nComponents,
&(fetiData.tmp_vec_lagrange)); &(fetiData.tmp_vec_lagrange));
VecCreateMPI(PETSC_COMM_WORLD, VecCreateMPI(PETSC_COMM_WORLD,
primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
primalDofMap[feSpace].nOverallDofs * nComponents,
&(fetiData.tmp_vec_primal)); &(fetiData.tmp_vec_primal));
MatCreateShell(PETSC_COMM_WORLD, MatCreateShell(PETSC_COMM_WORLD,
lagrangeMap[feSpace].nRankDofs * nComponents, lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
lagrangeMap[feSpace].nRankDofs * nComponents, lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
lagrangeMap[feSpace].nOverallDofs * nComponents,
lagrangeMap[feSpace].nOverallDofs * nComponents,
&fetiData, &mat_feti); &fetiData, &mat_feti);
MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti); MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);
...@@ -990,9 +975,9 @@ namespace AMDiS { ...@@ -990,9 +975,9 @@ namespace AMDiS {
int nRowsRankB = localDofMap.getRankDofs(); int nRowsRankB = localDofMap.getRankDofs();
int nRowsOverallB = localDofMap.getOverallDofs(); int nRowsOverallB = localDofMap.getOverallDofs();
int nRowsRankPrimal = primalDofMap.getRankDofs(feSpaces); int nRowsRankPrimal = primalDofMap.getRankDofs();
int nRowsOverallPrimal = primalDofMap.getOverallDofs(feSpaces); int nRowsOverallPrimal = primalDofMap.getOverallDofs();
int nRowsDual = dualDofMap.getRankDofs(feSpaces); int nRowsDual = dualDofMap.getRankDofs();
int nRowsInterior = nRowsRankB - nRowsDual; int nRowsInterior = nRowsRankB - nRowsDual;
MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL, MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
...@@ -1341,14 +1326,14 @@ namespace AMDiS { ...@@ -1341,14 +1326,14 @@ namespace AMDiS {
VecCreateMPI(PETSC_COMM_WORLD, VecCreateMPI(PETSC_COMM_WORLD,
localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &f_b); localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &f_b);
VecCreateMPI(PETSC_COMM_WORLD, VecCreateMPI(PETSC_COMM_WORLD,
primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
primalDofMap[feSpace].nOverallDofs * nComponents, &f_primal); &f_primal);
for (int i = 0; i < nComponents; i++) { for (int i = 0; i < nComponents; i++) {
DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS); DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) { for (dofIt.reset(); !dofIt.end(); ++dofIt) {
int index = dofIt.getDOFIndex(); int index = dofIt.getDOFIndex();
if (primalDofMap[feSpace].isSet(index)) { if (isPrimal(feSpace, index)) {
index = primalDofMap[feSpace][index] * nComponents + i; index = primalDofMap[feSpace][index] * nComponents + i;
double value = *dofIt; double value = *dofIt;
VecSetValues(f_primal, 1, &index, &value, ADD_VALUES); VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
...@@ -1408,203 +1393,9 @@ namespace AMDiS { ...@@ -1408,203 +1393,9 @@ namespace AMDiS {
{ {
FUNCNAME("PetscSolverFeti::solveFetiMatrix()"); FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
#if 0 // The function was removed when FETI-DP was reformulated for mixed finite
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0); // elements. To reconstruct this function, check for revision number 2024.
ERROR_EXIT("This function was removed!\n");
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,