tabulate_coordinates missing in python

Asked by Mikael Mortensen

Hi,

It seems that dofmap's tabulate_coordinates method is missing in dolfin 0.9.11 (I notice this i probably since it has changed to take a boost::multi_array as argument and not double**). Is there another (better?) way of accessing these coordinates now?

Thanks for any help

Mikael

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Johan Hake (johan-hake) said :
#1

I tried adding some typemaps for boost::multi_array to NumPy arrays, but I got
into the same trouble as with std::vector typemaps. There is no way we can
create a multi_array from existing data without copying.

If this is fine I can probably add this typemap very soonish.

Johan

> Hi,
>
> It seems that dofmap's tabulate_coordinates method is missing in dolfin
> 0.9.11 (I notice this i probably since it has changed to take a
> boost::multi_array as argument and not double**). Is there another
> (better?) way of accessing these coordinates now?
>
> Thanks for any help
>
> Mikael

Revision history for this message
Mikael Mortensen (mikael-mortensen) said :
#2

Pardon my ignorance, but couldn't you just keep the
tabulate_coordinates(double** coordinates, const Cell& cell) as an
additional option on the C++ side and not ignore this with swig? Seems like
ufc takes the double** anyway? Otherwise I have no problem with copying.

Mikael

On 6 August 2011 06:45, Johan Hake <email address hidden>wrote:

> Your question #167027 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/167027
>
> Status: Open => Answered
>
> Johan Hake proposed the following answer:
> I tried adding some typemaps for boost::multi_array to NumPy arrays, but I
> got
> into the same trouble as with std::vector typemaps. There is no way we can
> create a multi_array from existing data without copying.
>
> If this is fine I can probably add this typemap very soonish.
>
> Johan
>
> > Hi,
> >
> > It seems that dofmap's tabulate_coordinates method is missing in dolfin
> > 0.9.11 (I notice this i probably since it has changed to take a
> > boost::multi_array as argument and not double**). Is there another
> > (better?) way of accessing these coordinates now?
> >
> > Thanks for any help
> >
> > Mikael
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/dolfin/+question/167027/+confirm?answer_id=0
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/dolfin/+question/167027
>
> You received this question notification because you asked the question.
>

Revision history for this message
Garth Wells (garth-wells) said :
#3

On 06/08/11 15:26, Mikael Mortensen wrote:
> Question #167027 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/167027
>
> Status: Answered => Open
>
> Mikael Mortensen is still having a problem:
> Pardon my ignorance, but couldn't you just keep the
> tabulate_coordinates(double** coordinates, const Cell& cell) as an
> additional option on the C++ side and not ignore this with swig? Seems like
> ufc takes the double** anyway?

The double pointer should be removed from UFC.

> Otherwise I have no problem with copying.

I think that copying is fine, for now.

Garth

>
> Mikael
>
> On 6 August 2011 06:45, Johan Hake
> <email address hidden>wrote:
>
>> Your question #167027 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/167027
>>
>> Status: Open => Answered
>>
>> Johan Hake proposed the following answer:
>> I tried adding some typemaps for boost::multi_array to NumPy arrays, but I
>> got
>> into the same trouble as with std::vector typemaps. There is no way we can
>> create a multi_array from existing data without copying.
>>
>> If this is fine I can probably add this typemap very soonish.
>>
>> Johan
>>
>>> Hi,
>>>
>>> It seems that dofmap's tabulate_coordinates method is missing in dolfin
>>> 0.9.11 (I notice this i probably since it has changed to take a
>>> boost::multi_array as argument and not double**). Is there another
>>> (better?) way of accessing these coordinates now?
>>>
>>> Thanks for any help
>>>
>>> Mikael
>>
>> --
>> If this answers your question, please go to the following page to let us
>> know that it is solved:
>> https://answers.launchpad.net/dolfin/+question/167027/+confirm?answer_id=0
>>
>> If you still need help, you can reply to this email or go to the
>> following page to enter your feedback:
>> https://answers.launchpad.net/dolfin/+question/167027
>>
>> You received this question notification because you asked the question.
>>
>

Revision history for this message
Johan Hake (johan-hake) said :
#4

On Saturday August 6 2011 12:46:01 Garth Wells wrote:
> Question #167027 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/167027
>
> Status: Open => Answered
>
> Garth Wells proposed the following answer:
>
> On 06/08/11 15:26, Mikael Mortensen wrote:
> > Question #167027 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/167027
> >
> > Status: Answered => Open
> >
> > Mikael Mortensen is still having a problem:
> > Pardon my ignorance, but couldn't you just keep the
> > tabulate_coordinates(double** coordinates, const Cell& cell) as an
> > additional option on the C++ side and not ignore this with swig? Seems
> > like ufc takes the double** anyway?
>
> The double pointer should be removed from UFC.
>
> > Otherwise I have no problem with copying.
>
> I think that copying is fine, for now.

Ok, then I have a shot at it.

Johan

> Garth
>
> > Mikael
> >
> > On 6 August 2011 06:45, Johan Hake
> >
> > <email address hidden>wrote:
> >> Your question #167027 on DOLFIN changed:
> >> https://answers.launchpad.net/dolfin/+question/167027
> >>
> >> Status: Open => Answered
> >>
> >> Johan Hake proposed the following answer:
> >> I tried adding some typemaps for boost::multi_array to NumPy arrays, but
> >> I got
> >> into the same trouble as with std::vector typemaps. There is no way we
> >> can create a multi_array from existing data without copying.
> >>
> >> If this is fine I can probably add this typemap very soonish.
> >>
> >> Johan
> >>
> >>> Hi,
> >>>
> >>> It seems that dofmap's tabulate_coordinates method is missing in dolfin
> >>> 0.9.11 (I notice this i probably since it has changed to take a
> >>> boost::multi_array as argument and not double**). Is there another
> >>> (better?) way of accessing these coordinates now?
> >>>
> >>> Thanks for any help
> >>>
> >>> Mikael
> >>
> >> --
> >> If this answers your question, please go to the following page to let us
> >> know that it is solved:
> >> https://answers.launchpad.net/dolfin/+question/167027/+confirm?answer_id
> >> =0
> >>
> >> If you still need help, you can reply to this email or go to the
> >> following page to enter your feedback:
> >> https://answers.launchpad.net/dolfin/+question/167027
> >>
> >> You received this question notification because you asked the question.

Revision history for this message
Johan Hake (johan-hake) said :
#5

On Saturday August 6 2011 20:31:17 Johan Hake wrote:
> Question #167027 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/167027
>
> Johan Hake proposed the following answer:
>
> On Saturday August 6 2011 12:46:01 Garth Wells wrote:
> > Question #167027 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/167027
> >
> > Status: Open => Answered
> >
> > Garth Wells proposed the following answer:
> >
> > On 06/08/11 15:26, Mikael Mortensen wrote:
> > > Question #167027 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/167027
> > >
> > > Status: Answered => Open
> > >
> > > Mikael Mortensen is still having a problem:
> > > Pardon my ignorance, but couldn't you just keep the
> > > tabulate_coordinates(double** coordinates, const Cell& cell) as an
> > > additional option on the C++ side and not ignore this with swig? Seems
> > > like ufc takes the double** anyway?
> >
> > The double pointer should be removed from UFC.
> >
> > > Otherwise I have no problem with copying.
> >
> > I think that copying is fine, for now.
>
> Ok, then I have a shot at it.

Done!

I change the order of the arguments, making it optional to pass a NumPy array
for filling.

Johan

> Johan
>
> > Garth
> >
> > > Mikael
> > >
> > > On 6 August 2011 06:45, Johan Hake
> > >
> > > <email address hidden>wrote:
> > >> Your question #167027 on DOLFIN changed:
> > >> https://answers.launchpad.net/dolfin/+question/167027
> > >>
> > >> Status: Open => Answered
> > >>
> > >> Johan Hake proposed the following answer:
> > >> I tried adding some typemaps for boost::multi_array to NumPy arrays,
> > >> but I got
> > >> into the same trouble as with std::vector typemaps. There is no way we
> > >> can create a multi_array from existing data without copying.
> > >>
> > >> If this is fine I can probably add this typemap very soonish.
> > >>
> > >> Johan
> > >>
> > >>> Hi,
> > >>>
> > >>> It seems that dofmap's tabulate_coordinates method is missing in
> > >>> dolfin 0.9.11 (I notice this i probably since it has changed to take
> > >>> a boost::multi_array as argument and not double**). Is there another
> > >>> (better?) way of accessing these coordinates now?
> > >>>
> > >>> Thanks for any help
> > >>>
> > >>> Mikael
> > >>
> > >> --
> > >> If this answers your question, please go to the following page to let
> > >> us know that it is solved:
> > >> https://answers.launchpad.net/dolfin/+question/167027/+confirm?answer_
> > >> id =0
> > >>
> > >> If you still need help, you can reply to this email or go to the
> > >> following page to enter your feedback:
> > >> https://answers.launchpad.net/dolfin/+question/167027
> > >>
> > >> You received this question notification because you asked the
> > >> question.

Revision history for this message
Torbjørn Bækø Ness (torbjorn-ness) said :
#6

Hi,

I'm not able to understand how you are supposed to use this. Say I want to get out the coordinates of all points below, how do I do that?

from dolfin import *
mesh = UnitCube(6, 4, 5)
V = FunctionSpace(mesh, "CG", 2)
dist = 'pow(pow(x[0],2) + pow(x[1],2) + pow(x[2],2) + 0.1, -0.5)'
phiExp = dolfin.Expression(dist)
phi = dolfin.interpolate(phiExp, V)

Revision history for this message
Mikael Mortensen (mikael-mortensen) said :
#7

Not quite sure what you mean, but how about something like

xyz = interpolate(Expression(("x[0]", "x[1]", "x[2]"), value_shape=3), MixedFunctionSpace([V, V, V]))

?

Revision history for this message
Torbjørn Bækø Ness (torbjorn-ness) said :
#8

That was an interesting solution :-)
It seems to work fine. I want to compare the error between two solutions as a function of distance.
Below is an example comparing the difference between P1 and P2 elements for the toy example.

import numpy as np
from dolfin import *
import pylab as pl

mesh = UnitCube(2, 2, 2)
V1 = FunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 2)
dist = 'pow(pow(x[0],2) + pow(x[1],2) + pow(x[2],2) + 0.1, -0.5)'
phiExp = Expression(dist)
coor = mesh.coordinates()

phi1 = interpolate(phiExp, V1).vector().array()
phi2 = interpolate(phiExp, V2).vector().array()
xyz1 = interpolate(Expression(("x[0]", "x[1]", "x[2]"), value_shape=3), MixedFunctionSpace([V1, V1, V1])).vector().array()
xyz2 = interpolate(Expression(("x[0]", "x[1]", "x[2]"), value_shape=3), MixedFunctionSpace([V2, V2, V2])).vector().array()
xyz1 = xyz1.reshape(3,len(xyz1)/3)
xyz2 = xyz2.reshape(3,len(xyz2)/3)
x1, y1, z1 = xyz1[:,:]
x2, y2, z2 = xyz2[:,:]

dist1 = [np.sqrt(x1[i]**2 + y1[i]**2 + z1[i]**2) for i in xrange(len(phi1))]
dist2 = [np.sqrt(x2[i]**2 + y2[i]**2 + z2[i]**2) for i in xrange(len(phi2))]

pl.plot(dist1, phi1, 'Db')
pl.plot(dist2, phi2, 'or')
pl.show()

Revision history for this message
Johan Hake (johan-hake) said :
#9

DofMap.tabulate_coordinates is exposed in the Python interface. Have
look at:

  test/unit/fem/python/DofMap.py

For an example. Note however, that this is a low level interface which
works on a cell basis.

Johan

On 10/10/2012 02:36 PM, Torbjørn Bækø Ness wrote:
> Question #167027 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/167027
>
> Torbjørn Bækø Ness posted a new comment:
> Hi,
>
> I'm not able to understand how you are supposed to use this. Say I want
> to get out the coordinates of all points below, how do I do that?
>
> from dolfin import *
> mesh = UnitCube(6, 4, 5)
> V = FunctionSpace(mesh, "CG", 2)
> dist = 'pow(pow(x[0],2) + pow(x[1],2) + pow(x[2],2) + 0.1, -0.5)'
> phiExp = dolfin.Expression(dist)
> phi = dolfin.interpolate(phiExp, V)
>

Revision history for this message
Torbjørn Bækø Ness (torbjorn-ness) said :
#10

Hi,

Has there been any changes to FEniCS recently that would lead to that the program I supplied above would not give meaningful results anymore?

Spesifically it seems to me that,
xyz2 = interpolate(Expression(("x[0]", "x[1]", "x[2]"), value_shape=3), MixedFunctionSpace([V2, V2, V2])).vector().array()
xyz2 = xyz2.reshape(3,len(xyz2)/3)
does no longer return the coordinates corresponding to the FEM solution?

Torbjørn

Revision history for this message
Mikael Mortensen (mikael-mortensen) said :
#11

Yes,

The UnitCube should now be UnitCubeMesh.

Mikael

Revision history for this message
Mikael Mortensen (mikael-mortensen) said :
#12

Otherwise I think it still works. However, the dofs are not ordered such that you can reshape and get the dofs nicely aligned. The dofs of the first V2 space are not necessarily the first V2.dim() dofs in xyz2.vector().array().

Mikael

Revision history for this message
Torbjørn Bækø Ness (torbjorn-ness) said :
#13

Okey, but how do I find the coordinates of each value in phi2 then, corresponding to mesh.coordinates() for V1?

Revision history for this message
Mikael Mortensen (mikael-mortensen) said :
#14

You need the sub_dofmaps. Try this

phi1 = interpolate(phiExp, V1).vector().array()
phi2 = interpolate(phiExp, V2).vector().array()
VVV = MixedFunctionSpace([V1, V1, V1])
xyz = interpolate(Expression(("x[0]", "x[1]", "x[2]"), value_shape=3), VVV).vector().array()
dofmap = VVV.dofmap()
x = xyz[dofmap.extract_sub_dofmap(np.array([0], "I"), mesh).collapse(mesh)[1].values()]
y = xyz[dofmap.extract_sub_dofmap(np.array([1], "I"), mesh).collapse(mesh)[1].values()]
z = xyz[dofmap.extract_sub_dofmap(np.array([2], "I"), mesh).collapse(mesh)[1].values()]

Revision history for this message
Torbjørn Bækø Ness (torbjorn-ness) said :
#15

Thank you that worked.

Revision history for this message
Mikael Mortensen (mikael-mortensen) said :
#16

Great. This also works by the way

x = xyz[VVV.sub(0).collapse(mesh)[1].values()]
y = xyz[VVV.sub(1).collapse(mesh)[1].values()]
z = xyz[VVV.sub(2).collapse(mesh)[1].values()]

Can you help with this problem?

Provide an answer of your own, or ask Mikael Mortensen for more information if necessary.

To post a message you must log in.