# Poisson equation with a derivative jump

**URL:** <https://fenicsproject.discourse.group/t/poisson-equation-with-a-derivative-jump/19646>\
**Category:** General\
**Created:** [April 10, 2026, 9:09am UTC](https://fenicsproject.discourse.group/t/poisson-equation-with-a-derivative-jump/19646 "2026-04-10T09:09:56Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![mekong](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mekong/32/7729_2.png) [@mekong](https://fenicsproject.discourse.group/u/mekong)\
**Post date:** [April 10, 2026, 9:09am UTC](https://fenicsproject.discourse.group/t/poisson-equation-with-a-derivative-jump/19646/1 "2026-04-10T09:09:56Z")

</div>

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

```auto
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

```auto
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