A priori error estimates - finite element method - mixed boundary conditions

Consider the problem
$$ left{begin{array} {rcl}
-Delta u & = 0 & text{ in } Omega
u & = 0 & text{ on } Gamma_D
frac{partial u}{partial n} &= g &text{ on } Gamma_R
frac{partial u}{partial n} &= 0 &text{ on } Gamma_N
end{array}right. $$
Where $g$ is a trace of a $H^{1+varepsilon}(Omega)$ function. An example of such a configuration in is shown in the Figure below.

I am aware that depending on the function $g$ the solution to this problem might be less regular than $H^2(Omega)$. I am interested in reading about the a priori error estimates that can be obtained when solving this problem with triangular piecewise linear finite elements (P1, Lagrange).

My questions are:

• Is it possible to indicate a reference which treats in detail this classical case?
• What is the best convergence result that can be expected in this case in terms of the regularity properties of $u$? (i.e. what is the best exponent one can get in the estimate $|u-u_h|_{H^1(Omega)}leq Ch^gamma$ when $u$ is not necessarily in $H^2$)