Biharmonic equation with Dirichlet and Neumann BCs

Question

I'm interested in solving the biharmonic equation
$$\Delta^{2}u=0 \quad\text{on}\quad\Omega$$
with the boundary conditions
$$u=0 \quad\text{on}\quad\partial\Omega$$
$$\frac{\partial u}{\partial n}=0 \quad\text{on}\quad\partial\Omega$$
where $n$ is the outward unit normal to the boundary. I have a finite difference solver implemented for this but am interested in obtaining an FEM solution with fenics for validation purposes and to see if fenics is something I might use for future problems. Being new to fenics, and admittedly being a bit rusty on FEM methods, I was hoping to adapt the biharmonic demo code to these boundary conditions (for now I have left the source as is in the demo). In a comment to this question here, nate references a paper by Georgoulis and Houston which goes through a derivation and analysis of an interior penalty method for the more general boundary conditions
$$u=g_{D} \quad\text{on}\quad\partial\Omega$$
$$\frac{\partial u}{\partial n}=g_{N} \quad\text{on}\quad\partial\Omega$$.
Having skimmed through this I have changed the bilinear form in the demo code to

alpha = Constant(8.0)
beta = Constant(8.0)
a = inner(div(grad(u)), div(grad(v)))*dx \
  - inner(avg(div(grad(u))), jump(grad(v), n))*dS \
  - inner(jump(grad(u), n), avg(div(grad(v))))*dS \
  + beta*inner(jump(grad(u),n), jump(grad(v),n))*dS \
  + inner(avg(dot(grad(div(grad(u))),n)), jump(v))*dS \
  + inner(jump(u), avg(dot(grad(div(grad(v))),n)))*dS \
  + alpha*inner(jump(u), jump(v))*dS 

but the Neumann boundary condition doesn't appear to be enforced correctly. Actually it seems the boundary condition from the demo
$$\Delta u=0 \quad\text{on}\quad\partial\Omega$$
is still being enforced. I feel like I am probably missing something that would be fairly obvious to an experienced FEM/fenics user. Anything pointing me in the right direction would be appreciated.

Comments

What are your exterior boundary terms? And what function space are you employing?

---

I was using the same 2nd degree 'CG' space as in the demo, but am also interested in trying 3rd degree elements when I couple this to another PDE at a later stage.

I think I've worked it out actually, I hadn't realised that dS only applied to the interior facets, which I think is what you are getting at in your first question. After going through the paper in more detail and adding the appropriate ds terms I am getting a much better result. I ended up with the following bilinear form which is passed to the solver with no explicit boundary conditions

a = inner(div(grad(u)), div(grad(v)))*dx \
  - inner(avg(div(grad(u))), jump(grad(v), n))*dS \
  - inner(jump(grad(u), n), avg(div(grad(v))))*dS \
  + beta/h_avg*inner(jump(grad(u),n), jump(grad(v),n))*dS \
  + inner(avg(dot(grad(div(grad(u))),n)), jump(v))*dS \
  + inner(jump(u), avg(dot(grad(div(grad(v))),n)))*dS \
  + alpha/(h_avg**3)*inner(jump(u), jump(v))*dS \
  - inner(div(grad(u)),dot(grad(v),n))*ds \
  - inner(dot(grad(u),n),div(grad(v)))*ds \
  + beta/h*inner(dot(grad(u),n),dot(grad(v),n))*ds \
  + inner(grad(div(grad(u))), v*n)*ds \
  + inner(u*n, grad(div(grad(v))))*ds \
  + alpha/(h**3)*inner(u*n,v*n)*ds 

I could also potentially get rid of terms relating to the $u=0$ boundary condition and pass that explicitly to solve, do you think this makes much difference?

As a side note, after going through this I feel that perhaps the biharmonic demo could perhaps be a little better documented, particularly with respect to how the boundary conditions are enforced. For example, it would appear to me that perhaps $\Delta u=0$ on the boundary is being enforced automatically/by default in the demo much like a Neumann boundary condition is automatic/default in a second order PDE, is this correct?