A new numerical approach is proposed for the simulation of coupled three-dimensional and one-dimensional elliptic equations (3D-1D coupling) arising from dimensionality reduction of 3D-3D problems with thin inclusions. The method is based on a well posed mathematical formulation and results in a numerical scheme with high robustness and flexibility in handling geometrical complexities. This is achieved by means of a three-field approach to split the 1D problems from the bulk 3D problem, and then resorting to the minimization of a properly designed functional to impose matching conditions at the interfaces. Thanks to the structure of the functional, the method allows the use of independent meshes for the various subdomains.