max/tug
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/tug.git synced 2025-12-13 09:28:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

still an error inside indexing of sparse matrices

Browse Source
This commit is contained in:
Max Lübke 2021-10-26 13:30:55 +02:00
parent 4353cdc133
commit 72f527a59b
1 changed files with 32 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
src/diffusion.cpp
View File

@ -16,29 +16,47 @@ void BTCS2D(int x, int y, std::vector<double> &c, std::vector<double> &alpha,
double dx = 1. / x;
double dy = 1. / y;
int size = x * y;
int size = (x - 2) * (y - 2);
int local_x = x - 2;
Eigen::VectorXd b = Eigen::VectorXd::Constant(size, 0);
Eigen::VectorXd x_out(size);
std::vector<T> tripletList;
tripletList.reserve(((x - 1) * (y - 1) * 5));
tripletList.reserve(size * 5);
for (int i = 1; i < y - 1; i++) {
for (int j = 1; j < x - 1; j++) {
double sx = (alpha[i * y + j] * timestep) / (dx * dx);
double sy = (alpha[i * y + j] * timestep) / (dy * dy);
for (int i = 0; i < (size-local_x); i++) {
int current = local_x + i;
double sx = (alpha[current] * timestep) / (dx * dx);
double sy = (alpha[current] * timestep) / (dy * dy);
tripletList.push_back(T((i * x) + j, (i * x) + j, (1 + 2 * sx + 2 * sy)));
tripletList.push_back(T((i * x) + j, ((i * x) + j) + x, sy));
tripletList.push_back(T((i * x) + j, ((i * x) + j) - x, sy));
tripletList.push_back(T((i * x) + j, ((i * x) + j) + 1, sx));
tripletList.push_back(T((i * x) + j, ((i * x) + j) - 1, sx));
tripletList.push_back(T(i, current, (1 + 2 * sx + 2 * sy)));
tripletList.push_back(T(i, current + local_x, sy));
tripletList.push_back(T(i, current - local_x, sy));
tripletList.push_back(T(i, current + 1, sx));
tripletList.push_back(T(i, current - 1, sx));
b[(i * x) + j] = -c[(i * x) + j];
}
std::cout << current << std::endl;
b[i] = -c[x+i];
}
Eigen::SparseMatrix<double> A(size, size);
// for (int i = 0; i < y; i++) {
// for (int j = 0; j < x; j++) {
// double sx = (alpha[i * y + j] * timestep) / (dx * dx);
// double sy = (alpha[i * y + j] * timestep) / (dy * dy);
// tripletList.push_back(T((i * x) + j, (i * x) + j, (1 + 2 * sx + 2 *
// sy))); tripletList.push_back(T((i * x) + j, ((i * x) + j) + x, sy));
// tripletList.push_back(T((i * x) + j, ((i * x) + j) - x, sy));
// tripletList.push_back(T((i * x) + j, ((i * x) + j) + 1, sx));
// tripletList.push_back(T((i * x) + j, ((i * x) + j) - 1, sx));
// b[(i * x) + j] = -c[(i * x) + j];
// }
// }
Eigen::SparseMatrix<double> A(size, (x*y)-4);
A.setFromTriplets(tripletList.begin(), tripletList.end());
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>>