Neumann boundary condition in 2D

[Arijit-Majumdar]

Hi,
I am trying to solve a 2D heat poisson equation in square domain with holes. This is how the domain looks like

Here, the Neumann BC are implemented at the outer edges ($\Gamma_{h1}$ to $\Gamma_{h4}$) and Dirichilet BC at the holes ($\Gamma_{g1}$ and $\Gamma_{g4}$). Thus the differential equation with BC are as follows
$-k\nabla^2 u = f$
$u = g1$ at $\Gamma_{g1}$
$u = g2$ at $\Gamma_{g2}$
$k\nabla u \cdot \hat{n} = hi$ at $\Gamma_{hi}$, for $i = 1,2,3,4$

The variational equation for the above set of equations are
$\int_{\Omega}k\nabla u \cdot \nabla v dx = \int_{\Omega}fvdx + \int_{\Gamma_{h1}}h1vdS + ... + \int_{\Gamma_{h4}}h4vdS $.

I tried implementing this in Firedrake, but it seems that the Neumann BC are not being implemented correctly. I changed the values of h1 to h4, but the final result seems unchanged which leads me to believe that the Neumann fluxes are being zero for some reason. I am attaching my code
https://colab.research.google.com/drive/1SOaM9EuobCSF15r2t3nEyYb6tPQgYW_u?usp=sharing

Here is how my linear and bilinear functional look like
a = 0.0002*dot(grad(v),grad(u))*dx
L = Constant(fe)*v*dx + Constant(h1)*v('+')*dS(1) + Constant(h2)*v('-')*dS(2) + Constant(h3)*v('+')*dS(3) + Constant(h4)*v('-')*dS(4)

where u and v are the trial and test functions respectively. Any suggestions on how to implement Neumann BC in 2D? Thanks for the help.

Arijit

[colinjcotter]

Your variational form that you wrote is the the correct one. Most likely you are not getting the boundary tags correctly created or read in. Easiest way to test this is to create a Dirichlet boundary condition of 1 on each boundary separately, then create a zero field and use bc.apply to set the boundary condition, and visualise the field. Then you can see if the tagging is correct.

[Arijit-Majumdar]

Thanks for the suggestion. As a check, I replaced one of the Neumann BC (at $\Gamma_{h1}$) with a Dirichilet BC ($u = g1$) and the solution changed. So the boundary tags in the mesh seems to be correct.
Here is the solution with all the edges of the square as Neumann

And here it is with one of Neumann BC replaced with a Dirichilet BC

I am wondering if there is something more that I need to specify in my script for Neumann BC over the edges?

[danshapero]

You can also do

fig, axes = plt.subplots()
firedrake.triplot(mesh, axes=axes)
axes.legend()

and the legend will show the colors of the different boundary markers.

[Arijit-Majumdar]

I tried the above lines and here is how the mesh looks like

The boundary labels look fine to me. I am not sure why the Neumann fluxes are then coming out to be zero?

[colinjcotter]

many of your h values *are * zero

[ksagiyam]

I think you want to use ds instead of dS: the former represents the exterior facet integrations and the latter the interior.

[Arijit-Majumdar]

Thanks! It worked.