Extracting data from specific points

If you are solving with a MixedFunctionSpace ala TaylorHood you can extract
the sub function for the pressure from the mixed solution function. Something
like:

  P = u.split(deepcopy=True)[1]

This assumes the pressure is the last sub function of two. The degrees of
freedom (the solved pressure values) can then be read out directly from the
vector of the functions.

  values = P.vector()

I hope this answered your question.

Johan

On Wednesday May 18 2011 02:55:54 Kyle wrote:
> New question #158012 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/158012
>
> I am looking to extract data from specific nodes on a solved mesh in FEniCS
> (used stokes, want to extract pressure). I want to do this so I can derive
> stresses and apply it to a structure mesh and create a simplistic FSI
> model. I have been looking through fsiproblem.py in cbc.solve to get some
> ideas but have not had any luck.
>
> Thanks!
>
> Kyle

Hi Johan,

Thanks this helped me a bit, so now I have

solver.solve(problem,s.vector())
(U, P) = s.split(deepcopy=True)

values = P.vector()

This is working fine it seems, but I am still unsure how to extract specific data points from 'values' say from x[0.3 to 0.5] & y[0.0 to 0.3] as an example

On 19 May 2011 14:15, Kyle <email address hidden> wrote:
> Question #158012 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/158012
>
>    Status: Answered => Open
>
> Kyle is still having a problem:
> Hi Johan,
>
> Thanks this helped me a bit, so now I have
>
> solver.solve(problem,s.vector())
> (U, P) = s.split(deepcopy=True)
>
> values = P.vector()
>
> This is working fine it seems, but I am still unsure how to extract
> specific data points from 'values' say from x[0.3 to 0.5] & y[0.0 to
> 0.3] as an example

I think you need to:
  get the dofmap of P,
  loop the cells in the mesh,
    tabulate the dofs on each cell
    tabulate the coordinates of each dof on the cell
    check which dofs satisfy your criterion and then extract the
values from the P vector using the dofs

Kristian

> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp
>

On Thu, May 19, 2011 at 02:30:18PM +0200, Kristian Ølgaard wrote:
> On 19 May 2011 14:15, Kyle <email address hidden> wrote:
> > Question #158012 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/158012
> >
> >    Status: Answered => Open
> >
> > Kyle is still having a problem:
> > Hi Johan,
> >
> > Thanks this helped me a bit, so now I have
> >
> > solver.solve(problem,s.vector())
> > (U, P) = s.split(deepcopy=True)
> >
> > values = P.vector()
> >
> > This is working fine it seems, but I am still unsure how to extract
> > specific data points from 'values' say from x[0.3 to 0.5] & y[0.0 to
> > 0.3] as an example
>
> I think you need to:
> get the dofmap of P,
> loop the cells in the mesh,
> tabulate the dofs on each cell
> tabulate the coordinates of each dof on the cell
> check which dofs satisfy your criterion and then extract the
> values from the P vector using the dofs

Why not just call U(x), P(x) with the points in question?

--
Anders

If I used P(x) how could I extract points specifically at the nodes on the boundary at x=1 for example?

On Thu, May 19, 2011 at 03:05:58PM -0000, Kyle wrote:
> Question #158012 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/158012
>
> Status: Answered => Open
>
> Kyle is still having a problem:
> If I used P(x) how could I extract points specifically at the nodes on
> the boundary at x=1 for example?

Do you want the vertices only or the nodes for higher order elements?

--
Anders

On May 19, 2011, at 5:31, Kristian B. Ølgaard<email address hidden> wrote:

> Question #158012 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/158012
>
> Status: Open => Answered
>
> Kristian B. Ølgaard proposed the following answer:
> On 19 May 2011 14:15, Kyle <email address hidden> wrote:
>> Question #158012 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/158012
>>
>> Status: Answered => Open
>>
>> Kyle is still having a problem:
>> Hi Johan,
>>
>> Thanks this helped me a bit, so now I have
>>
>> solver.solve(problem,s.vector())
>> (U, P) = s.split(deepcopy=True)
>>
>> values = P.vector()
>>
>> This is working fine it seems, but I am still unsure how to extract
>> specific data points from 'values' say from x[0.3 to 0.5] & y[0.0 to
>> 0.3] as an example
>
> I think you need to:
> get the dofmap of P,
> loop the cells in the mesh,
> tabulate the dofs on each cell
> tabulate the coordinates of each dof on the cell
> check which dofs satisfy your criterion and then extract the
> values from the P vector using the dofs

I think such a function would be very useful. Would it be possible to put a free function of this into DOLIN. It could take a FunctionSpace and a subdomain/MeshFunction as input and return an Array<uint> of dofs indices.

This would then hopefully also just work in parallel.

Johan

> Kristian
>
>> --
>> You received this question notification because you are a member of
>> DOLFIN Team, which is an answer contact for DOLFIN.
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~dolfin
>> Post to : <email address hidden>
>> Unsubscribe : https://launchpad.net/~dolfin
>> More help : https://help.launchpad.net/ListHelp
>>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp

I want the pressure and velocity values at all nodes along a boundary so that I can compute the stress tensor and multiply it by the stress normal. With the computed stress I would like to apply them to their corresponding node on the solid mesh to solve the elasticity equation.

On Thu, May 19, 2011 at 03:46:03PM -0000, Kyle wrote:
> Question #158012 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/158012
>
> Status: Answered => Open
>
> Kyle is still having a problem:
> I want the pressure and velocity values at all nodes along a boundary so
> that I can compute the stress tensor and multiply it by the stress
> normal. With the computed stress I would like to apply them to their
> corresponding node on the solid mesh to solve the elasticity equation.

Yes, but which type of function space do you have? If you have
piecewise linears, you only need to get the values at the vertices
(which is easy). If you have quadratics, then you need the edges,
which is also easy. If you have higher order, then it gets more tricky
and you would have to follow the recipee by Kristian.

If you're interested, you can have a look in the cbc.solve code
(available on Launchpad) where there is a simple FSI solver that does
exactly this: computing the stress from a P2-P1 Taylor-Hood
approximation of the fluid field and applying the force to a P1
structure problem.

--
Anders

Oh sorry!

Yes I am using the Taylor-Hood elements. Is this the code located in fsiproblem.py in sandbox/fsi? I looked through it but it looked a little larger than I thought it was going to be.

Kyle

On Thu, May 19, 2011 at 04:35:57PM -0000, Kyle wrote:
> Question #158012 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/158012
>
> Status: Answered => Open
>
> Kyle is still having a problem:
> Oh sorry!
>
> Yes I am using the Taylor-Hood elements.

Then just do something like this:

for v in vertices(mesh)
  ... = U(v.point())

for e in edges(mesh)
  ... = U(e.midpoint())

I noticed now that function evaluation does not work for the Point
class. It requires instead the coordinates as input so you need to do
something like

  p = v.point()
  ... = U((v.x(), v.y()))

instead.

Note to Hake: Would it be easy to add the possibility of supplying a
Point argument to __call__?

> Is this the code located in
> fsiproblem.py in sandbox/fsi? I looked through it but it looked a little
> larger than I thought it was going to be.

Take a look in the file subproblems.py and search for "transfer". You
might also want to look in fsiproblem.py for how the mapping between
vertex numbers and edge numbers in the two meshes are computed.

--
Anders

Sure!

I just added the request as comment to the bug:

 "Lacking check of input for function evaluation"

So I remember it when I fix that one :)

Johan

On Thursday May 19 2011 10:17:43 Anders Logg wrote:
> On Thu, May 19, 2011 at 04:35:57PM -0000, Kyle wrote:
> > Question #158012 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/158012
> >
> > Status: Answered => Open
> >
> > Kyle is still having a problem:
> > Oh sorry!
> >
> > Yes I am using the Taylor-Hood elements.
>
> Then just do something like this:
>
> for v in vertices(mesh)
> ... = U(v.point())
>
> for e in edges(mesh)
> ... = U(e.midpoint())
>
> I noticed now that function evaluation does not work for the Point
> class. It requires instead the coordinates as input so you need to do
> something like
>
> p = v.point()
> ... = U((v.x(), v.y()))
>
> instead.
>
> Note to Hake: Would it be easy to add the possibility of supplying a
> Point argument to __call__?
>
> > Is this the code located in
> > fsiproblem.py in sandbox/fsi? I looked through it but it looked a little
> > larger than I thought it was going to be.
>
> Take a look in the file subproblems.py and search for "transfer". You
> might also want to look in fsiproblem.py for how the mapping between
> vertex numbers and edge numbers in the two meshes are computed.
>
> --
> Anders
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp

Great!

--
Anders

On Thu, May 19, 2011 at 10:37:24AM -0700, Johan Hake wrote:
> Sure!
>
> I just added the request as comment to the bug:
>
> "Lacking check of input for function evaluation"
>
> So I remember it when I fix that one :)
>
> Johan
>
>
> On Thursday May 19 2011 10:17:43 Anders Logg wrote:
> > On Thu, May 19, 2011 at 04:35:57PM -0000, Kyle wrote:
> > > Question #158012 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/158012
> > >
> > > Status: Answered => Open
> > >
> > > Kyle is still having a problem:
> > > Oh sorry!
> > >
> > > Yes I am using the Taylor-Hood elements.
> >
> > Then just do something like this:
> >
> > for v in vertices(mesh)
> > ... = U(v.point())
> >
> > for e in edges(mesh)
> > ... = U(e.midpoint())
> >
> > I noticed now that function evaluation does not work for the Point
> > class. It requires instead the coordinates as input so you need to do
> > something like
> >
> > p = v.point()
> > ... = U((v.x(), v.y()))
> >
> > instead.
> >
> > Note to Hake: Would it be easy to add the possibility of supplying a
> > Point argument to __call__?
> >
> > > Is this the code located in
> > > fsiproblem.py in sandbox/fsi? I looked through it but it looked a little
> > > larger than I thought it was going to be.
> >
> > Take a look in the file subproblems.py and search for "transfer". You
> > might also want to look in fsiproblem.py for how the mapping between
> > vertex numbers and edge numbers in the two meshes are computed.
> >