I'd like to create a Neumann boundary condition that is different depending on the magnitude of the unknown variable.
Now, if u is the unknown of the system, I want something like
Q = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(Q)
v = TestFunction(Q)
u_m = interpolate(Expression("pow(x[0],2)"), Q)
g_b = conditional( lt(u, u_m), 200.0, 0.0)
F = inner(grad(v), grad(u))*dx + g_b*v*ds
solve(lhs(F) == rhs(F))
Of course, when solving this, I get the error
raise ArityMismatch("Integrand arguments {0} differ from form arguments {1}.".format(args, arguments))
Is there some other _smart_ way of doing this (besides fixed-point iteration) that anyone knows about?
Thanks!