Use PointSource as a source term (f) for Poisson equation

Legacy DOLFIN
1

Hello,

I am trying to use a Dirac source term with the PointSource function, based on a post in the fenics forum PointSource. Here is my code :

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.0, 0.0))
g = Constant(0.0)
a = inner(grad(u), grad(v))*dx
L = f*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, b)

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

But I keep getting an error on line f = PointSource(V, Point(0.0, 0.0)) saying that object of type ‘dolfin.cpp.geometry.Point’ has no len()

Is it possible that the PointSource function is used in another way for dolfin 2019.1.0 ?

Thank you.

2

Hi

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)

Cheers

2 Likes
3

Thank you again, you saved me huge time ! :slight_smile:

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.

Have a great day !

4

True.

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.

6

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:

0.7842512789122973
0.08417347307356629
0.13157524801413645
0.6151517482997229
0.17422806288686354
0.2106201888134136

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.

Cheers

1 Like