FilamentIO

FilamentIO

Module defining tools for reading and writing filament data.

write_vtkhdf(
    [func::Function],
    filename::AbstractString,
    fs::AbstractVector{<:AbstractFilament};
    refinement = 1,
    dataset_type = :PolyData,
    parametrisation = true,
    periods = (nothing, nothing, nothing),
)

Write new VTK HDF file containing a list of filaments.

A VTK HDF file is an HDF5 file organised in such a way that it can be readily visualised in tools such as ParaView. Data in the file can be readily accessed using HDF5 tools and libraries.

Following the VTK HDF specification, this function creates a "VTKHDF" group on top of the HDF5 file. Then, it creates the datasets allowing to describe the filaments as an unstructured grid.

See also the VTK documentation for details on the VTK HDF format.

Extended help

Relevant datasets

Some relevant datasets which are written are:

• /VTKHDF/Points: contains the coordinates of all filament nodes. Points are represented as an array of dimensions (3, Np) where Np is the total number of nodes. Note that the points include the endpoint, meaning that for a closed filament the initial coordinate appears twice in the file. This is done to disambiguate between closed and infinite but unclosed filaments (see end_to_end_offset).

• /VTKHDF/CellData/FilamentIds: contains the filament id associated to each VTK cell. It is a dataset of length Nc, where Nc is the number of cells. In our case a cell is a spatial curve made of discrete points. The ids are Int32 values from 1 to Nf, where Nf is the number of filaments. The values are always sorted increasingly (FilamentIds[i + 1] ≥ FilamentIds[i] for all i). Note that, in general, one cell corresponds to one filament (so that Nc == Nf), but this may not be the case when the periods argument is used and a filament is broken onto multiple curves (see Periodic wrapping of filaments below for details).

Attaching extra data

The optional func argument can be used to attach other data (such as velocity vectors or the current time) to the generated file. This is most conveniently done using the do block syntax. See further below for some examples.

Optional keyword arguments

• refinement::Int = 1: allows to output more than 1 point for each filament segment. This is mostly useful for producing nice visualisations. The level of refinement is written to the /VTKHDF/RefinementLevel dataset, which allows to read back the data skipping intra-segment nodes.

• parametrisation::Bool = true: if true (default), write curve parametrisation values. This allows read_vtkhdf to reconstruct the exact same curve that was written (assuming the same discretisation method is used, e.g. cubic splines), even if the curve used some non-default curve parametrisation. Values are written to the /VTKHDF/PointData/Parametrisation dataset.

• dataset_type::Symbol: can be either :PolyData (default) or :UnstructuredGrid. There's usually no reason to change this.

Periodic wrapping of filaments

When using periodic boundary conditions, one can use the optional periods argument to periodically wrap the filaments. In this case, one should pass a tuple periods = (Lx, Ly, Lz) with the period in each direction (one can pass nothing if there are non-periodic directions).

In this case, filaments outside the main periodic box will be translated to fit in the periodic box. Moreover, if a filament locally goes out of the periodic box, it will be broken onto multiple curves so that they all fit within the domain. One can then look at the /VTKHDF/CellData/FilamentIds dataset to determine which curves belong to the same filament.

Typical usage

# Note: the file extension is arbitrary, but ParaView prefers ".vtkhdf" (or ".hdf") if one
# wants to use the files for visualisation.
write_vtkhdf("filaments.vtkhdf", fs; refinement = 2, periods = (2π, 2π, 2π)) do io
    io["Velocity"] = vs  # adds velocity as VTK point data, assuming vs is a VectorOfVectors
    io["Curvatures"] = CurvatureVector()  # convenient syntax for writing geometric quantities along filaments
    io["Time"] = 0.3     # adds time as VTK field data, since it's a scalar
    # one can add other fields here...
end

Base.setindex!(io::VTKHDFFile, vs::AbstractVector{<:AbstractVector}, name::AbstractString)

Attach data to filament nodes.

One generally wants to use the syntax io[name] = vs which calls this function.

This can be used to write fields defined at filament nodes (for instance, the velocity of each node).

The data is written to the dataset /VTKHDF/PointData/$name.

For vector fields (such as velocity), the written dataset has dimensions (3, Np) where Np is the total number of filament nodes (including endpoints). The format is exactly the same as for the Points dataset as detailed in write_vtkhdf. As also explained there, the Offsets dataset can be used to recover the values associated to each filament.

Writing geometric quantities

It is also possible to write geometric quantities such as unit tangents or curvature vectors along filaments (see GeometricQuantity for a list). For this, one can pass the wanted quantity as the vs argument.

For example, one can simply do

io["Curvatures"] = CurvatureVector()

to write curvature vectors along each filament.

Base.setindex!(io::VTKHDFFile, data, name::AbstractString)

Write data as VTK field data to VTK HDF file.

In VTK, field data refers to data which is not directly attached to the geometry. This is typically small datasets or simple values, such as the current time or simulation parameters.

Note that scalar data (such as time) is always written as a single-element vector, since otherwise it cannot be parsed by VTK.

This function interprets everything that is not a vector of vectors as field data.

read_vtkhdf(
    [func::Function],
    filename::AbstractString,
    ::Type{T},
    method::DiscretisationMethod,
) where {T}

Read filament locations from VTK HDF file.

This function loads filaments based on the datasets /VTKHDF/Points and /VTKHDF/Offsets as written by the write_vtkhdf function.

Returns a vector of filaments, where each filament is discretised according to the chosen method. See Filaments.init for possible options.

One can also read other datasets using read and read!, as shown and explained below.

Extended help

Typical usage

local vs, t  # make sure these variables still exist after the `do` block

# The returned `fs` is a list of filaments.
fs = read_vtkhdf("filaments.vtkhdf", Float64, CubicSplineMethod()) do io
    vs = read(io, "Velocity", PointData())            # here `vs` contains one velocity vector per filament node
    t = only(read(io, "Time", FieldData(), Float64))  # note: field data is always written as an array
    # one can read other fields here...
end

The available reading functions are:

• read(io, name::AbstractString, ::PointData) for reading point data (e.g. a velocity field);

• read(io, name::AbstractString, ::FieldData, ::Type{T}) for reading field data (i.e. data not attached to filament nodes, such as the current time);

• read!(io, vs::AbstractVector{<:AbstractVector}, name::AbstractString) for reading point data onto a preallocated vector of vectors.

Accessing filament data in do block

When using the do block syntax as in the above example, one may want to have access to the filament locations fs from within the do block. For instance, this can be useful if one has a preallocated vector of velocities vs which needs to be resized to match the number of filaments and filament nodes, before reading values using read!.

In fact, fs can be easily obtained from the io object:

read_vtkhdf("filaments.vtkhdf", Float64, CubicSplineMethod()) do io
    fs = io.fs  # this is the vector of filaments
    # Say we want to resize an existent vector of velocities (must have the right type...):
    resize!(vs, length(fs))
    for (v, f) ∈ zip(vs, fs)
        resize!(v, length(nodes(f)))  # resize velocities of a single filament
    end
    # Now we can read the velocities
    read!(io, vs, "Velocity")
end

TimeSeriesFile() -> TimeSeriesFile

Initialise a time series file compatible with ParaView.

A ParaView time series file is a JSON which points to a series of VTK files representing a time series. To each VTK file one can associate a (simulation) time.

One can use this to describe a series of VTKHDF files generated by write_vtkhdf. Note that, in this case, the extension of the time series file must be .vtkhdf.series.

Typical usage

# Initialise TimeSeriesFile before starting a simulation
tsf = TimeSeriesFile()

# Write a VTKHDF file associated to a single time
timestep = 4200
time = 24.0
filename = "rings_$(timestep).vtkhdf"
write_vtkhdf(filename, filaments, etc...)  # write VTKHDF file
tsf[time] = filename                       # add VTKHDF file to time series file

# Save the file at the end of the simulation (or at an intermediate time, that's OK too)
save("rings.vtkhdf.series", tsf)  # note: the extension must be .vtkhdf.series to make ParaView happy

Base.setindex!(tsf::TimeSeriesFile, filename::AbstractString, time::Real)

Add file to time series file.

Typical usage

timestep = 4200
time = 24.0
filename = "rings_$(timestep).vtkhdf"
tsf[time] = filename

Base.empty!(tsf::TimeSeriesFile)

Reset TimeSeriesFile, removing all entries.

save(filename::AbstractString, tsf::TimeSeriesFile)
save(io::IO, tsf::TimeSeriesFile)

Write TimeSeriesFile to a file.

The filename must have extension .VTKFORMAT.series. In particular, if writing VTKHDF files, the format must be .vtkhdf.series.

write_to_text(io::IO, fs::AbstractVector{<:AbstractVector})
write_to_text(filename::AbstractString, fs::AbstractVector{<:AbstractVector})

Dump filament locations or other filament data to text file.

Here fs can be a list of filaments or a list of values (e.g. velocities) on filament nodes.

The output is a text file where:

• each line represents a single filament,

• node coordinates (or vector values such as velocity) are flattened so that they can be written in a single line,

• the end node is included. For a closed filament, this means that the first and last written points are identical. This also means that we write N + 1 values for a filament of length N (= number of independent nodes).

One can use read_from_text to generate filaments from files written by this function.

read_from_text(io::IO, ::Type{T}, method::DiscretisationMethod) -> fs
read_from_text(filename::AbstractString, ::Type{T}, method::DiscretisationMethod) -> fs

Read filament locations from text file.

This function can be used to read filament locations written by write_to_text.

The type T corresponds to the wanted numerical precision. It is usually Float64 or Float32.