Some issues with PointSource
Affects  Status  Importance  Assigned to  Milestone  

 DOLFIN 
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(
The function space no longer exists when apply() is called.
Additionally, I suggest that PointSource inherits BoundaryCondition (with noop LHS apply). This would allow PointSource to be used with the highlevel 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(
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 byhand 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 :  #1 
Marie Rognes (megsimula) wrote :  #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?
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:/
>
> Title:
> Some issues with PointSource
>
> To manage notifications about this bug go to:
> https:/
Anders Logg (logg) wrote :  #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 codesnippet)
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_
interior_
dirac_integral / point_integral / ???

Anders
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 codesnippet)
>
> 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_
> interior_
> dirac_integral / point_integral / ???
>
> 
> Anders
>
> _______
> Mailing list: https:/
> Post to : <email address hidden>
> Unsubscribe : https:/
> More help : https:/
>
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 codesnippet)
>
> 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_
> interior_
> dirac_integral / point_integral / ???
Saved as blueprint:
https:/

Anders
Joachim Haga (jobh) wrote :  #7 
Good! I'll rephrase the bug report to contain issues 1&3 only, those are quite concrete and easy.
description:  updated 
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):q.vector())
"""Unit area delta function in discrete space V"""
q = Function(V)
PointSource(V, pt).apply(
q /= assemble(q*dx)
return q
...
q_s = delta(S, Point(0.0))