L = f*v*dx + g("+")*v("+")*dSs(0)
works when g is a Constant or an Expression.
I'd now like to make this internal Neumann bc proportional to the difference between the solution function on either side of the boundary facets. I've tried the following code
V = FunctionSpace(mesh, "DG", 0)
...
L = f*v*dx + jump(u)*avg(v)*dSs(0)
but get the error
All terms in form must have same rank.
Is jump() the correct function to use here? If so, what's the correct usage in this instance?
