I assume, that you mean the dirac distribution by delta(x-x0).
You can approximate delta by a field that has the value 1/V in the grid cell containing x0 and 0 elsewhere, resulting in:
integral(f(x)delta(x-x0))dV=f(x0).

For a more accurate solution you may need to implement a jump boundary condition and split your mesh on each "side" of x0.