# Dirac point source class implementation for Poisson equation

Hello,

I have been trying to adapt the proposed Delta (source point) class by “MiroK” in this post Delta variational formulation, but I keep getting a shape mismatch error. My code is :

``````class Delta(UserExpression):
def __init__(self, eps, x0, **kwargs):
self.eps = eps
self.x0 = x0
UserExpression.__init__(self, **kwargs)
def eval(self, values, x):
eps = self.eps
values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
values[1] = 0
values[2] = 0

def value_shape(self): return (3, )
``````
``````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)
delta = Delta(eps=1E-4, x0=np.array([0.5, 0, 0]), degree=5)
f = delta
magnitude = Constant(100)
g = Constant(0.0)
L = inner(magnitude*f, v)*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

plot(u)
plt.colorbar(plot(u))
``````

The class was initialy done for a 3 dimensions problem. Since i’m working with a 2D problem, how can I modify the class, or maybe adapt my code to it.

The error I keep getting is on line : `L = inner(magnitude*f, v)*dx`
Saying : shapes do not match.

Thank you.

Hi

The error comes from the fact that your are solving a 2D problem and using a scalar source. You have to redefine the Class to be scalar-type. The following works and impose a point source at (0.5,0.5)

``````from mshr import *
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

class Delta(UserExpression):
def __init__(self, eps, x0, **kwargs):
self.eps = eps
self.x0 = x0
UserExpression.__init__(self, **kwargs)

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

def value_shape(self): return ()

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)
delta = Delta(eps=1E-4, x0=np.array([0.5, 0.5]), degree=5)
f = delta
magnitude = Constant(100)
g = Constant(0.0)