We have developed a first-order stable Cartesian grid discretization that uses only interior grid points for inhomogeneous anisotropic elliptic operators subject to Neumann boundary conditions on a bounded nonrectangular geometry in three dimensions. For this discretization method, a necessary and sufficient condition depending on the mesh size h for the existence of this first-order stable scheme at a regular (i.e., interior) grid point is found in terms of the anisotropy matrix. For this discretization method, a way to analyze the existence of a first-order stable scheme at an irregular (i.e., boundary) grid point is also given. The arguments are identical to those for the two-dimensional case [M. A. Dumett and J. P. Keener, A Numerical Method for Solving Anisotropic Elliptic Boundary Value Problems on an Irregular Domain in 2D, manuscript]; only the details change. Unlike in [M. A. Dumett and J, P. Keener, A Numerical Method for Solving Anisotropic Elliptic Boundary Value Problems on an Irregular Domain in 2D, manuscript], a discussion of Dirichlet and Robin boundary conditions is also included. In particular, it is shown that the Gerschgorin condition does not impose sign restrictions on irregular grid points stencil coefficients as in the Neumann case.