The reservoir simulation of the complex reservoirs with anisotropic permeability,which includes faults and non-orthogonal grids, with a fully discontinuous permeability tensor in the discretization is a major challenge. Several methods have already been developed and implemented within industry standard reservoir simulators for non-orthogonal grids (e.g., Multi-Point Flux Approximation (MPFA) “O” method). However, it has been noticed that some of the numerical methods for elliptic/parabolic equations may violate the maximum principle (i.e., lead to spurious oscillations), especially when the anisotropy is particularly strong. It has been found that the oscillations are closely related to the poor approximation of the pressure gradient in the flux computation. Therefore, proposed methods must correctly approximate underlying operators, satisfy a discrete maximum principle and have coercivity properties. Furthermore, the method must be robust and efficient. This paper presents the meshless multi-point flux approximation of second order elliptic operators containing a tensor coefficient. The method is based on a pressure gradient approximation commonly used in meshless methods (or Smoothed Particle Hydrodynamics method—SPH method). The proposed discretization schemes can be written as a sum of sparse positive semidefinite matrix and perturbation matrix. We show that convergence rates are retained as for finite difference methods O(hα), 1 ≤ α < 2, where h denotes the maximum particle spacing. The results are presented, discussed and future studies are outlined.