Projecting a Gaussian distribution

On 14 September 2011 13:05, Josef Höök
<email address hidden> wrote:
> New question #171106 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/171106
>
> I would like to project a particle cloud represented as a sum of dirac masses on a grid using PointSource. I have tried the following code but it gives me a set of points on the mesh instead of a nice Gaussian shaped function, any ideas why?
> Perhaps I'm not using the PointSource correctly?

I think it depends on how PointSource should work. Reading the
docstring of PointSource.h I think it should work as you expect, but
this means that the line:

  b.set(&values[0], V->element().space_dimension(), &dofs[0]);

should be changed to

  b.add(&values[0], V->element().space_dimension(), &dofs[0]);

(which makes your code create a nice distribution)

Let's hear Anders' opinion on PointSource functionality.

Kristian

> Cheers
> /J
>
>
> from dolfin import *
> import random
>
> mesh = Rectangle(-10,-10, 10,10,50,50)
> V = FunctionSpace(mesh, 'CG', 1)
> u = TrialFunction(V)
> v = TestFunction(V)
> a = inner(u,v)*dx
> A=assemble(a)
> L= v*dx
> b = assemble(L)
> n = 100000
> for i in range(n):
>    p= Point(random.normalvariate(0,1), random.normalvariate(0,1))
>    delta = PointSource(V, p, 1)
>    delta.apply(b)
>
> u = Function(V)
> solve(A, u.vector(), b)
>
>
>
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Hi,

Modifying PointSource.cpp, line 92, from set(..) to add(..) solves the problem (see attached patch).

Cheers
/Josef

________________________________________
Från: <email address hidden> [<email address hidden>] f&#246;r Josef Höök [<email address hidden>]
Skickat: den 14 september 2011 13:05
Till: Josef Höök
Ämne: [Question #171106]: Projecting a Gaussian distribution

New question #171106 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/171106

I would like to project a particle cloud represented as a sum of dirac masses on a grid using PointSource. I have tried the following code but it gives me a set of points on the mesh instead of a nice Gaussian shaped function, any ideas why?
Perhaps I'm not using the PointSource correctly?

Cheers
/J

from dolfin import *
import random

mesh = Rectangle(-10,-10, 10,10,50,50)
V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u,v)*dx
A=assemble(a)
L= v*dx
b = assemble(L)
n = 100000
for i in range(n):
    p= Point(random.normalvariate(0,1), random.normalvariate(0,1))
    delta = PointSource(V, p, 1)
    delta.apply(b)

u = Function(V)
solve(A, u.vector(), b)

--
You received this question notification because you asked the question.

Hi,

It would be nice to have this feature in Dolfin since projection of particle cloud is very common in for instance particle-in-cell simulations.

Cheers
/J

On Wednesday September 14 2011 04:45:40 Kristian B. Ølgaard wrote:
> Question #171106 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171106
>
> Status: Open => Answered
>
> Kristian B. Ølgaard proposed the following answer:
> On 14 September 2011 13:05, Josef Höök
>
> <email address hidden> wrote:
> > New question #171106 on DOLFIN:
> > https://answers.launchpad.net/dolfin/+question/171106
> >
> > I would like to project a particle cloud represented as a sum of dirac
> > masses on a grid using PointSource. I have tried the following code but
> > it gives me a set of points on the mesh instead of a nice Gaussian
> > shaped function, any ideas why? Perhaps I'm not using the PointSource
> > correctly?
>
> I think it depends on how PointSource should work. Reading the
> docstring of PointSource.h I think it should work as you expect, but
> this means that the line:
>
> b.set(&values[0], V->element().space_dimension(), &dofs[0]);
>
> should be changed to
>
> b.add(&values[0], V->element().space_dimension(), &dofs[0]);
>
> (which makes your code create a nice distribution)

I think this is the natural thing to do. Otherwise one cannot combine a
pointsource with any other source term.

Johan

> Let's hear Anders' opinion on PointSource functionality.
>
> Kristian
>
> > Cheers
> > /J
> >
> >
> > from dolfin import *
> > import random
> >
> > mesh = Rectangle(-10,-10, 10,10,50,50)
> > V = FunctionSpace(mesh, 'CG', 1)
> > u = TrialFunction(V)
> > v = TestFunction(V)
> > a = inner(u,v)*dx
> > A=assemble(a)
> > L= v*dx
> > b = assemble(L)
> > n = 100000
> > for i in range(n):
> > p= Point(random.normalvariate(0,1), random.normalvariate(0,1))
> > delta = PointSource(V, p, 1)
> > delta.apply(b)
> >
> > u = Function(V)
> > solve(A, u.vector(), b)
> >
> >
> >
> >
> >
> > --
> > You received this question notification because you are a member of
> > DOLFIN Team, which is an answer contact for DOLFIN.

I've push the fix to the development version.

B0;268;0cOn Wed, Sep 14, 2011 at 11:45:40AM -0000, Kristian B. Ølgaard wrote:
> Question #171106 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/171106
>
> Status: Open => Answered
>
> Kristian B. Ølgaard proposed the following answer:
> On 14 September 2011 13:05, Josef Höök
> <email address hidden> wrote:
> > New question #171106 on DOLFIN:
> > https://answers.launchpad.net/dolfin/+question/171106
> >
> > I would like to project a particle cloud represented as a sum of dirac masses on a grid using PointSource. I have tried the following code but it gives me a set of points on the mesh instead of a nice Gaussian shaped function, any ideas why?
> > Perhaps I'm not using the PointSource correctly?
>
> I think it depends on how PointSource should work. Reading the
> docstring of PointSource.h I think it should work as you expect, but
> this means that the line:
>
> b.set(&values[0], V->element().space_dimension(), &dofs[0]);
>
> should be changed to
>
> b.add(&values[0], V->element().space_dimension(), &dofs[0]);
>
> (which makes your code create a nice distribution)
>
> Let's hear Anders' opinion on PointSource functionality.

Looks good. I see that it has been fixed already.

--
Anders

> Kristian
>
> > Cheers
> > /J
> >
> >
> > from dolfin import *
> > import random
> >
> > mesh = Rectangle(-10,-10, 10,10,50,50)
> > V = FunctionSpace(mesh, 'CG', 1)
> > u = TrialFunction(V)
> > v = TestFunction(V)
> > a = inner(u,v)*dx
> > A=assemble(a)
> > L= v*dx
> > b = assemble(L)
> > n = 100000
> > for i in range(n):
> >    p= Point(random.normalvariate(0,1), random.normalvariate(0,1))
> >    delta = PointSource(V, p, 1)
> >    delta.apply(b)
> >
> > u = Function(V)
> > solve(A, u.vector(), b)
> >
> >
> >
> >
> >
> >
>