In contrast to the non-relativistic approaches, three-dimensional (3D) mesh calculations for the {it relativistic} density functional theory have not been realized because of the challenges of variational collapse and fermion doubling. We overcome these difficulties by developing a novel method based on the ideas of Wilson fermion as well as the variational principle for the inverse Hamiltonian. We demonstrate the applicability of this method by applying it to $^{16}$O, $^{24}$Mg, and $^{28}$Si nuclei, providing detailed explanation on the formalism and verification of numerical implementation.