We present an efficient solver for large-scale linear elasticity problems that uses a GPU-based, matrix-free geometric multigrid (GMG) preconditioner combined with voxel-based finite element discretizations and immersed boundary methods (IBM). By using structured voxel meshes, which eliminate complex meshing, and cut-cell-based IBM techniques, we leverage a matrix-free GMG preconditioner on a GPU to dramatically reduce memory overhead and enable simulations with billions of degrees of freedom (DOF). Our implementation is built on NanoVDB, a GPU-friendly sparse voxel data structure that enables efficient memory access and parallelization. Benchmark results show the solver’s efficiency. For industrial scenarios, our approach solves billion-scale DOF problems in minutes on a single GPU, whereas commercial state-of-the-art solvers require extensive meshing and significantly more time, even for much smaller models. The solver achieves accuracy comparable to conforming mesh solutions. By eliminating the meshing bottleneck and using an efficient CUDA implementation, our approach dramatically reduces both hardware requirements and the level of expertise needed for high-fidelity simulation. This enables the simulation of complex geometries, such as lattice structures and parts with intricate details, which were previously challenging for traditional solvers.