Computing Green's function with FEniCS

Question

Using the FEniCS library in python I wish to find a numerical solution on a 2D domain to an equation such as

Lf(x) = d(x)

with Dirichlet boundary conditions, where L is a linear differential operator and d(x) is the Dirac delta function centered at (0,0).

The following code works (but doesn't solve the equation above):

from dolfin import *

mesh = Mesh("mesh_file.xml")
V = FunctionSpace(mesh, "Lagrange", 1)

class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

bc = DirichletBC(V, Constant(0), DirichletBoundary())

u = TrialFunction(V)
v = TestFunction(V)

f = Expression("-exp(x[0]*x[0] + x[1]*x[1])")
equation = inner(grad(u), grad(v))*dx == f*v*dx

w = Function(V)
solve(equation, w, bc)

---

I am having trouble implementing the right-hand term of the equation. In the FEniCS documentation, I found the PointSource class but I don't understand how to use it.

I tried replacing the value of f in the snippet above with

f = PointSource(V, Point(0,0))

but this returns TypeError: unsupported operand type(s) for *: 'PointSource' and 'Argument'

I also considered using a function f such that f(x)=1/3 if x is a node of the cell that contains the center of the delta function, and f(x)=0 elsewise, but I can't figure out how to implement this idea so that it can be used with FEniCS.

---

How can I put a Dirac delta function in the right-hand term of my equation ?

Answer 1

PointSource is applied to a vector, similar to DirichletBC. For example,

a = inner(grad(u), grad(v)) * dx
L = Constant(0.) * v * dx
f = PointSource(V, p)

A = assemble(a)
bc.apply(A)

b = assemble(L)
f.apply(b)
bc.apply(b)

uh = Function(V)
x = uh.vector()

solve(A, x, b)

For your other question, concerning a function f that assumes values 1./3 at the vertices of the element containing p, consider the following example.

cell_id, _ = mesh.bounding_box_tree().compute_closest_entity(p)
cell_dofs = V.dofmap().cell_dofs(cell_id)

f = Function(V)
f.vector()[cell_dofs] += 1./3

This of course assumes that V is a P1 space.