Poisson equation with a derivative jump

General
1

Hello,
Consider a 2d domain given by a square with a disk in it. Both are meshed. Let me denote by \Omega_\text{O} the disk and \Omega_\text{S} the region between the circle and the square boundary (\partial \Omega_\text{S}).
I have the Poisson equation
\nabla^2 u = f_\text{O} in \Omega_\text{O}
\nabla^2 u = f_\text{S} in \Omega_\text{S}

with boundary conditions (BCs)
u = g on \partial \Omega_\text{S}
and
n_{\text{O} i} \partial_i u + n_{\text{S} i} \partial_i u = g on \partial \Omega_\text{S}, where n_\text{O} and n_\text{S} are the unit normals pointing outwards \Omega_\text{O} and \Omega_\text{S}, respectively.

Note that here u is continuous, but \nabla u can be discontinuous

How can I solve this in Fenics?

I tried by defining a case where I know the exact solution

u_\text{O} = 1+ x^2 + 2y^2
u_\text{S} = 1+ x^2 + 2y^2 + [(x-c_x)^2 + (y-c_y)^2 - r^2]

Without going into the details for the code (ask if needed), I set laplacian_S and laplacian_O to their respective profiles from the expressions above, and similarly for grad_S and grad_O. Scalars and vectors are defined in the function spaces

Q = FunctionSpace(mesh, 'P', 4)
V = VectorFunctionSpace(mesh, 'P', 4)

define dx_O and dx_S and the surface measures, and dS and the interior boundary (the circle), and set

F = (u.dx(i) * nu.dx(i) + laplacian_O * nu) * dx_O \
    + (u.dx(i) * nu.dx(i) + laplacian_S * nu) * dx_S\
    - facet_normal('-')[i] * grad_S[i] * nu('-') * dS\
    - facet_normal('+')[i] * grad_O[i] * nu('+') * dS

with Dirichlet BCs on \partial \Omega_{\text{S}}.

I implemented the BC on the gradient as a natural BC. However, the solution disagrees with the exact one. Do you have any clue on why?

Thank you