# Implement point source in the variational formulation

Dear FEniCS community,

I’m having a hard time trying to implement a point load on a 3D problem in 1D mesh. I have read old fenics post, i.e., https://fenicsproject.org/qa/14258/how-to-apply-pointsource-to-linearvariationalproblem/ or https://fenicsproject.org/qa/7941/applying-dirac-function-using-fenics-pointsource-function/ , with no luck.

The sample code is as simple as the following:

``````DOF = 50

# Create mesh and define function space
L=4
mesh = IntervalMesh(DOF, 0, L)

degree = 1
U = VectorFunctionSpace(mesh, 'CG', degree, dim=3)

u = TrialFunction(U)

v = TestFunction(U)

# Dirichlet on left boundary
def l_boundary(x, on_boundary):
return on_boundary and near(x, 0.0)

zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(U, zero, l_boundary)

# Variational formulation

# I want the following load f to be a point load on the endpoint L=4
L = inner(f, v)*dx

# Solution to the problem
u_h = Function(U)
solve(A == L,u_h,bc)
``````

From the post related above, a suggested solution was that an approximation of the Dirac delta function can do the work, but I have no idea on how to do this on FEniCS. I tried the solution from Mirok but get an Maximum recursion depth exceeded error.

Thanks in advance for any help!.

The `Expression` approach should be modified as follows

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

def value_shape(self): return (3, )
``````
3 Likes

Dear Mirok,

Thank your for your solution. The code runs succesfully if I set delta=Delta(1E-4) and then

L=inner(delta,v)*dx

But I don’t know how this works, i.e, ¿how to set the magnitude and location of the source given that location is x=4 and source is f=(0,0,4)?

MiroK has outlined the code for a mollified Dirac delta function. See here.

For the precise formulation of a point source in your variational formulation, and FEniCS code, use the `PointSource` class documented here.

Dear nate,

I missed to point out in my question that I’m trying NOT to use the PointSource class in order to avoid the assemble() and apply() part. That will mess up my code…

However, thanks for the documentations.

Changing location mounts to

``````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 = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
values = 0
values = 0

def value_shape(self): return (3, )

delta = Delta(eps=1E-4, x0=np.array([0.5, 0, 0]), degree=5)
``````

You can control magnitude by scaling, i.e. `L = inner(Constant(100)*delta, v)*dx`

4 Likes

Dear MiroK,

Thank you for your valuable help. The code works successfully. For a 1D mesh, we can call Delta as:

``delta=Delta(eps=1E-4, x0=0.5, degree=5)``

Hello @MiroK Mirok, what is the magnitude of the force and location in this code?
I want to add a point load at the free end of a cantilever, is this code useful for that?
Thanks.

Hi Mirok, I have a question to ask you. When I run this code, it returns this error. Is this a problem with my “fenics_adjoint”, or is it something else?

Using user-expressions in dolfin_adjoint requires Extra attention, see

You also need to actually post your code, encapsulated with 3x`, as

`````````python
from fenics import *
```
``````

and make sure that the code is a minimal reproducible example.

1 Like

Thank you very much for your help, I will read these and find my mistakes.

``````deltas = [Delta(eps=1E-4, x0=np.array([0.5, 0, 0]), degree=5),