Commit 9ce826f5 authored by Thomas Witkowski's avatar Thomas Witkowski

Mh, seems to work.

parent 56e9d5d4
......@@ -47,7 +47,8 @@ namespace AMDiS {
void GlobalDofMap::addOffset(int offset)
{
for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin(); it != dofMap.end(); ++it)
it->second.local += offset;
it->second.global = it->second.local + offset;
// it->second.local += offset;
}
void GlobalDofMap::computeNonLocalIndices()
......@@ -60,7 +61,7 @@ namespace AMDiS {
for (DofComm::Iterator it(*sendDofs, feSpace); !it.end(); it.nextRank())
for (; !it.endDofIter(); it.nextDof())
if (dofMap.count(it.getDofIndex()) && !nonRankDofs.count(it.getDofIndex()))
stdMpi.getSendData(it.getRank()).push_back(dofMap[it.getDofIndex()].local);
stdMpi.getSendData(it.getRank()).push_back(dofMap[it.getDofIndex()].global);
stdMpi.updateSendDataSize();
......@@ -84,7 +85,7 @@ namespace AMDiS {
int i = 0;
for (; !it.endDofIter(); it.nextDof())
if (nonRankDofs.count(it.getDofIndex()))
dofMap[it.getDofIndex()].local = stdMpi.getRecvData(it.getRank())[i++];
dofMap[it.getDofIndex()].global = stdMpi.getRecvData(it.getRank())[i++];
}
}
......
......@@ -369,8 +369,10 @@ namespace AMDiS {
typedef map<DegreeOfFreedom, MultiIndex>::iterator ItType;
for (ItType it = dofMap.begin(); it != dofMap.end(); ++it) {
if (data[feSpaces[i]].isRankDof(it->first)) {
int globalMatIndex =
it->second.local - data[feSpaces[i]].rStartDofs + offset;
int globalMatIndex = it->second.local + offset;
// int globalMatIndex =
// it->second.local - data[feSpaces[i]].rStartDofs + offset;
dofToMatIndex.add(i, it->first, globalMatIndex);
}
}
......@@ -400,9 +402,9 @@ namespace AMDiS {
stdMpi.startCommunication();
{
int counter = 0;
for (DofComm::Iterator it(*recvDofs, feSpaces[i]);
!it.end(); it.nextRank()) {
int counter = 0;
for (; !it.endDofIter(); it.nextDof()) {
if (dofMap.count(it.getDofIndex())) {
DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
......@@ -465,6 +467,13 @@ namespace AMDiS {
data[*it].setDofComm(pSend, pRecv);
}
inline int getMatIndex(int ithComponent, DegreeOfFreedom d)
{
FUNCNAME("FeSpaceData::getMatIndex()");
return dofToMatIndex.get(ithComponent, d);
}
private:
MPI::Intracomm* mpiComm;
......
......@@ -181,7 +181,6 @@ namespace AMDiS {
PetscSolverFeti::PetscSolverFeti()
: PetscSolver(),
nComponents(-1),
schurPrimalSolver(0)
{
FUNCNAME("PetscSolverFeti::PetscSolverFeti()");
......@@ -408,8 +407,6 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::createMatLagrange()");
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
// === Create distributed matrix for Lagrange constraints. ===
MatCreateMPIAIJ(PETSC_COMM_WORLD,
......@@ -425,35 +422,34 @@ namespace AMDiS {
// === m == r, than the rank sets -1.0 for the corresponding ===
// === constraint. ===
map<DegreeOfFreedom, MultiIndex> &dualMap = dualDofMap[feSpace].getMap();
for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
("Should not happen!\n");
for (unsigned int k = 0; k < feSpaces.size(); k++) {
map<DegreeOfFreedom, MultiIndex> &dualMap =
dualDofMap[feSpaces[k]].getMap();
// Global index of the first Lagrange constriant for this node.
int index = lagrangeMap[feSpace][it->first].local;
// 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 =
(localDofMap[feSpace].rStartDofs +
localDofMap[feSpace][it->first].local) * nComponents + k;
for (map<DegreeOfFreedom, MultiIndex>::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.getMatIndex(k, 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) {
int colIndex = localDofMap.mapGlobal(it->first, k);
double value = (W[i] == mpiRank ? 1.0 : -1.0);
MatSetValue(mat_lagrange, rowIndex, colIndex, value,
MatSetValue(mat_lagrange, index, colIndex, value,
INSERT_VALUES);
}
index++;
}
index++;
}
}
}
......@@ -767,8 +763,6 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::recoverSolution()");
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
// === Get local part of the solution for B variables. ===
PetscScalar *localSolB;
......@@ -778,19 +772,25 @@ namespace AMDiS {
// === Create scatter to get solutions of all primal nodes that are ===
// === contained in rank's domain. ===
vector<const FiniteElemSpace*> feSpaces = getFeSpaces(&vec);
vector<PetscInt> globalIsIndex, localIsIndex;
globalIsIndex.reserve(primalDofMap.getLocalDofs());
localIsIndex.reserve(primalDofMap.getLocalDofs());
{
int counter = 0;
for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
it != primalDofMap[feSpace].getMap().end(); ++it) {
for (int i = 0; i < nComponents; i++) {
globalIsIndex.push_back(it->second.local * nComponents + i);
localIsIndex.push_back(counter++);
int cnt = 0;
for (unsigned int i = 0; i < feSpaces.size(); i++) {
map<DegreeOfFreedom, MultiIndex>& dofMap =
primalDofMap[feSpaces[i]].getMap();
for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin();
it != dofMap.end(); ++it) {
globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
localIsIndex.push_back(cnt++);
}
}
TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs())("Should not happen!\n");
}
IS globalIs, localIs;
......@@ -825,21 +825,19 @@ namespace AMDiS {
// === And copy from PETSc local vectors to the DOF vectors. ===
for (int i = 0; i < nComponents; i++) {
int cnt = 0;
for (int i = 0; i < vec.getSize(); i++) {
DOFVector<double>& dofVec = *(vec.getDOFVector(i));
for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpace].getMap().begin();
it != localDofMap[feSpace].getMap().end(); ++it) {
for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpaces[i]].getMap().begin();
it != localDofMap[feSpaces[i]].getMap().end(); ++it) {
int petscIndex = localDofMap.mapLocal(it->first, i);
dofVec[it->first] = localSolB[petscIndex];
}
int counter = 0;
for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
it != primalDofMap[feSpace].getMap().end(); ++it) {
dofVec[it->first] = localSolPrimal[counter * nComponents + i];
counter++;
}
for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpaces[i]].getMap().begin();
it != primalDofMap[feSpaces[i]].getMap().end(); ++it)
dofVec[it->first] = localSolPrimal[cnt++];
}
VecRestoreArray(vec_sol_b, &localSolB);
......@@ -852,13 +850,10 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::fillPetscMatrix()");
vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
nComponents = mat->getSize();
// === Create all sets and indices. ===
vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
primalDofMap.init(mpiComm, feSpaces, true, true);
dualDofMap.init(mpiComm, feSpaces, false, false);
lagrangeMap.init(mpiComm, feSpaces, true, true);
......@@ -943,6 +938,7 @@ namespace AMDiS {
// === Traverse all sequentially created matrices and add the values to ===
// === the global PETSc matrices. ===
int nComponents = mat->getSize();
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
if (!(*mat)[i][j])
......@@ -955,7 +951,7 @@ namespace AMDiS {
for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()),
cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
bool rowPrimal = isPrimal(feSpace, *cursor);
bool rowPrimal = isPrimal(feSpaces[i], *cursor);
cols.clear();
colsOther.clear();
......@@ -972,7 +968,7 @@ namespace AMDiS {
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
bool colPrimal = primalDofMap[feSpace].isSet(col(*icursor));
bool colPrimal = isPrimal(feSpaces[j], col(*icursor));
if (colPrimal) {
if (rowPrimal) {
......@@ -996,29 +992,29 @@ namespace AMDiS {
// === For preconditioner ===
if (!rowPrimal && !colPrimal) {
int rowIndex = localDofMap[feSpace][*cursor].local;
int colIndex = localDofMap[feSpace][col(*icursor)].local;
int rowIndex = localDofMap[feSpaces[i]][*cursor].local;
int colIndex = localDofMap[feSpaces[j]][col(*icursor)].local;
if (rowIndex < nLocalInterior[feSpace]) {
if (colIndex < nLocalInterior[feSpace]) {
if (rowIndex < nLocalInterior[feSpaces[i]]) {
if (colIndex < nLocalInterior[feSpaces[j]]) {
int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
colsLocal.push_back(colIndex2);
valuesLocal.push_back(value(*icursor));
} else {
int colIndex2 =
(localDofMap[feSpace][col(*icursor)].local - nLocalInterior[feSpace]) *
(localDofMap[feSpaces[j]][col(*icursor)].local - nLocalInterior[feSpaces[j]]) *
nComponents + j;
colsLocalOther.push_back(colIndex2);
valuesLocalOther.push_back(value(*icursor));
}
} else {
if (colIndex < nLocalInterior[feSpace]) {
if (colIndex < nLocalInterior[feSpaces[j]]) {
int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
colsLocalOther.push_back(colIndex2);
valuesLocalOther.push_back(value(*icursor));
} else {
int colIndex2 = (localDofMap[feSpace][col(*icursor)].local - nLocalInterior[feSpace]) *
int colIndex2 = (localDofMap[feSpaces[j]][col(*icursor)].local - nLocalInterior[feSpaces[j]]) *
nComponents + j;
colsLocal.push_back(colIndex2);
......@@ -1034,11 +1030,9 @@ namespace AMDiS {
// === Set matrix values. ===
if (rowPrimal) {
int rowIndex =
primalDofMap[feSpace][*cursor].local * nComponents + i;
int rowIndex = primalDofMap.getMatIndex(i, *cursor);
for (unsigned int k = 0; k < cols.size(); k++)
cols[k] = primalDofMap[feSpace][cols[k]].local * nComponents + j;
cols[k] = primalDofMap.getMatIndex(j, cols[k]);
MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES);
......@@ -1061,10 +1055,8 @@ namespace AMDiS {
if (colsOther.size()) {
int globalRowIndex = localDofMap.mapGlobal(*cursor, i);
for (unsigned int k = 0; k < colsOther.size(); k++)
colsOther[k] =
primalDofMap[feSpace][colsOther[k]].local * nComponents + j;
colsOther[k] = primalDofMap.getMatIndex(j, colsOther[k]);
MatSetValues(mat_b_primal, 1, &globalRowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
......@@ -1077,10 +1069,10 @@ namespace AMDiS {
switch (fetiPreconditioner) {
case FETI_DIRICHLET:
{
int rowIndex = localDofMap[feSpace][*cursor].local;
int rowIndex = localDofMap[feSpaces[i]][*cursor].local;
if (rowIndex < nLocalInterior[feSpace]) {
int rowIndex2 = localDofMap[feSpace][*cursor].local * nComponents + i;
if (rowIndex < nLocalInterior[feSpaces[i]]) {
int rowIndex2 = localDofMap[feSpaces[i]][*cursor].local * nComponents + i;
MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
......@@ -1090,7 +1082,7 @@ namespace AMDiS {
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
} else {
int rowIndex2 =
(localDofMap[feSpace][*cursor].local - nLocalInterior[feSpace]) * nComponents + i;
(localDofMap[feSpaces[i]][*cursor].local - nLocalInterior[feSpaces[i]]) * nComponents + i;
MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
......@@ -1105,11 +1097,11 @@ namespace AMDiS {
case FETI_LUMPED:
{
int rowIndex = localDofMap[feSpace][*cursor].local;
int rowIndex = localDofMap[feSpaces[i]][*cursor].local;
if (rowIndex >= nLocalInterior[feSpace]) {
if (rowIndex >= nLocalInterior[feSpaces[i]]) {
int rowIndex2 =
(localDofMap[feSpace][*cursor].local - nLocalInterior[feSpace]) * nComponents + i;
(localDofMap[feSpaces[i]][*cursor].local - nLocalInterior[feSpaces[i]]) * nComponents + i;
MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
......@@ -1187,8 +1179,6 @@ namespace AMDiS {
FUNCNAME("PetscSolverFeti::fillPetscRhs()");
vector<const FiniteElemSpace*> feSpaces = getFeSpaces(vec);
const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
int nComponents = vec->getSize();
VecCreateMPI(PETSC_COMM_WORLD,
localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &f_b);
......@@ -1196,14 +1186,13 @@ namespace AMDiS {
primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
&f_primal);
for (int i = 0; i < nComponents; i++) {
for (unsigned int i = 0; i < feSpaces.size(); i++) {
DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
int index = dofIt.getDOFIndex();
if (isPrimal(feSpace, index)) {
index = primalDofMap[feSpace][index].local * nComponents + i;
double value = *dofIt;
VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
if (isPrimal(feSpaces[i], index)) {
index = primalDofMap.getMatIndex(i, index);
VecSetValue(f_primal, index, *dofIt, ADD_VALUES);
} else {
index = localDofMap.mapGlobal(index, i);
VecSetValue(f_b, index, *dofIt, INSERT_VALUES);
......
......@@ -151,9 +151,6 @@ namespace AMDiS {
}
protected:
/// Number of components in the PDE to be solved.
int nComponents;
/// Mapping from primal DOF indices to a global index of primals.
FeSpaceData<GlobalDofMap> primalDofMap;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment