Some issues with PointSource

[Joachim Haga]

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)]

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)]

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)]

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.
> https://bugs.launchpad.net/bugs/1010850
>
> Title:
>  Some issues with PointSource
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1010850/+subscriptions

[Anders Logg (logg)]

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

This would require adding an addition integral to UFC:

  cell_integral
  exterior_facet_integral
  interior_facet_integral
  dirac_integral / point_integral / ???

--
Anders

[Joachim Haga (jobh)]

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
>
> This would require adding an addition integral to UFC:
>
> cell_integral
> exterior_facet_integral
> interior_facet_integral
> dirac_integral / point_integral / ???
>
> --
> Anders
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

[Anders Logg (logg)]

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
>
> This would require adding an addition integral to UFC:
>
> cell_integral
> exterior_facet_integral
> interior_facet_integral
> dirac_integral / point_integral / ???

Saved as blueprint:

https://blueprints.launchpad.net/dolfin/+spec/dirac-measure-point-source

--
Anders

[Joachim Haga (jobh)]

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