The mathematical modeling of mixed-dimensional problems, such as coupled 3D-1D and 3D-2D systems, often leads to computational domains with highly complex geometries. In practical applications, these geometric challenges usually arise from heterogeneous features such as stones, impermeable barriers, fractures, or layered soil structures. These difficulties give rise to intricate geometrical configurations that significantly complicate mesh generation and numerical discretization. To address these challenges, the Virtual Element Method (VEM) has emerged as a particularly suitable numerical framework. A key advantage of VEM is its ability to operate on highly general meshes composed of arbitrary polygonal and polyhedral elements, enabling aligned edges and faces and offering great flexibility in the treatment of complex geometries. For realistic large-scale simulations, however, the computational cost associated with highly refined meshes can become prohibitive. As a result, adaptive strategies aimed at efficiently reducing the number of degrees of freedom in the resulting linear systems are essential. While refinement techniques, corresponding to the REFINE step in the classical adaptive paradigm, are well established for simplicial or structured meshes, their extension to fully general polytopal meshes remains a challenging open problem. This presentation addresses these issues by introducing a novel refinement strategy for three-dimensional polyhedral meshes. The proposed approach is designed to preserve, and in some cases improve, the quality of the initial coarse mesh while avoiding the creation of poorly shaped elements, such as excessively small edges, faces, or polyhedra. At the same time, the method remains aware to the lower-dimensional (2D and 1D) coupled features inherited from the original mixed-dimensional problem. Numerical examples are presented to demonstrate the effectiveness and robustness of the proposed strategy.