This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

solve a PDE with an integral term

+2 votes

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.

asked Sep 18, 2014 by bseguin FEniCS Novice (650 points)

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.

1 Answer

+1 vote

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.

answered Sep 20, 2014 by Christian Waluga FEniCS Expert (12,310 points)

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)).

...