What's the best way to read in many fields from a VTU file?

I'm just shooting from the hip here, but it seems the best would be to
create read the VTU file and create a function based on it. I haven't
done this in years, but it use to be called a DiscreteFunction and you
would just assign vector values. Now, you can create your own mesh,
function space, and function then interpolate on top of your mesh.

Perhaps if you had a small example of what you are doing there would
be more to say.

-- Andy

On Tue, Oct 9, 2012 at 5:11 PM, Patrick Farrell
<email address hidden> wrote:
> New question #210788 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/210788
>
> Hi,
>
> Suppose I have P1 scalar fields stored in a VTU file, and I'd like to read these in to DOLFIN to do some computations.
>
> I can do something like the following (and I am):
>
> class MyField(Expression):
> def eval(self, values, x):
> ... do some VTK magic to interpolate the value of the VTU at the point x
>
> but this is extremely slow.
>
> What's the recommended way to read in a mesh and P1 fields from a VTU file?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

We don't have input from VTK. My answer would be not to save the
solution to VTK other than for visualization. Instead, store the dofs
in a TimeSeries, XML file or the new XDMF file format.

--
Anders

On Wed, Oct 10, 2012 at 02:31:05PM -0000, Andy R Terrel wrote:
> Question #210788 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210788
>
> Status: Open => Answered
>
> Andy R Terrel proposed the following answer:
> I'm just shooting from the hip here, but it seems the best would be to
> create read the VTU file and create a function based on it. I haven't
> done this in years, but it use to be called a DiscreteFunction and you
> would just assign vector values. Now, you can create your own mesh,
> function space, and function then interpolate on top of your mesh.
>
> Perhaps if you had a small example of what you are doing there would
> be more to say.
>
> -- Andy
>
> On Tue, Oct 9, 2012 at 5:11 PM, Patrick Farrell
> <email address hidden> wrote:
> > New question #210788 on DOLFIN:
> > https://answers.launchpad.net/dolfin/+question/210788
> >
> > Hi,
> >
> > Suppose I have P1 scalar fields stored in a VTU file, and I'd like to read these in to DOLFIN to do some computations.
> >
> > I can do something like the following (and I am):
> >
> > class MyField(Expression):
> > def eval(self, values, x):
> > ... do some VTK magic to interpolate the value of the VTU at the point x
> >
> > but this is extremely slow.
> >
> > What's the recommended way to read in a mesh and P1 fields from a VTU file?
> >
>
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Anders: the data I'm trying to read wasn't generated by dolfin, it is climatological data. I agree that VTK shouldn't be a format for the storage or transmission of data, but I had to work with what I had.

In any case, I worked around it: if you turn off parameters["reorder_dofs"], then you can read them in index-by-index (in serial), which is very fast. Then I can save them in a dolfin-friendly format, and read them in again in the main solver.

I have a script which more or less does this. It should be compatible with both vtu and pvd files. If you are interested:

  http://pastebin.com/UQhSxkzX

for the vtkread.py file and a copy of the tvtk file array_handler.py

  http://pastebin.com/4X7NAHwg

Then you should be able to do stuff like:

  from vtkread import VTKToDOLFIN

  # Should work with single vtu file too
  reader = VTKToDOLFIN("your_data.pvd", mesh=None) # Provide mesh if you already have it

  for i, u in enumerate(reader):
      if isinstance(u, tuple):
          t, u = u # Split out time if that is stored in the pvd file.
      # u is here a dolfin.Function

Not tested with latest trunk though...

Johan

On Wed, Oct 10, 2012 at 06:50:58PM -0000, Johan Hake wrote:
> Question #210788 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210788
>
> Johan Hake posted a new comment:
> I have a script which more or less does this. It should be compatible
> with both vtu and pvd files. If you are interested:
>
> http://pastebin.com/UQhSxkzX
>
> for the vtkread.py file and a copy of the tvtk file array_handler.py
>
> http://pastebin.com/4X7NAHwg
>
> Then you should be able to do stuff like:
>
> from vtkread import VTKToDOLFIN
>
> # Should work with single vtu file too
> reader = VTKToDOLFIN("your_data.pvd", mesh=None) # Provide mesh if you already have it
>
> for i, u in enumerate(reader):
> if isinstance(u, tuple):
> t, u = u # Split out time if that is stored in the pvd file.
> # u is here a dolfin.Function
>
> Not tested with latest trunk though...

This could go into the utils directory.

--
Anders

Could be. However, it turns out that the array_handler script I thought
was a standalone one, actually depends on enthought's tvtk. The latter
is a wrapper around vtk, making it easier for vtk->numpy data exchange.

One could probably clean it to not be dependent on tvtk. Not sure I am
up for the task though...

Johan

On 10/12/2012 05:01 AM, Anders Logg wrote:
> Question #210788 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210788
>
> Anders Logg posted a new comment:
> On Wed, Oct 10, 2012 at 06:50:58PM -0000, Johan Hake wrote:
>> Question #210788 on DOLFIN changed:
>> https://answers.launchpad.net/dolfin/+question/210788
>>
>> Johan Hake posted a new comment:
>> I have a script which more or less does this. It should be compatible
>> with both vtu and pvd files. If you are interested:
>>
>> http://pastebin.com/UQhSxkzX
>>
>> for the vtkread.py file and a copy of the tvtk file array_handler.py
>>
>> http://pastebin.com/4X7NAHwg
>>
>> Then you should be able to do stuff like:
>>
>> from vtkread import VTKToDOLFIN
>>
>> # Should work with single vtu file too
>> reader = VTKToDOLFIN("your_data.pvd", mesh=None) # Provide mesh if you already have it
>>
>> for i, u in enumerate(reader):
>> if isinstance(u, tuple):
>> t, u = u # Split out time if that is stored in the pvd file.
>> # u is here a dolfin.Function
>>
>> Not tested with latest trunk though...
>
> This could go into the utils directory.
>
> --
> Anders
>