# Some issues with PointSource

Bug #1010850 reported by Joachim Haga on 2012-06-09
This bug affects 1 person
Affects Status Importance Assigned to Milestone
Undecided
Unassigned

### Bug Description

Edit: The original report is below; some issues in the original are converted to a blueprint.

PointSource has a memory bug, when a temporary is used for the function space:
qu = PointSource(W.sub(1), ...)
The function space no longer exists when apply() is called.

Additionally, I suggest that PointSource inherits BoundaryCondition (with no-op LHS apply). This would allow PointSource to be used with the high-level solve(a==L, bcs=...) interface.

================== original ===================

I have tried to use the PointSource class. However, I've found some problems with it:

1/ It doesn't work properly for mixed (sub)spaces in python:
qu = PointSource(W.sub(1), Point(0.0))
fails with some assertion error, but
Q = W.sub(1); qu=PointSource(Q, Point(0.0))
seems to work.

2/ It is often desired to scale the point source so that it integrates to 1.0 with the given basis (like the delta function does in the continuous case). This can be done by hand, but it would be extra nice to have it done automatically. Especially since such by-hand scaling would require multiple steps (create trial source; assemble to find scale; create scaled source) in the general case.

3/ It can not be used with the solve(A==b) interface in python (I think). Would it be a good idea to let it inherit from BoundaryCondition, so that one can say solve(A==b, bcs=qu)? (but see next question)

4/ Even nicer would be if it could be used in forms (if it was a Function or Expression). That's how it's usually written in mathematical notation. Hm. I guess I can make such a function myself easily enough with apply() on an f.vector(). Still: Is making PointSource a Function, instead of something which is called through apply(), something you would consider?

 Joachim Haga (jobh) wrote on 2012-06-09: #1

Update (for reference): I have solved my problems using method 4 (above), with the small python function below. So I'm happy. Problems 1 and 3 may still be worth fixing, though.
==

def delta(V, pt):
"""Unit area delta function in discrete space V"""
q = Function(V)
PointSource(V, pt).apply(q.vector())
q /= assemble(q*dx)
return q

...
q_s = delta(S, Point(0.0))

 Marie Rognes (meg-simula) wrote on 2012-06-10: #2

Response to 4/

(Caveat: I've never used PointSource)

How about adding a new measure: a point measure, corresponding to evaluation at a given point,
say

dp

Could that resolve the issues at hand?

 Garth Wells (garth-wells) wrote on 2012-06-10: Re: [Bug 1010850] Re: Some issues with PointSource #3

On 10 June 2012 14:21, Marie Rognes <email address hidden> wrote:
> Response to 4/
>
> (Caveat: I've never used PointSource)
>
> How about adding a new measure: a point measure, corresponding to evaluation at a given point,
> say
>
>  dp
>

Supporting the Delta function would be neater. I don't know what the
implementation should look like - it's a little tricky since the user
would typically provide the location in terms of the physical
coordinates.

Garth

> Could that resolve the issues at hand?
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
>
> Title:
>  Some issues with PointSource
>

 Anders Logg (logg) wrote on 2012-06-11: #4

On Sun, Jun 10, 2012 at 01:21:57PM -0000, Marie Rognes wrote:
> Response to 4/
>
> (Caveat: I've never used PointSource)
>
> How about adding a new measure: a point measure, corresponding to evaluation at a given point,
> say
>
> dp
>
> Could that resolve the issues at hand?

I think that would be neat.

We set (formally)

dp(\bar{x}) = \delta_{\bar{x}} dx

and then the generated code would need to

1. Check whether point is in the current triangle (easily handled with
an appropriate code-snippet)

2. Evaluate the integrand at the point in question and add to the
local tensor

cell_integral
exterior_facet_integral
interior_facet_integral
dirac_integral / point_integral / ???

--
Anders

 Joachim Haga (jobh) wrote on 2012-06-11: Re: [Dolfin] [Bug 1010850] Re: Some issues with PointSource #5

I won't pretend that I understand the issues here, but if it works
predictably (even when hitting the vertex coordinate of a discontinuous
space) then I'm all for it ;)

-j.

On 11 June 2012 14:42, Anders Logg <email address hidden> wrote:

> On Sun, Jun 10, 2012 at 01:21:57PM -0000, Marie Rognes wrote:
> > Response to 4/
> >
> > (Caveat: I've never used PointSource)
> >
> > How about adding a new measure: a point measure, corresponding to
> evaluation at a given point,
> > say
> >
> > dp
> >
> > Could that resolve the issues at hand?
>
> I think that would be neat.
>
> We set (formally)
>
> dp(\bar{x}) = \delta_{\bar{x}} dx
>
> and then the generated code would need to
>
> 1. Check whether point is in the current triangle (easily handled with
> an appropriate code-snippet)
>
> 2. Evaluate the integrand at the point in question and add to the
> local tensor
>
>
> cell_integral
> exterior_facet_integral
> interior_facet_integral
> dirac_integral / point_integral / ???
>
> --
> Anders
>
> _______________________________________________
> Post to : <email address hidden>
>

 Anders Logg (logg) wrote on 2012-06-13: Re: [Bug 1010850] Re: Some issues with PointSource #6

On Mon, Jun 11, 2012 at 02:42:51PM +0200, Anders Logg wrote:
> On Sun, Jun 10, 2012 at 01:21:57PM -0000, Marie Rognes wrote:
> > Response to 4/
> >
> > (Caveat: I've never used PointSource)
> >
> > How about adding a new measure: a point measure, corresponding to evaluation at a given point,
> > say
> >
> > dp
> >
> > Could that resolve the issues at hand?
>
> I think that would be neat.
>
> We set (formally)
>
> dp(\bar{x}) = \delta_{\bar{x}} dx
>
> and then the generated code would need to
>
> 1. Check whether point is in the current triangle (easily handled with
> an appropriate code-snippet)
>
> 2. Evaluate the integrand at the point in question and add to the
> local tensor
>
>
> cell_integral
> exterior_facet_integral
> interior_facet_integral
> dirac_integral / point_integral / ???

Saved as blueprint:

--
Anders

 Joachim Haga (jobh) wrote on 2012-06-13: #7

Good! I'll rephrase the bug report to contain issues 1&3 only, those are quite concrete and easy.

Joachim Haga (jobh) on 2012-06-13
 description: updated