  A new numerical scheme is proposed for flow computation in complex discrete
fracture networks. The method is based on a three-field formulation of the
Darcy law for the description of the hydraulic head on the fractures and uses a
cost functional to enforce the required coupling condition at fracture
intersections. The resulting method can handle non conforming meshes,
independently built on each geometrical object of the computational domain, and
ensures local mass conservation properties at fracture intersections. An
iterative solver is devised for the method, ready for parallel implementation
on parallel computing architectures.
