Is it possible to apply many point sources to b without a python for loop?

Question

I need to do the python loop below for each time step in a transient solution with many >100 point sources. If there a way of applying them all at the same time given a list of points without having to use point.apply(b) for each one? Could I use c++ code in the python script?

mesh_obj = UnitSquareMesh(100,100)
V = FunctionSpace( mesh_obj, "CG", 1 )
v = TestFunction( V )
L = Constant(0.0) * v * dx
b  = assemble ( L )

for pt in list_of_points:
    PointSource( V, Point( pt ), 1.0 ).apply( b )
    # Use this for some calculation.

Answer 1

Yes it is possible to apply the point sources all at the same time. You need to provide a list of tuples each which contain the point and magnitude of the point source and then apply.

Here is some sample code.

mesh_obj = UnitSquareMesh(10,10)
V = FunctionSpace( mesh_obj, "CG", 1 )
v = TestFunction( V )
L = Constant(0.0) * v * dx
b  = assemble ( L )
point_sources = [(Point(0.5, 0.5), 1.0), (Point(0.3, 0.7), 1.0)]

ps = PointSource(V, point_sources)
ps.apply(b)

If you are running in parallel, you need to make sure that the points are either all provided on one processor or only locally on the processor that own the points. This stops them being added twice.

Comments

Thank you Ettie thats exactly what I was looking for. For the parallel point source issue, do you have a sample code?

---

If you just have a list and don't know which processor owns them, the easiest way is just to choose one of the processors to add it on, in this example it is the rank 0 processor. The function then works out which processor owns them to add them to the correct part of the matrix.

mesh_obj= UnitSquareMesh(1,1)
V = FunctionSpace( mesh_obj, "CG", 1 )
v = TestFunction( V )
L = Constant(0.0) * v * dx
b  = assemble ( L )  
point_sources = [];

if MPI.rank(mesh_obj.mpi_comm()) == 0:
     point_sources = [(Point(0.5, 0.5), 1.0), (Point(0.3, 0.7), 1.0)]

ps = PointSource(V, point_sources)
ps.apply(b)

If you know which ones are on each processor you can have a different list for each of them. If you look in the unit tests for the PointSource function there are some tests showing some more examples of how to add them in parallel.

---

Hi Ettie, I'm trying to put all the points and magnitudes for the point source on a list like this:

point_sources = []
nSources = int(l/l_i)
edge     = (l-l_i*nSources)/2.0

for nx in range(0, nSources+1): 
    for ny in range(0, nSources+1): 
    for nz in range(0, nSources+1):

        x = -l/2.0 + edge + nx*l_i
        y = -l/2.0 + edge + ny*l_i
        z = -l/2.0 + edge + nz*l_i

        if np.sqrt(x*x + y*y + z*z) < d_i/2.0:
        point_sources.append((Point(x, y, z), ps_mag))  

delta = PointSource(V, point_sources)

But I get the following error:

File "main.py", line 89, in
delta = PointSource(V, point_sources)
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/fem.py", line 1744, in init
_fem.PointSource_swiginit(self, _fem.new_PointSource(V, p, magnitude))
TypeError: in method 'new_PointSource', argument 2 of type 'dolfin::Point const &'
Aborted (core dumped)

Do you know what the problem is?

---

The code you posted is not working, I run the exact same code and I get the same error:

Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
File "test.py", line 23, in
ps = PointSource(V, point_sources)
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/fem.py", line 1744, in init
_fem.PointSource_swiginit(self, _fem.new_PointSource(V, p, magnitude))
TypeError: in method 'new_PointSource', argument 2 of type 'dolfin::Point const &'
Aborted (core dumped)

I'm using dolfin version 2016.2.0

---

This feature was only added to the development version in the past month or so. If you are using the stable release it won't be there until a new release is made. I would suggest trying the development version.

---

So you mean the PointSource function accepting a list is not in the stable release? It works for one point source. I will try the development version thank you

---

You need to use the development version of DOLFIN available from the bitbucket repository.