Thin wall boundary condition

1

Is it possible to impose a boundary condition in fenics that is called thin wall and is useful on MHD flows? The boundary condition is:

-df/dn = c d2f/dt2,

where df/dn is the wall normal derivative of a scalar electric potential f, c a constant, and d2f/dt2 the partial Laplacian of f at the wall (partial means that the Laplacian has all other components except of the one normal to the wall, d2f/dn2).
Thank you for the help!

2

If the normal derivative is the flux corresponding to a differential equation governing f on the interior (which it sounds like it is, if f is an electric potential satisfying the Poisson equation), then maybe something like the following would work (with derivation and notes in comments):

from dolfin import *
mesh = UnitSquareMesh(100,100)
V = FunctionSpace(mesh,"CG",1)
u = TrialFunction(V)
v = TestFunction(V)
uh = Function(V)

c = Constant(1.0e0)
x = SpatialCoordinate(mesh)
f = sin(2.0*pi*x[0])*cos(2.0*pi*x[1])
n = FacetNormal(mesh)

# Tangential gradient operator:
def grad_t(u):
    return grad(u) - dot(grad(u),n)*n

# Residual of Poisson problem with the specified BC:
res = inner(grad(u),grad(v))*dx - c*inner(grad_t(u),grad_t(v))*ds \
      - inner(f,v)*dx

# Derivation:
#
# 1. Integrate -div(grad(u))*v*dx by parts to get
#
#      inner(grad(u),grad(v))*dx - dot(grad(u),n)*v*ds
#
# 2. Replace -dot(grad(u),n) with c times the "partial Laplacian", which I'm
#    interpreting as the Laplace--Beltrami operator on the boundary:
#
#      inner(grad(u),grad(v))*dx + c*LaplaceBeltrami(u)*v*ds
#
# 3. Integrate by parts on the boundary itself, to get the weak form of the
#    Laplace--Beltrami operator, in terms of the tangential gradient.
#
# NOTE: Check the signs; I attempted to follow "-df/dn = c d2f/dt2" from the
# original post, but the boundary term looks like it may not be stable in
# general, based on coercivity.

# Pin down one corner to get rid of constant mode:
bc = DirichletBC(V,Constant(0.0),
                 "near(x[0],0) && near(x[1],0)","pointwise")

# Take a look at the solution:
solve(lhs(res)==rhs(res),uh,[bc,])

from matplotlib import pyplot as plt
plot(uh)
plt.show()
3

Thank you very much for the quick response!!