Commit a9ba2f66 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* UMFPack solver can now deal with muliple componets having different number of nodes

parent 493f4b49
...@@ -44,8 +44,22 @@ namespace AMDiS { ...@@ -44,8 +44,22 @@ namespace AMDiS {
// Number of systems. // Number of systems.
int nComponents = m->getSize(); int nComponents = m->getSize();
// Size of the new composed matrix.
int newMatrixSize = ((*m)[0][0])->getFESpace()->getAdmin()->getUsedSize() * nComponents; // Calculate size of the new composed matrix.
int newMatrixSize = 0;
std::vector<int> matricesSize(nComponents);
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
if ((*m)[i][j]) {
matricesSize[i] = ((*m)[i][j])->getFESpace()->getAdmin()->getUsedSize();
newMatrixSize += matricesSize[i];
break;
}
}
}
// int newMatrixSize = ((*m)[0][0])->getFESpace()->getAdmin()->getUsedSize() * nComponents;
// The new matrix has to be stored in compressed col format, therefore // The new matrix has to be stored in compressed col format, therefore
// the cols are collected. // the cols are collected.
...@@ -54,34 +68,44 @@ namespace AMDiS { ...@@ -54,34 +68,44 @@ namespace AMDiS {
// Counter for the number of non-zero elements in the new matrix. // Counter for the number of non-zero elements in the new matrix.
int nElements = 0; int nElements = 0;
for (int stencilRow = 0; stencilRow < nComponents; stencilRow++) { for (int stencilRow = 0, shiftRow = 0; stencilRow < nComponents; stencilRow++) {
int shiftCol = 0;
for (int stencilCol = 0; stencilCol < nComponents; stencilCol++) { for (int stencilCol = 0; stencilCol < nComponents; stencilCol++) {
if (!(*m)[stencilRow][stencilCol]) { if (!(*m)[stencilRow][stencilCol]) {
shiftCol += matricesSize[stencilCol];
continue; continue;
} }
DOFMatrix::Iterator matrixRow((*m)[stencilRow][stencilCol], USED_DOFS); DOFMatrix::Iterator matrixRow((*m)[stencilRow][stencilCol], USED_DOFS);
int rowIndex = 0; int rowIndex = 0;
for (matrixRow.reset(); !matrixRow.end(); matrixRow++, rowIndex++) { for (matrixRow.reset(); !matrixRow.end(); matrixRow++, rowIndex++) {
for (int i = 0; i < static_cast<int>((*matrixRow).size()); i++) { for (int i = 0; i < static_cast<int>((*matrixRow).size()); i++) {
if ((*matrixRow)[i].col >= 0) { if ((*matrixRow)[i].col >= 0) {
MatEntry me; MatEntry me;
me.entry = (*matrixRow)[i].entry; me.entry = (*matrixRow)[i].entry;
// The col field is used to store the row number of the new element. // The col field is used to store the row number of the new element.
me.col = (rowIndex * nComponents) + stencilRow; me.col = rowIndex + shiftRow;
// And save the new element in the corresponding column. // And save the new element in the corresponding column.
cols[((*matrixRow)[i].col * nComponents) + stencilCol].push_back(me); cols[(*matrixRow)[i].col + shiftCol].push_back(me);
nElements++; nElements++;
} }
}
}
}
}
shiftCol += matricesSize[stencilCol];
} }
shiftRow += matricesSize[stencilRow];
} }
// Data fields for UMFPack. // Data fields for UMFPack.
int *Ap = (int*)malloc(sizeof(int) * (newMatrixSize + 1)); int *Ap = (int*)malloc(sizeof(int) * (newMatrixSize + 1));
int *Ai = (int*)malloc(sizeof(int) * nElements); int *Ai = (int*)malloc(sizeof(int) * nElements);
...@@ -90,15 +114,15 @@ namespace AMDiS { ...@@ -90,15 +114,15 @@ namespace AMDiS {
double *xvec = (double*)malloc(sizeof(double) * newMatrixSize); double *xvec = (double*)malloc(sizeof(double) * newMatrixSize);
// Resort the right hand side of the linear system. // Resort the right hand side of the linear system.
for (int i = 0; i < b->getSize(); i++) { for (int i = 0, counter = 0; i < b->getSize(); i++) {
DOFVector<double>::Iterator it(b->getDOFVector(i), USED_DOFS); DOFVector<double>::Iterator it(b->getDOFVector(i), USED_DOFS);
int counter = 0; for (it.reset(); !it.end(); ++it) {
for (it.reset(); !it.end(); ++it, counter++) { bvec[counter++] = *it;
bvec[counter * nComponents + i] = *it;
} }
} }
// Create fields Ap, Ai and Ax. // Create fields Ap, Ai and Ax.
int elCounter = 0; int elCounter = 0;
Ap[0] = 0; Ap[0] = 0;
...@@ -164,12 +188,11 @@ namespace AMDiS { ...@@ -164,12 +188,11 @@ namespace AMDiS {
// Copy and resort solution. // Copy and resort solution.
for (int i = 0; i < x->getSize(); i++) { for (int i = 0, counter = 0; i < x->getSize(); i++) {
DOFVector<double>::Iterator it(x->getDOFVector(i), USED_DOFS); DOFVector<double>::Iterator it(x->getDOFVector(i), USED_DOFS);
int counter = 0; for (it.reset(); !it.end(); it++) {
for (it.reset(); !it.end(); it++, counter++) { *it = xvec[counter++];
*it = xvec[counter * nComponents + i];
} }
} }
......
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