# PointSource Method Malfunctioning

**URL:** <https://fenicsproject.discourse.group/t/pointsource-method-malfunctioning/14465>\
**Category:** General\
**Created:** [April 27, 2024, 6:23pm UTC](https://fenicsproject.discourse.group/t/pointsource-method-malfunctioning/14465 "2024-04-27T18:23:36Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![zaccandels](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/zaccandels/32/6204_2.png) [@zaccandels](https://fenicsproject.discourse.group/u/zaccandels)\
**Post date:** [April 27, 2024, 6:23pm UTC](https://fenicsproject.discourse.group/t/pointsource-method-malfunctioning/14465/1 "2024-04-27T18:23:36Z")

</div>

Hi guys,

I’m using FEniCS 2019.1.0 to test the convergence of finite element methods for problems with singular data. Specifically, I’m looking at the Poisson equation with a delta source - \nabla^2u = \delta(x).

Using a mollified delta, I’m able to achieve convergence to the analytical solution u = -\frac{1}{2\pi}\ln(r) for any point except the origin. However, when using the built-in FEniCS PointSource method, I can’t get convergence at all. Does anyone have any idea why this is? Find attached the code I’m using.

```auto
import mshr as mshr
import fenics as fe
import matplotlib.pyplot as plt
import numpy as np

plt.close('all')
class Delta(fe.UserExpression):
    def __init__ (self, eps, x0, **kwargs):
        self.eps = eps
        self.x0 = x0
        fe.UserExpression. __init__ (self, **kwargs)

    def eval(self, values, x):
        eps = self.eps
        values[0] = eps/fe.pi/(np.linalg.norm(x-self.x0) **2 + eps** 2)

    def value_shape(self): return ()

resolution = [5, 10, 20, 40]
error_moll_delta = np.zeros(4)
error_pt_source = np.zeros(4)
h = np.zeros(4)
order_moll_delta = np.zeros(4)
order_pt_source = np.zeros(4)
for i in range(len(resolution)):
    domain_geometry = mshr.Circle(fe.Point(0., 0.), 1 + 0.01)
    mesh = mshr.generate_mesh(domain_geometry, resolution[i])
    h[i] = mesh.hmax()
    V = fe.FunctionSpace(mesh, "Lagrange", 1)
    
    
    def boundary(x, on_boundary):
        return on_boundary
    
    
    # Define boundary condition
    u0 = fe.Constant(0.0)
    bc = fe.DirichletBC(V, u0, boundary)
    
    #%% Using Continuous Delta Approximation
    
    # Define variational problem
    u = fe.TrialFunction(V)
    v = fe.TestFunction(V)
    delta = Delta(eps=1E-4, x0=np.array([0., 0.]), degree=5)
    f = delta
    magnitude = fe.Constant(1)
    
    a = fe.dot( fe.grad(u), fe.grad(v) )*fe.dx
    L = magnitude*f*v*fe.dx
    
    u = fe.Function(V)
    
    assembler = fe.SystemAssembler(a, L, bc)
    solver = fe.LUSolver("mumps")
    A = fe.Matrix()
    solver.set_operator(A)
    assembler.assemble(A)
    
    b = fe.assemble(L)
    bc.apply(b)
    fe.solve(A, u.vector(), b)
    
    
    u_e = fe.Expression('-1/(2*3.14159)*'
                        + 'std::log( pow( pow(x[0],2) + pow(x[1],2), 1/2) )',
                        degree = 6)
    V = u.function_space()
    mesh = V.mesh()
    degree = V.ufl_element().degree()
    W = fe.FunctionSpace(mesh, "P", degree + 3)
    u_e_w = fe.interpolate(u_e, W)
    u_w = fe.interpolate(u, W)
    e_w = fe.Function(W)
    e_w.vector()[:] = u_e_w.vector().get_local() - u_w.vector().get_local()
    error = e_w**2*fe.dx
    
    E = fe.errornorm(u_e, u, 'L2')
    print("Error using mollified delta with resolution ",
          resolution[i], " is ", E)
    error_moll_delta[i] = E
    
    #%% Using PointSource Method
    
    # Define variational problem
    u = fe.TrialFunction(V)
    v = fe.TestFunction(V)
    f = fe.Constant(0)
    
    a = fe.dot( fe.grad(u), fe.grad(v) )*fe.dx
    L = f*v*fe.dx
    
    u = fe.Function(V)
    
    assembler = fe.SystemAssembler(a, L, bc)
    solver = fe.LUSolver("mumps")
    A = fe.Matrix()
    solver.set_operator(A)
    assembler.assemble(A)
    
    b = fe.assemble(L)
    delta = fe.PointSource(V, fe.Point(0,0), 1)
    delta.apply(b)
    bc.apply(b)
    fe.solve(A, u.vector(), b)
    
    u_e = fe.Expression('-1/(2*3.14159)*'
                        + 'std::log( pow( pow(x[0],2) + pow(x[1],2), 1/2) )',
                        degree = 6)
    V = u.function_space()
    mesh = V.mesh()
    degree = V.ufl_element().degree()
    W = fe.FunctionSpace(mesh, "P", degree + 3)
    u_e_w = fe.interpolate(u_e, W)
    u_w = fe.interpolate(u, W)
    e_w = fe.Function(W)
    e_w.vector()[:] = u_e_w.vector().get_local() - u_w.vector().get_local()
    error = e_w**2*fe.dx
    
    # Also plot analytical solution and error function
    
    E = fe.errornorm(u_e, u, 'L2')
    error_pt_source[i] = E
    print("Error using point source with resolution ", resolution[i], " is ", E)
    
for j in range(1, len(resolution)):
    h_ratio = h[j]/h[j-1]
    error_ratio_pt_source = error_pt_source[j]/error_pt_source[j-1]
    
    error_ratio_moll_delta = error_moll_delta[j]/error_moll_delta[j-1]
    
    order_moll_delta[j] = np.log(error_ratio_moll_delta)/np.log(h_ratio)
    order_pt_source[j] = np.log(error_ratio_pt_source)/np.log(h_ratio)

```

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [April 28, 2024, 5:25am UTC](https://fenicsproject.discourse.group/t/pointsource-method-malfunctioning/14465/2 "2024-04-28T05:25:25Z")

</div>

Seems to me a duplicate of [Non-Convergence When Using 'PointSource' Method](https://fenicsproject.discourse.group/t/non-convergence-when-using-pointsource-method/14221), hence I am closing. If I am wrong, please get in touch with me.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [April 28, 2024, 5:25am UTC](https://fenicsproject.discourse.group/t/pointsource-method-malfunctioning/14465/3 "2024-04-28T05:25:40Z")

</div>