In the following, I've addressed this from a numerical perspective, as we're on scicomp.SE.

For the Poisson equation with homogeneous Dirichlet boundary conditions, there is a unique solution for any right-hand side. Once discretized, the equation can be written in the form $Ax = b$, where $A$ is the standard discretization of the 3D Laplacian operator with Dirichlet boundary and $b$ is a standard discretization of $f$. Since $A$ is positive definitive, it is invertible, and the system will have a unique solution for any $b$, and thus any $f$.

There are, of course, the usual issues with undersampling; if two different values of $f$ give rise to the same discretization due to use of coarse grid or discontinuities in $f$, there may be some ambiguity about what system is actually being solved. But provided the discretization of $f$ is well-behaved, a meaningful unique solution will exist.

The situation is slightly more complicated in the case of periodic boundary conditions because the standard discretization of the 3D Laplacian operator with periodic boundary is positive semidefinitive, and has a one-dimensional kernel $K$ comprising solutions of the form $x \equiv C$ with $C$ constant. Because $A$ is still symmetric, we have that $\operatorname{Range} A = K^\perp$, and so $Ax = b$ will not have a solution unless $\mathbf{1} \cdot b = 0$, where $\mathbf{1}$ is the vector consisting of all 1s. This provides the consistency condition for the right-hand side, at least in numerical form.

From a numerical perspective, it's perhaps easiest to discuss the discretizations directly.

For the Poisson equation with homogeneous Dirichlet boundary conditions, there is a unique solution for any right-hand side. Once discretized, the equation can be written in the form $Ax = b$, where $A$ is the standard discretization of the 3D Laplacian operator with Dirichlet boundary and $b$ is a standard discretization of $f$. Since $A$ is positive definitive, it is invertible, and the system will have a unique solution for any $b$, and thus any $f$.

There are, of course, the usual issues with undersampling; if two different values of $f$ give rise to the same discretization due to use of a coarse grid or due to discontinuities in $f$, there may be some ambiguity about what system is actually being solved. But provided the discretization of $f$ is well-behaved, a meaningful unique solution will exist.

The situation is slightly more complicated in the case of periodic boundary conditions because the standard discretization of the 3D Laplacian operator with periodic boundary is positive semidefinite and has a one-dimensional kernel $K$ comprising solutions of the form $x \equiv C$ with $C$ constant. 

Because $A$ is still symmetric in the periodic case, we have that $\operatorname{Range} A = K^\perp$, and so $Ax = b$ will not have a solution unless $\mathbf{1} \cdot b = 0$, where $\mathbf{1}$ is the vector consisting of all 1s. This provides the consistency condition for the right-hand side in discretized form.

Note that, analytically, there is a somewhat simpler way to look at it. Recall that, for $\phi \in C^2(\Omega)$, where $\Omega$ is our domain, we have $$\int_\Omega \Delta \phi \, dx = \int_{\partial \Omega} \frac{\partial \phi}{\partial \mathbf{n}}\, dS.$$ If we stipulate a periodic boundary condition on $\phi$, then the boundary term in the right-hand side disappears, and we're left with $$\int_\Omega \Delta\phi\, dx = 0.$$ So if $\phi$ satisfies $\Delta \phi = f$, it immediately follows that we must have $$\int_\Omega f\, dx = 0.$$ This is the analytical analog to $\mathbf{1} \cdot b = 0$, as both are expressing the fact that the average value of $f$, and thus $b$, over the domain must be zero.