solve a PDE with an integral term

Question

I want to use FEniCS to solve a PDE with an integral of the unknown in either the PDE or in a boundary condition. A simple example of this would be something of the form

$-\Delta u + \int_\Omega u = f$ in $\Omega$

with zero Dirichlet boundary conditions. How can the weak form of the left-hand side of the PDE be written as a bilinear form that FEniCS can understand? I naively tried

a = inner(grad(u), grad(v))*dx + assemble(u*dx)*v*dx

but this didn't work. Any help would be appreciated.

Comments

This is something that I am also currently trying to do myself, and I have not been successful. I will be very interested in knowing how to resolve this - if some one on this thread is able to help.

Answer 1

We could do an optimization with constraints, i.e., introduce the mean of u by introducing an additional globally constant Lagrange multiplier m. I did not test this or check if this is well-posed and all, but check this out and see if it helps...

from dolfin import *

mesh = UnitSquareMesh(8, 8)
f = Constant(1)

V = FunctionSpace(mesh, 'Lagrange', 1)
R = FunctionSpace(mesh, 'R', 0)
X = V*R

u, m = TrialFunctions(X)
v, r = TestFunctions(X)

F = dot(grad(u),grad(v))*dx + m*v*dx + (u-m)*r*dx - f*v*dx

x = Function(X)
bc = DirichletBC(X.sub(0), Constant(0), 'on_boundary')
solve(lhs(F) == rhs(F), x, bc)

uh, mh = x.split()
plot(uh)
plot(mh)
interactive()

Testing with (v,r) = (0,1) you see that if you obtain a solution, then m is in fact the mean of u (you may need to be careful for non-unit domains). The first equation (test with (v,0)) for any v in H^1_0 then results from deriving the usual weak form and replacing the mean of u with m.

Comments

That worked well, thanks! Also, it is not hard to change it so that it works for non-unit domains.

---

Yes, basically you just need to compute the area. This can be done by looping over the mesh and summing up element volumes or simply using assemble(1*dx(mesh)).