x
inline int LinearIndex(const int i, const int j, const int k) { return ((i * YDIM) + j) * ZDIM + k; }CSRMatrix BuildLaplacianMatrix() { static constexpr int matSize = XDIM * YDIM * ZDIM; CSRMatrixHelper matrixHelper(matSize); for (int i = 1; i < XDIM - 1; i++) for (int j = 1; j < YDIM - 1; j++) for (int k = 1; k < ZDIM - 1; k++) { matrixHelper(LinearIndex(i, j, k), LinearIndex(i, j, k)) = -6.; matrixHelper(LinearIndex(i, j, k), LinearIndex(i + 1, j, k)) = 1.; matrixHelper(LinearIndex(i, j, k), LinearIndex(i - 1, j, k)) = 1.; matrixHelper(LinearIndex(i, j, k), LinearIndex(i, j + 1, k)) = 1.; matrixHelper(LinearIndex(i, j, k), LinearIndex(i, j - 1, k)) = 1.; matrixHelper(LinearIndex(i, j, k), LinearIndex(i, j, k + 1)) = 1.; matrixHelper(LinearIndex(i, j, k), LinearIndex(i, j, k - 1)) = 1.; } return matrixHelper.ConvertToCSRMatrix();}CSRMatrix BuildPreconditionerMatrix() { static constexpr int matSize = XDIM * YDIM * ZDIM; CSRMatrixHelper matrixHelper(matSize); // A cheap but effective approximation of an Incomplete Cholesky preconditioner // NOTE : You don't need to explain *why* this ends up being an effective preconditioner for (int i = 1; i < XDIM - 1; i++) for (int j = 1; j < YDIM - 1; j++) for (int k = 1; k < ZDIM - 1; k++) { matrixHelper(LinearIndex(i, j, k), LinearIndex(i, j, k)) = 3.05; // Basically, 3+epsilon, where epsilon is small, // but not "too" small. epsilon = .05 works well // for grid sizes in the order of 128^3 - 256^3 // This will be a lower-triangular matrix ... if (i > 1) matrixHelper(LinearIndex(i, j, k), LinearIndex(i - 1, j, k)) = -1.; if (j > 1) matrixHelper(LinearIndex(i, j, k), LinearIndex(i, j - 1, k)) = -1.; if (k > 1) matrixHelper(LinearIndex(i, j, k), LinearIndex(i, j, k - 1)) = -1.; } // Need to put part of the "identity" matrix on nodes of the margin for (int i = 0; i < XDIM; i++) for (int j = 0; j < YDIM; j++) for (int k = 0; k < ZDIM; k++) if (i == 0 || i == XDIM - 1 || j == 0 || j == YDIM - 1 || k == 0 || k == ZDIM - 1) matrixHelper(LinearIndex(i, j, k), LinearIndex(i, j, k)) = 1.; return matrixHelper.ConvertToCSRMatrix();}void ComputeLaplacian(CSRMatrix& laplacianMatrix, const float (&u)[XDIM][YDIM][ZDIM], float (&Lu)[XDIM][YDIM][ZDIM]) { // Treat the arrays u & Lu as flattened vectors, and apply matrix-vector multiplication MatVecMultiply(laplacianMatrix, &u[0][0][0], &Lu[0][0][0]);}