  The objective of this article is to address the discretisation of
fractured/faulted poromechanical models using 3D polyhedral meshes in order to
cope with the geometrical complexity of faulted geological models.
  A polytopal scheme is proposed for contact-mechanics, based on a mixed
formulation combining a fully discrete space and suitable reconstruction
operators for the displacement field with a face-wise constant approximation of
the Lagrange multiplier accounting for the surface tractions along the
fracture/fault network. To ensure the inf--sup stability of the mixed
formulation, a bubble-like degree of freedom is included in the discrete space
of displacements (and taken into account in the reconstruction operators). It
is proved that this fully discrete scheme for the displacement is equivalent to
a low-order Virtual Element scheme, with a bubble enrichment of the VEM space.
This $\mathbb{P}^1$-bubble VEM--$\mathbb{P}^0$ mixed discretization is combined
with an Hybrid Finite Volume scheme for the Darcy flow. All together, the
proposed approach is adapted to complex geometry accounting for network of
planar faults/fractures including corners, tips and intersections; it leads to
efficient semi-smooth Newton solvers for the contact-mechanics and preserve the
dissipative properties of the fully coupled model. Our approach is investigated
in terms of convergence and robustness on several 2D and 3D test cases using
either analytical or numerical reference solutions both for the stand alone
static contact mechanical model and the fully coupled poromechanical model.
