Interpolation and PVD Animation

Overview

Most geometry so far is static: a model is built, meshed, and inspected at a single pose. This tutorial adds motion. Given a mesh and a final pose — a rotation matrix R and a translation vector t — Rasmah can interpolate the mesh from its current pose to that pose, producing a sequence of frames, and write them out as a ParaView Data (.pvd) animation that ParaView plays as a time series.

Two pieces:

• interpolate(mesh, R, t, n) — builds the n-frame sequence.
• write_pvd(meshes, path) — writes the sequence (one .vtp/.vtu per frame plus the .pvd collection).

Under the hood the rotation is interpolated with slerp (spherical linear interpolation), which we tour first.

Spherical linear interpolation

A rotation is a point on a sphere of orientations, and the correct “straight line” between two rotations is a spherical path — not a component-wise average of the matrices (which would not stay a rotation). Rasmah provides the matrix logarithm and exponential to build that path:

using Rasmah
using LinearAlgebra

R0 = Matrix(1.0I, 3, 3)            # identity
R1 = rotation_z(π / 2)             # quarter turn about z

slerp(R0, R1, 0.0)                 # R0
slerp(R0, R1, 0.5)                 # halfway → rotation_z(π/4)
slerp(R0, R1, 1.0)                 # R1

3×3 Matrix{Float64}:
 1.11022e-16  -1.0          0.0
 1.0           1.11022e-16  0.0
 0.0           0.0          1.0

slerp(R0, R1, s) = R0 · exp(s · log(R0ᵀ R1)), where rotation_log returns the skew-symmetric (axis-angle) matrix and rotation_exp reverses it via Rodrigues’ formula. The round-trip is exact:

R = rotation_x(0.7) * rotation_y(1.3)
println("exp(log(R)) == R: ", rotation_exp(rotation_log(R)) ≈ R)

exp(log(R)) == R: true

The mesh to animate

We animate a small bracket — a plate with a drilled hole — so the motion is easy to read:

model = difference(box(4.0, 2.0, 0.5), cylinder(0.5, 1.0))
mesh = surface_mesh(model; n=48)

println(size(mesh.vertices, 2), " vertices, ", size(mesh.faces, 2), " triangles")

vtk(mesh)

13464 vertices, 26928 triangles

Interpolating the pose

interpolate morphs the mesh from its current pose to the final pose v ↦ R v + t. The first frame is the original, the last is the full transform, and the frames in between are smoothly in between:

R = rotation_z(π / 2)
t = [2.0, 1.0, 0.0]

frames = interpolate(mesh, R, t, 9)          # 9 frames, kind=:rigid by default

println("frame 1 vertices == original: ", frames[1].vertices == mesh.vertices)
println("frame 9 vertices == R v + t: ",
        frames[9].vertices ≈ R * mesh.vertices .+ t)
println("frame 5 (mid) == rotation_z(π/4) v + t/2: ",
        frames[5].vertices ≈ rotation_z(π / 4) * mesh.vertices .+ 0.5 .* t)

vtk(frames[5])

frame 1 vertices == original: true
frame 9 vertices == R v + t: true
frame 5 (mid) == rotation_z(π/4) v + t/2: true

trans = interpolate(mesh, R, t, 5; kind=:translation)
rot   = interpolate(mesh, R, t, 5; kind=:rotation)

println("translation-only final == v + t: ",
        trans[5].vertices ≈ mesh.vertices .+ t)
println("rotation-only final == R v: ",
        rot[5].vertices ≈ R * mesh.vertices)

translation-only final == v + t: true
rotation-only final == R v: true

disp = [norm(rot[end].vertices[:, j] - mesh.vertices[:, j]) for j in axes(mesh.vertices, 2)]
vtk(rot[end]; field=disp, fieldname="displacement")

Easing

easing selects how the interpolation parameter s ∈ [0, 1] is sampled. :linear spaces the frames uniformly; :smooth (smoothstep, s²(3 − 2s)) eases in and out, so the motion starts and ends gently. A translation-only motion makes the pacing easy to read: both sequences end at the same offset, but the smooth one lags early on:

linear = interpolate(mesh, R, t, 7; kind=:translation, easing=:linear)
smooth = interpolate(mesh, R, t, 7; kind=:translation, easing=:smooth)

# frame 2 is one step in (s = 1/6): the x-offset reveals the pacing
println("linear frame 2 x-offset: ", round(linear[2].vertices[1, 1] - mesh.vertices[1, 1]; digits=3))
println("smooth frame 2 x-offset: ", round(smooth[2].vertices[1, 1] - mesh.vertices[1, 1]; digits=3))
println("final x-offset (both):  ", round(smooth[7].vertices[1, 1] - mesh.vertices[1, 1]; digits=3))

linear frame 2 x-offset: 0.333
smooth frame 2 x-offset: 0.148
final x-offset (both):  2.0

The uniform sequence is already at 2·(1/6) ≈ 0.333 by frame 2, while the smoothstep sequence is still easing in at 2·(1/6)²·(3−1/3) ≈ 0.148 — yet both reach the same 2.0 at the end.

Writing the PVD collection

write_pvd writes each frame to a .vtp (triangle) / .vtu (tetrahedron) file next to the .pvd, and the .pvd itself references them so ParaView plays them as a time series:

mktempdir() do dir
    path = joinpath(dir, "spin.pvd")
    write_pvd(frames, path)

    println("wrote:")
    for f in sort(readdir(dir))
        println("  ", f)
    end

    println("\n.pvd collection:")
    print(read(path, String))
end

wrote:
  spin.pvd
  spin_0000.vtp
  spin_0001.vtp
  spin_0002.vtp
  spin_0003.vtp
  spin_0004.vtp
  spin_0005.vtp
  spin_0006.vtp
  spin_0007.vtp
  spin_0008.vtp

.pvd collection:
<?xml version="1.0"?>
<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">
  <Collection>
    <DataSet timestep="0" group="" part="0" file="spin_0000.vtp"/>
    <DataSet timestep="1" group="" part="0" file="spin_0001.vtp"/>
    <DataSet timestep="2" group="" part="0" file="spin_0002.vtp"/>
    <DataSet timestep="3" group="" part="0" file="spin_0003.vtp"/>
    <DataSet timestep="4" group="" part="0" file="spin_0004.vtp"/>
    <DataSet timestep="5" group="" part="0" file="spin_0005.vtp"/>
    <DataSet timestep="6" group="" part="0" file="spin_0006.vtp"/>
    <DataSet timestep="7" group="" part="0" file="spin_0007.vtp"/>
    <DataSet timestep="8" group="" part="0" file="spin_0008.vtp"/>
  </Collection>
</VTKFile>

times overrides the default 0:n-1 timesteps, and format passes through to the VTK writer (:ascii, :binary, or :appended):

mktempdir() do dir
    write_pvd(frames, joinpath(dir, "spin.pvd"); times=collect(0.0:0.5:4.0))
    println("timesteps: ", join(0.0:0.5:4.0, ", "))
end

timesteps: 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0

Volume meshes work the same way — interpolate on a TetMesh writes a .vtu sequence:

vol = tetrahedralize_box(2.0, 1.0, 0.5, 6, 3, 2)
volframes = interpolate(vol, rotation_z(π / 4), [0.0, 0.0, 1.0], 4)

mktempdir() do dir
    write_pvd(volframes, joinpath(dir, "vol.pvd"))
    println("volume frames: ", filter(f -> endswith(f, ".vtu"), readdir(dir)))
end

volume frames: ["vol_0000.vtu", "vol_0001.vtu", "vol_0002.vtu", "vol_0003.vtu"]

Open the .pvd in ParaView to play the animation.

Reading the PVD manifest back

read_pvd(path) is the complement of write_pvd: it parses the .pvd collection’s <DataSet> list and returns the frame directory plus the ordered (file, timestep) entries — the same information the animated viewer uses to load a .pvd written by ParaView or by write_pvd:

mktempdir() do dir
    path = joinpath(dir, "spin.pvd")
    write_pvd(frames, path; times=collect(0.0:0.5:4.0))

    fdir, entries = read_pvd(path)
    println("frame directory: ", fdir)
    for (file, t) in entries
        println("  t = ", t, "  →  ", file)
    end
end

frame directory: /tmp/jl_T39JQM
  t = 0.0  →  spin_0000.vtp
  t = 0.5  →  spin_0001.vtp
  t = 1.0  →  spin_0002.vtp
  t = 1.5  →  spin_0003.vtp
  t = 2.0  →  spin_0004.vtp
  t = 2.5  →  spin_0005.vtp
  t = 3.0  →  spin_0006.vtp
  t = 3.5  →  spin_0007.vtp
  t = 4.0  →  spin_0008.vtp

Next steps

• Differentiable Surface Meshing — the frozen-topology meshes these transforms move.
• VTK: Legacy and XML — the underlying VTK writers write_pvd builds on.