  In order to construct regularizations of continuous linear functionals acting
on Sobolev spaces such as $W_0^{1,q}(\Omega)$, where $1<q<\infty$ and $\Omega$
is a Lipschitz domain, we propose a projection method in negative Sobolev
spaces $W^{-1,p}(\Omega)$, $p$ being the conjugate exponent satisfying $p^{-1}
+ q^{-1} = 1$. Our method is particularly useful when one is dealing with a
rough (irregular) functional that is a member of $W^{-1,p}(\Omega)$, though not
of $L^1(\Omega)$, but one strives for a regular approximation in $L^1(\Omega)$.
We focus on projections onto discrete finite element spaces $G_n$, and consider
both discontinuous as well as continuous piecewise-polynomial approximations.
  While the proposed method aims to compute the best approximation as measured
in the negative (dual) norm, for practical reasons, we will employ a
computable, discrete dual norm that supremizes over a discrete subspace $V_m$.
We show that this idea leads to a fully discrete method given by a mixed
problem on $V_m\times G_n$. We propose a discontinuous as well as a continuous
lowest-order pair, prove that they are compatible, and therefore obtain
quasi-optimally convergent methods.
  We present numerical experiments that compute finite element approximations
to Dirac delta's and line sources. We also present adaptively generate meshes,
obtained from an error representation that comes with the method. Finally, we
show how the presented projection method can be used to efficiently compute
numerical approximations to partial differential equations with rough data.
