This work is intended to provide a convenient tool for the mathematical analysis of a particular kind of finite volume approximation which can be used, for instance, in the context of nonlinear and/or anisotropic diffusion operators in 3D. Following the approach developed by F. Hermeline and by K.~Domelevo and P. Omnès in 2D, we consider a ``double'' covering $\Tau$ of a three-dimensional domain by a rather general primal mesh and by a well-chosen ``dual'' mesh. The associated discrete divergence operator $\div^{\ptTau}$ is obtained by the standard finite volume approach. A simple and consistent discrete gradient operator $\grad^\ptTau$ is defined by local affine interpolation that takes into account the geometry of the double mesh. Under mild geometrical constraints on the choice of the dual volumes, we show that $-\div^{\ptTau}$, $\grad^\ptTau$ are linked by the ``discrete duality property'', which is an analogue of the integration-by-parts formula. The primal mesh need not be conformal, and its interfaces can be general polygons. We give several numerical examples for anisotropic linear diffusion problems; good convergence properties are observed. The sequel [3] of this paper will summarize some key discrete functional analysis tools for DDFV schemes and give applications to proving convergence of DDFV schemes for several nonlinear degenerate parabolic PDEs.Read less <