  In the present work we introduce a novel refinement algorithm for
two-dimensional elliptic partial differential equations discretized with
Virtual Element Method (VEM). The algorithm improves the numerical solution
accuracy and the mesh quality through a controlled refinement strategy applied
to the generic polygonal elements of the domain tessellation. The numerical
results show that the outlined strategy proves to be versatile and applicable
to any two-dimensional problem where polygonal meshes offer advantages. In
particular, we focus on the simulation of flow in fractured media, specifically
using the Discrete Fracture Network (DFN) model. A residual a-posteriori error
estimator tailored for the DFN case is employed. We chose this particular
application to emphasize the effectiveness of the algorithm in handling complex
geometries. All the numerical tests demonstrate optimal convergence rates for
all the tested VEM orders.
