The present work deals with the numerical treatment of coupled sub-surface flow and geomechanical deformation in the presence of fractures. This is relevant for a large variety of practical applications, such as the assessment of the safety and of the environmental impact of human activities concerning the subsoil. Fractures are represented as 2-dimensional planar polygonal regions embedded in a 3-dimensional porous matrix. This setting results in a hybrid dimensional computational domain with a possibly discontinuous pressure solution across the 2-D interfaces. The main challenges in the above setting lie in the size and geometrical complexity of the domain. Such complexities are here addressed though the use of a domain decomposition approach on non conforming meshes. The method is based on the introduction of additional interface variables and on the minimization of a quadratic functional to achieve the coupling. The effectiveness of the methodology is discussed through the presentation of some numerical tests.