I find myself in the position of having to solve
$-\Delta u = f$ on a subset of Cartesian grid points that don't necessarily form a cuboid domain subject to a homogenious Dirichlet boundary condition ($u = 0$). The Laplacian filter I require is the standard 7-point stencil for 3D.
Now, for a cuboid domain, I know this can be efficiently solved using a type-I discrete sine transform (DST) since the DST diagonalizes the equation system.
I wonder, whether it's possible to leverage the fast DST for domains that are not cuboid. It's still a Cartesian grid, but only a subset which is not necessarily convex and might include holes.
Are there fast algorithms that can exploit the Cartesian grid structure? Or do I have to use a more general solver (like an algebraic multigrid preconditioned conjugate gradient method or something like this)?
Maybe there is an efficient method to reduce my problem to another with an emcompassing cuboid domain so I still get to use the DST? As far as I can tell, I could introduce source terms in $f$ on my actual boundary so that the DST-based approach will reconstruct a boundary of zero. But the problem with this is that determining those source terms seems to require solving a big and "ugly" linear equation system (dense matrix which is not trivially diagonalizable).
What would be your recommended solver for this problem? The number of unknowns are in the tens of millions.