You are invoking PointSource wrong. You need to call PointSource(V,Point,value). Also, as the documentation says here you need to use apply() to add PointSource to the RHS. There is no point of using it on L.

The following runs with no errors.

from mshr import *
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
resolution = 15
domain_geometry = Circle(dolfin.Point(0.5, 0.5), 0.5+0.01)
mesh = generate_mesh(domain_geometry, resolution)
V = FunctionSpace(mesh, "Lagrange", 1)
def boundary(x, on_boundary):
return on_boundary
# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = PointSource(V, Point(0.4, 0.5), 1)
g = Constant(0.0)
a = inner(grad(u), grad(v))*dx
L = Constant(0.0)*v*dx
A = assemble(a)
bc.apply(A)
b = assemble(L)
f.apply(b)
bc.apply(b)
# Compute solution
u = Function(V)
solve(A, u.vector(), b)

Out of curiosity : Why did you need to use the Delta class instead of PointSource ? It looks like the PointSource is more precise to produce a dirac than the Delta class.

At that moment, I only needed to have an explicit approximation of the Dirac delta to play with the tolerance and observe the behavior in the limit. It was enough for my simulation.

I think it is because the coordinates are not part of the mesh nodes (point source contribution is distributed over neighbor elements). Note that the max values your provided are not the only ones:

Hence there will be an interpolation error when approximating the point source. You could extract the mesh coordinates and use them as point sources. Also you can use a sufficiently refined mesh or use a uniform mesh to have better control over the error.