Filaments

Filaments

Module for dealing with the discretisation of curves in 3D space.

AbstractFilament{T} <: AbstractVector{Vec3{T}}

Abstract type representing a curve in 3D space.

The curve coordinates are parametrised as $\bm{s}(t)$ with $t ∈ ℝ$.

The curve is discretised by a set of nodes (or discretisation points) $\bm{s}(t_i) = \bm{s}_i$ for $i ∈ \{1, 2, …, N\}$.

See ClosedFilament for a concrete implementation of AbstractFilament.

An AbstractFilament is treated as an AbstractVector of length N, in which each element is a discretisation point $\bm{s}_i$. Therefore, one can use the usual indexing notation to retrieve and to modify discretisation points. See ClosedFilament for some examples.

Extended help

Evaluating coordinates and derivatives

Different ways are proposed of evaluating filament coordinates and derivatives along a filament f, depending on whether one wants values on discretisation points or in-between them.

In short, square brackets f[...] should be used to evaluate on filament nodes, while round brackets f(...) to evaluate in-between nodes.

Values on discretisation points

Coordinates $\bm{s}_i$ of discretisation points can be simply obtained by indexing the filament object:

s⃗ = f[i]

Derivatives at discretisation points can be similarly obtained by doing:

s⃗′ = f[i, Derivative(1)]
s⃗″ = f[i, Derivative(2)]

(Note that this also works with Derivative(0), in which case it's the same as f[i].)

For convenience, other geometric quantities can be evaluated in a similar way:

ρ⃗ = f[i, CurvatureVector()]

(See GeometricQuantity for available quantities.)

Values in-between discretisation points

In this case, one wants to evaluate a value in-between two discretisation points, i.e. for $t ∈ [t_i, t_{i + 1}]$ for some index i. Two options are proposed:

• if one knows the parametrisation of the filament (see knots to obtain the parametrisation knots), then one can evaluate using a value of t:

  s⃗  = f(t)
  s⃗′ = f(t, Derivative(1))
  s⃗″ = f(t, Derivative(2))

• alternatively, if one knows the index i associated to the segment of interest, then one can do

  s⃗  = f(i, ζ)
  s⃗′ = f(i, ζ, Derivative(1))
  s⃗″ = f(i, ζ, Derivative(2))

  where ζ must be in $[0, 1]$, and the two limits correspond to knots $t_i$ and $t_{i + 1}$. This is convenient if one wants to evaluate, say, right in the middle between two discretisation points, in which case one would choose ζ = 0.5.

For convenience, other geometric quantities can be evaluated in a similar way:

ρ⃗ = f(t, CurvatureVector())
ρ⃗ = f(i, ζ, CurvatureVector())

ClosedFilament{T, D <: DiscretisationMethod} <: AbstractFilament{T}

Describes a closed curve (a loop) in 3D space.

It can also be used to represent infinite but unclosed curves described by a periodic function (such as an infinite straight line or a sinusoidal curve).

ClosedFilaments should be generally constructed using Filaments.init.

Extended help

Examples

The following examples use the CubicSplineMethod for representing filament curves, but other methods are also available.

Initialise filament from a set of discretisation points:

julia> f = Filaments.init(ClosedFilament, 16, CubicSplineMethod());

julia> θs = range(-1, 1; length = 17)[1:16]
-1.0:0.125:0.875

julia> @. f = Vec3(cospi(θs), sinpi(θs), 0);

julia> f[4]
3-element SVector{3, Float64} with indices SOneTo(3):
 -0.3826834323650898
 -0.9238795325112867
  0.0

julia> f[5] = (f[4] + 2 * f[6]) ./ 2
3-element SVector{3, Float64} with indices SOneTo(3):
  0.1913417161825449
 -1.38581929876693
  0.0

julia> update_coefficients!(f);

Note that update_coefficients! should be called whenever filament coordinates are changed, before doing other operations such as estimating derivatives.

Estimate derivatives at discretisation points:

julia> f[4, Derivative(1)]
3-element SVector{3, Float64} with indices SOneTo(3):
  0.9090457394297018
 -0.7273334611006509
  0.0

julia> f[4, Derivative(2)]
3-element SVector{3, Float64} with indices SOneTo(3):
  0.20911715113294102
 -2.09047051482799
  0.0

Estimate coordinates and derivatives in-between discretisation points:

julia> f(4, 0.32)
3-element SVector{3, Float64} with indices SOneTo(3):
 -0.16753415613203387
 -1.1324592487590195
  0.0

julia> Ẋ, Ẍ = f(4, 0.32, Derivative(1)), f(4, 0.32, Derivative(2))
([0.8947546127964856, -0.9527970723463657, 0.0], [-0.3303413370703831, 0.17798009799460934, 0.0])

julia> X′, X″ = f(4, 0.32, UnitTangent()), f(4, 0.32, CurvatureVector())
([0.6845546705034081, -0.7289615237390588, 0.0], [-0.050762951240829336, -0.047670575508846375, 0.0])

Curve parametrisation

The parametrisation knots $t_i$ are directly obtained from the interpolation point positions. A standard choice, which is used here by default, is for the knot increments to approximate the arc length between two interpolation points:

\[ℓ_{i} ≡ t_{i + 1} - t_{i} = |\bm{s}_{i + 1} - \bm{s}_i|,\]

which is a zero-th order approximation (and a lower bound) for the actual arc length between points $\bm{s}_i$ and $\bm{s}_{i + 1}$.

Vec3{T}

Three-element static vector, alias for SVector{3, T}.

Used to describe vectors and coordinates in 3D space.

Derivative{N}

Represents the $N$-th order derivative operator.

Used in particular to interpolate derivatives along filaments.

Filaments.init(
    ClosedFilament{T}, N::Integer, method::DiscretisationMethod;
    offset = zero(Vec3{T}),
    parametrisation = Filaments.default_parametrisation(method),
    nderivs = Val(2),
) -> ClosedFilament{T}

Allocate data for a closed filament with N discretisation points.

The element type T can be omitted, in which case the default T = Float64 is used.

The optional offset keyword argument allows to define a filament with a spatial offset between points f[i] and f[i + N]. By default the offset is zero, meaning that the filament is a closed loop. This can be used for defining infinite (so not really closed) filaments living in periodic domains.

Possible discretisation methods include CubicSplineMethod, QuinticSplineMethod and FiniteDiffMethod.

Filaments.init(
    ClosedFilament, points::AbstractVector{<:Vec3}, method::DiscretisationMethod;
    kwargs...,
)

Initialise new filament with the chosen discretisation points.

Note that update_coefficients! does not need to be called after using this variant (until node locations change, of course!).

Filaments.init(S::Function, ClosedFilament{T}, τs::AbstractVector, method::DiscretisationMethod)
Filaments.init(S::Function, ClosedFilament{T}, N::Integer, method::DiscretisationMethod)

Construct new filament from parametric function $S(τ): [0, 1] ↦ ℝ³$.

The parametric function $S(τ)$ must generate a 3D coordinate (for example as an SVector{3} or an NTuple{3}) for any $τ ∈ [0, 1]$. The curve described by $S$ is closed if and only if $S(0) = S(1)$. More generally, for infinite but unclosed curves, the end-to-end offset is obtained as $Δ = S(1) - S(0)$.

In particular, one can pass functions generated by PredefinedCurves.define_curve.

Two variants of this method are provided:

• in the first variant, the function $S$ is evaluated at user-provided values of $τ$, which should be in $[0, 1)$;

• in the second variant, one simply passes the desired number $N$ of nodes of the resulting filament, and the function chooses an equispaced list of $τ$ values.

The element type T can be omitted, in which case the default T = Float64 is used.

Some examples are given further below ("Initialising filament from parametric function"). See also from_vector_field to initialise a filament from an analytical vector field.

Extended help

Evaluating high-order derivatives

By default the filament discretisation allows to evaluate up to the second derivative of the curve with respect to its parametrisation. This means in particular than one has access to the tangent, curvature and binormal vectors.

One may want to evaluate higher order derivatives, for example to also have access to the curve torsion (which depends on the 3rd derivative). This is only possible for some high-order discretisation methods (QuinticSplineMethod, as well as FourierMethod with the limitation that one has to evaluate on discretisation points). Moreover, one needs to explicitly pass nderivs = Val(3) (or higher) to the init function.

Customising the filament parametrisation

By default the filament is parametrised as $\bm{s}(t)$ such that the parameter $t$ roughly corresponds to the arc length $ξ$.^[1] More precisely, the discrete values $t_i$ at the discretisation points are defined from the distances between discretisation points:

\[t_{i + 1} = t_{i} + |\bm{s}_{i + 1} - \bm{s}_i|, \quad t_1 = 0.\]

One can use the parametrisation keyword argument to change this. This may be a function parametrisation(Xs, i) which takes the vector Xs of discretisation points and the index i associated to the start of the current segment.

For instance, to parametrise filaments according to the $z$ increments between discretisation points, one can pass:

parametrisation = (Xs, i) -> abs(Xs[i + 1].z - Xs[i].z)

One can also pass a predefined parametrisation. See AbstractParametrisation for a list of options.

Examples

Initialising filament from parametric function

Initialise circular filament with N = 16 nodes:

julia> S(t) = (cos(2π * t), sin(2π * t), 0)  # define a circular ring with period T = 1
S (generic function with 1 method)

julia> N = 16;

julia> f = Filaments.init(S, ClosedFilament, N, CubicSplineMethod())
16-element ClosedFilament{SVector{3, Float64}, CubicSplineMethod}:
 [0.9807852804032304, 0.19509032201612825, 0.0]
 [0.8314696123025452, 0.5555702330196022, 0.0]
 [0.5555702330196023, 0.8314696123025452, 0.0]
 [0.19509032201612833, 0.9807852804032304, 0.0]
 [-0.1950903220161282, 0.9807852804032304, 0.0]
 [-0.555570233019602, 0.8314696123025453, 0.0]
 [-0.8314696123025453, 0.5555702330196022, 0.0]
 [-0.9807852804032304, 0.1950903220161286, 0.0]
 [-0.9807852804032304, -0.19509032201612836, 0.0]
 [-0.8314696123025455, -0.555570233019602, 0.0]
 [-0.5555702330196022, -0.8314696123025452, 0.0]
 [-0.19509032201612866, -0.9807852804032303, 0.0]
 [0.1950903220161283, -0.9807852804032304, 0.0]
 [0.5555702330196018, -0.8314696123025455, 0.0]
 [0.8314696123025452, -0.5555702330196022, 0.0]
 [0.9807852804032303, -0.19509032201612872, 0.0]

Same but choosing the locations τ:

julia> τs = range(0, 1; length = N + 1)[1:N]  # make sure the location τ = 1 is *not* included!
0.0:0.0625:0.9375

julia> f = Filaments.init(S, ClosedFilament, τs, CubicSplineMethod())
16-element ClosedFilament{SVector{3, Float64}, CubicSplineMethod}:
 [1.0, 0.0, 0.0]
 [0.9238795325112867, 0.3826834323650898, 0.0]
 [0.7071067811865476, 0.7071067811865475, 0.0]
 [0.38268343236508984, 0.9238795325112867, 0.0]
 [6.123233995736766e-17, 1.0, 0.0]
 [-0.3826834323650897, 0.9238795325112867, 0.0]
 [-0.7071067811865475, 0.7071067811865476, 0.0]
 [-0.9238795325112867, 0.3826834323650899, 0.0]
 [-1.0, 1.2246467991473532e-16, 0.0]
 [-0.9238795325112868, -0.38268343236508967, 0.0]
 [-0.7071067811865477, -0.7071067811865475, 0.0]
 [-0.38268343236509034, -0.9238795325112865, 0.0]
 [-1.8369701987210297e-16, -1.0, 0.0]
 [0.38268343236509, -0.9238795325112866, 0.0]
 [0.7071067811865474, -0.7071067811865477, 0.0]
 [0.9238795325112865, -0.3826834323650904, 0.0]

Using VortexPasta.PredefinedCurves:

julia> using VortexPasta.PredefinedCurves

julia> trefoil = define_curve(TrefoilKnot());

julia> f = Filaments.init(trefoil, ClosedFilament, N, CubicSplineMethod())
16-element ClosedFilament{SVector{3, Float64}, CubicSplineMethod}:
 [0.9604571867463079, -0.866973784619343, -0.5555702330196022]
 [2.4033292980421757, 0.06610274757236567, -0.9807852804032304]
 [2.679228677325119, 1.3209370977497819, -0.19509032201612828]
 [1.74615214513341, 2.042849387038702, 0.8314696123025452]
 [0.21541841567305087, 1.6526687430064453, 0.8314696123025452]
 [-1.0162894527200281, 0.20979663171057739, -0.19509032201612828]
 [-1.2921888320029713, -1.5968364770327248, -0.9807852804032304]
 [-0.5702765427140513, -2.828544345425804, -0.5555702330196022]
 [0.5702765427140513, -2.828544345425804, 0.5555702330196022]
 [1.2921888320029713, -1.5968364770327248, 0.9807852804032304]
 [1.0162894527200281, 0.20979663171057739, 0.19509032201612828]
 [-0.21541841567305087, 1.6526687430064453, -0.8314696123025452]
 [-1.74615214513341, 2.042849387038702, -0.8314696123025452]
 [-2.679228677325119, 1.3209370977497819, 0.19509032201612828]
 [-2.4033292980421757, 0.06610274757236567, 0.9807852804032304]
 [-0.9604571867463079, -0.866973784619343, 0.5555702330196022]

Filaments.from_vector_field(
    ClosedFilament, ωf::Function, s⃗₀, dτ, method::DiscretisationMethod;
    max_steps = 1000, nsubsteps = 1, redistribute = true,
) -> ClosedFilament

Initialise closed filament from vector field.

Here ωf is a function $\bm{ω}(\bm{x})$ which takes a 3D vector x⃗ and returns a vector value ω⃗. A filament will be created such that, for each point $\bm{s}$ on the filament, the tangent vector $\bm{s}'$ is parallel to the vector field at that point, $\bm{ω}(\bm{s})$. The field should be such that lines constructed this way are closed.

A possible application of this function is for constructing a filament which approximates vortex lines from a vorticity field $\bm{ω}(\bm{x})$, which are precisely defined in this way. In that case, if one actually knows the velocity field $\bm{v}(\bm{x})$ – such that $\bm{ω} = \bm{\nabla} × \bm{v}$ – and one is too lazy to analytically derive the corresponding vorticity field, one can pass ωf = Filaments.curl(vf) (see curl) where vf is a function defining the velocity field.

One may use Filaments.distance_to_field to verify the result of this function.

Positional arguments

• ωf::Function: function taking a 3D coordinate x⃗ and returning a vector value ω⃗;

• s⃗₀::Vec3: location where to start iterating. This point is guaranteed to be in the generated filament;

• dτ::Real: approximately distance between filament nodes. Determines the number of nodes in the resulting filament. The value of dτ/nsubsteps should be small enough to converge and so that the discretised lines are properly closed (see below for details);

• method::DiscretisationMethod: discretisation method to use (e.g. QuinticSplineMethod()).

Optional keyword arguments

• max_steps::Int = 1000: maximum number of steps before we stop iterating. This is also the maximum possible length of the returned filament;

• nsubsteps::Int = 1: number of solver substeps to perform for each spatial increment dτ. Larger values may be used to improve accuracy;

• redistribute = true: if true (default), redistribute_nodes! is called at the end to make sure that nodes are approximately distributed in a uniform way along the filament.

Extended help

Example

Generate a single filament aligned with the Taylor–Green vorticity field.

For convenience, we work with the Taylor–Green velocity field, and the vorticity is obtained via automatic differentiation (but we could directly work with the analytical vorticity instead).

julia> function taylor_green_velocity(x⃗::Vec3)
           x, y, z = x⃗
           Vec3(
               cos(x) * sin(y) * cos(z),
               -sin(x) * cos(y) * cos(z),
               0,
           )
       end
taylor_green_velocity (generic function with 1 method)

julia> ωf = Filaments.curl(taylor_green_velocity)  # vorticity field (via automatic differentiation)
#60 (generic function with 1 method)

julia> s⃗₀ = Vec3(0.1, 1.8, 0.42);  # starting point for creating the filament

julia> dτ = 0.1;  # pseudo time-step (has units of length; determines line resolution)

julia> nsubsteps = 4;  # this is just to ensure convergence (should be larger for larger dτ)

julia> f = Filaments.from_vector_field(ClosedFilament, ωf, s⃗₀, dτ, QuinticSplineMethod(); nsubsteps);

julia> summary(f)
"29-element ClosedFilament{SVector{3, Float64}, QuinticSplineMethod}"

julia> Filaments.distance_to_field(ωf, f)  # check that we're close to the actual vortex line
3.9561046185765664e-5

Implementation details

The filament is generated by numerically solving the ODE:

\[\frac{\mathrm{d}\bm{s}(τ)}{\mathrm{d}τ} = \hat{\bm{ω}}(\bm{s}), \quad \bm{s}(0) = \bm{s}_0,\]

where $τ$ denotes a "pseudo-time" (which actually has units of a length) and $\hat{\bm{ω}} = \bm{ω} / |\bm{ω}|$ is a unitary vector aligned with the vector field. The ODE is solved numerically using a standard 4th order Runge–Kutta scheme.

In this context, the $dτ$ argument is actually the "timestep" used when solving this ODE. It must be small enough so that the curve is accurately tracked.

Note that the curve will be automatically closed (and the ODE stopped) if we reach an $\bm{s}(τ)$ which is sufficiently close (closer than $dτ/2$) to the starting point $\bm{s}_0$. If that never happens, we stop after we have performed max_steps solver iterations.

DiscretisationMethod

Abstract type defining a filament discretisation method.

SplineMethod{k} <: DiscretisationMethod

Represents parametric curves using splines of order $k$.

A spline of order $k$ is a piecewise polynomial with polynomial degree $k - 1$ in-between interpolation nodes. For instance, cubic splines correspond to $k = 4$. On an interpolation node, the function has continuity $C^{k - 2}$.

See also CubicSplineMethod and QuinticSplineMethod.

CubicSplineMethod <: DiscretisationMethod

Represents curves using cubic splines.

This is an alias for SplineMethod{4} (see SplineMethod).

In the case of closed curves, periodic cubic splines are used.

QuinticSplineMethod <: DiscretisationMethod

Represents curves using quintic splines.

This is an alias for SplineMethod{6} (see SplineMethod).

A quintic spline is made of polynomials of degree 5 and has global continuity $C^4$.

In the case of closed curves, periodic quintic splines are used.

FiniteDiffMethod{M} <: DiscretisationMethod{M}
FiniteDiffMethod([M = 2], [interpolation = HermiteInterpolation(M)])

Estimation of curve derivatives at filament nodes using finite differences.

For now, only the case M = 2 (default) is implemented (4th order / 5-point finite differences), following the method used by Baggaley and Barenghi (2011) based on a paper by Gamet et al. (1999).

FiniteDiffMethod also requires specifying an interpolation scheme for evaluating coordinates and derivatives in-between discretisation points. By default, Hermite interpolations of continuity $C^M$ are used. For the case M = 2, that means quintic Hermite interpolations, which match the first two derivatives estimated by finite differences at the discretisation points.

References

• (Baggaley and Barenghi, 2011) Baggaley & Barenghi, Phys. Rev. B 83, 134509 (2011)
• (Gamet et al., 1999) Gamet et al., Int. J. Numer. Meth. Fluids 29, 2 (1999)

FourierMethod <: DiscretisationMethod
FourierMethod(interpolation = HermiteInterpolation(2))

Represents closed curves using Fourier series.

Derivatives at nodes are estimated using FFTs. To interpolate in-between nodes, a local interpolation method (typically HermiteInterpolation) is used, which is much faster (but less accurate) than evaluating Fourier series.

Note that using FFTs require the knot locations $t_i$ to be equidistant. The default parametrisation used by this method ensures this. However, this usually works best when the distance between discretisation points is more or less constant.

This method should only be used for simple settings and for verification of filament derivatives in other methods.

discretisation_method(f::AbstractFilament) -> DiscretisationMethod

Return the method used to discretise the filament based on its node locations.

interpolation_method(m::DiscretisationMethod)

Return the interpolation method associated to a DiscretisationMethod.

For FiniteDiffMethod and FourierMethod, this is usually a HermiteInterpolation (which requires derivatives at the interpolation nodes). For SplineMethods, since splines don't rely on a separate interpolation scheme, this simply returns m (i.e. the passed SplineMethod).

required_derivatives(m::DiscretisationMethod) -> Int

Return the number of derivatives on interpolation nodes required by the discretisation method.

This is generally larger than 0 for methods relying on Hermite interpolations (which require derivatives on interpolation nodes). This is 0 for other methods such as SplineMethod.

init_coefficients(method::DiscretisationMethod, ys::AbstractVector, [nderivs::Val]) -> DiscretisationCoefs
init_coefficients(method::DiscretisationMethod, cs::PaddedVector, cderivs::NTuple)  -> DiscretisationCoefs

Initialise interpolation coefficients.

In the first case, the ys vector is interpreted as a vector of values at interpolation nodes, and it will not be modified by the function. Note that this requires the allocation of a coefficient vector. Moreover, the nderivs argument indicates the number of derivatives that one wants to be able to compute. By default, nderiv is chosen as the minimum number of derivatives required by the method (see required_derivatives).

In the second case, cs is interpreted as a vector of interpolation coefficients, and is thus stored in the returned structure. Moreover, to compute one or more derivatives, one can pass cderivs vectors which should have the same type and length as cs.

compute_coefficients!(coefs::DiscretisationCoefs, [ys::AbstractVector], ts::PaddedVector)

Compute interpolation coefficients of parametric function $f(t)$.

Here ys contains values at interpolation nodes, and ts contains the values of the curve parameter $t$ at the nodes.

If ys is not passed, then the coefs.cs vector is expected to have the values at interpolation points (which will be overwritten with the interpolation coefficients).

This should be called before evaluating an interpolation using evaluate.

evaluate(coefs::DiscretisationCoefs, ts::PaddedVector, t::Real, [Derivative(0)]; [ileft = nothing])
evaluate(coefs::DiscretisationCoefs, ts::PaddedVector, i::Int, ζ::Real, [Derivative(0)])

Evaluate interpolation at a given location $t$ or $ζ$.

In the first case, $t$ corresponds to the global curve parametrisation. One can optionally indicate the segment index i (such that ts[i] ≤ t < ts[i + 1]) using the ileft keyword argument.

In the second case, $ζ ∈ [0, 1]$ corresponds to the local curve parametrisation, within the segment i.

LocalInterpolationMethod

Abstract type defining a local interpolation method for estimating curve properties (coordinates, derivatives) in-between interpolation points using local information.

HermiteInterpolation{M} <: LocalInterpolationMethod

Hermite interpolation of continuity $C^M$ at interpolation points.

Hermite interpolations are obtained using curve derivatives up to order $M$.

Allowed cases

• for $M = 0$ this is simply linear interpolation (note that curvatures cannot be estimated from linear interpolations);

• for $M = 1$ this is the standard Hermite interpolation (piecewise cubic polynomials, requires first derivatives at interpolation points);

• for $M = 2$ this is a quintic Hermite interpolation requiring first and second derivatives at interpolation points.

interpolate(
    method::LocalInterpolationMethod, derivative::Derivative, t::Number,
    values_at_interpolation_points...,
)

Interpolate coordinates $\bm{X}(t)$ or their derivatives at a given location $t$.

Under this form, the location $t$ must be normalised such that the interval of interest is in $t ∈ [0, 1]$. Note that input and output derivatives must also be normalised accordingly.

The values_at_interpolation_points depend on the interpolation method:

• for HermiteInterpolation, one should pass the coordinates $\bm{X}$ and the first $M$ derivatives at the two endpoints of the interval ($t = 0$ and $1$).

knots(f::AbstractFilament{T}) -> AbstractVector{T}

Return parametrisation knots $t_i$ of the filament.

Filaments are parametrised by $\bm{s}(t)$ for $t ∈ [0, T]$.

knotlims(f::AbstractFilament) -> (t_begin, t_end)

Return limits within which the filament can be evaluated.

end_to_end_offset(f::ClosedFilament{T}) -> Vec3{T}

Return the end-to-end offset Δ⃗ = f[end + 1] - f[begin] of a "closed" filament.

For actually closed filaments, the end-to-end offset is zero. However, ClosedFilament also supports the case of infinite (but unclosed) filaments, which infinitely extend along one or more Cartesian directions. The restriction imposed by ClosedFilament is that infinite filaments repeat themselves periodically, such that f[i + m * N] == f[i] + m * Δ⃗ where N is the length of the filament (i.e. the number of degrees of freedom, or the total number of independent filament nodes).

Base.minimum(node_distance, f::AbstractFilament) -> Real
Base.minimum(node_distance, fs::AbstractVector{<:AbstractFilament}) -> Real

Return the minimum distance $δ$ between neighbouring discretisation points of a filament.

If a vector of filaments is passed (variant 2), then the minimum among all filaments is returned.

Base.minimum(knot_increment, f::AbstractFilament) -> Real
Base.minimum(knot_increment, fs::AbstractVector{<:AbstractFilament}) -> Real

Return the minimum increment $Δt = t_{i + 1} - t_{i}$ between filament knots.

This is generally a good approximation for the minimum segment length, at least when the default ChordalParametrisation is used. If this is not the case, it is better to pass node_distance instead of knot_increment, which always returns the distance between two neighbouring nodes.

Base.maximum(node_distance, f::AbstractFilament) -> Real
Base.maximum(node_distance, fs::AbstractVector{<:AbstractFilament}) -> Real

Return the maximum distance $δ$ between neighbouring discretisation points of a filament.

If a vector of filaments is passed (variant 2), then the maximum among all filaments is returned.

Base.maximum(knot_increment, f::AbstractFilament) -> Real
Base.maximum(knot_increment, fs::AbstractVector{<:AbstractFilament}) -> Real

Return the maximum increment $Δt = t_{i + 1} - t_{i}$ between filament knots.

This is generally a good approximation for the maximum segment length, at least when the default ChordalParametrisation is used. If this is not the case, it is better to pass node_distance instead of knot_increment, which always returns the distance between two neighbouring nodes.

minimum_node_distance(f::AbstractFilament) -> Real
minimum_node_distance(fs::AbstractVector{<:AbstractFilament}) -> Real

Return the minimum distance $δ$ between neighbouring discretisation points of a filament.

This is equivalent to calling minimum(node_distance, f).

See also minimum and minimum_knot_increment.

minimum_knot_increment(f::AbstractFilament) -> Real
minimum_knot_increment(fs::AbstractVector{<:AbstractFilament}) -> Real

Return the minimum increment $Δt = t_{i + 1} - t_{i}$ between filament knots.

This is equivalent to calling minimum(knot_increment, f).

This is generally a good and fast approximation for the minimum segment length. However, this approximation is generally incorrect if one is using the FourierMethod discretisation method.

See also minimum, minimum_node_distance and maximum_knot_increment.

maximum_knot_increment(f::AbstractFilament) -> Real
maximum_knot_increment(fs::AbstractVector{<:AbstractFilament}) -> Real

Return the maximum increment $Δt = t_{i + 1} - t_{i}$ between filament knots.

This is equivalent to calling maximum(knot_increment, f).

This is generally a good approximation for the maximum segment length. See also maximum and minimum_knot_increment.

nodes(f::AbstractFilament{T}) -> AbstractVector{T}

Return the nodes (or discretisation points) $\bm{s}_i$ of the filament.

filament_length(f::AbstractFilament; quad = nothing) -> Real
filament_length(fs::AbstractVector{<:AbstractFilament}; quad = nothing) -> Real

Estimate the length of one or more filaments.

By default, the filament length is estimated using a straight segment approximation, which is fast but doesn't account for the actual curve geometry in-between discretisation points, underestimating the actual length. A quadrature rule may be optionally passed using quad (e.g. quad = GaussLegendre(4)) to obtain a more accurate result.

See also segment_length, which is used by this function.

Requirements

In the straight-segment implementation (quad = nothing), one needs the filaments to be closed, i.e. f[end + 1] should be set to the right value. This can be achieved either by update_coefficients! or close_filament!.

In the quadrature-based implementation (quad <: AbstractQuadrature), the interpolation coefficients must already have been computed via update_coefficients!.

Filaments.distance_to_field(ωf::Function, f::AbstractFilament) -> Real

Return an estimate of the normalised "distance" between a filament and a target vector field.

This function is meant to be used to verify the result of Filaments.from_vector_field, more specifically to verify that the filament is everywhere tangent to the objective vector field.

Returns 0 if the filament is perfectly tangent to the vector field at all discretisation points.

See Filaments.from_vector_field for more details.

update_coefficients!(f::AbstractFilament; knots = nothing)

Compute coefficients needed to perform inter-node interpolations and estimate derivatives.

Uses the current locations of the filament nodes. If nodes change, this function should be called to update the coefficients.

By default, this function also updates the parametrisation knots $t_i$ according to the current node positions. One can override this by passing a knots vector as a keyword argument. In particular, one can pass knots = knots(f) to keep the parametrisation knots unchanged.

This function will fail if the number of filament nodes is smaller than that required by the discretisation method. For instance, closed filaments discretised using cubic splines must have at least 3 nodes.

close_filament!(f::ClosedFilament) -> f

Set the filamend endpoint based on its starting point and its end-to-end offset.

This sets f[end + 1] = f[begin] + Δ⃗ where Δ⃗ is the end-to-end offset of the filament (see end_to_end_offset).

Calling this function is not needed if one has recently used update_coefficients!, which already closes the filaments. This function should be used as a cheaper alternative to update_coefficients!, when one only wants to close a filament (for instance to compute its length via filament_length). This function can also be used when the filament is no longer considered as "valid" by the chosen discretisation method, which happens when the number of nodes falls below some threshold (see check_nodes).

change_offset(f::ClosedFilament, offset::Vec3{<:Real}) -> f′

Change spatial offset between a filament start and endpoints.

See init for more details on optional spatial offsets.

This function is allocation-free. It returns a new filament which shares the same arrays as f, and only differs in the offset. Modifying nodes of the returned filament also modifies nodes of f.

fold_periodic!(Xs::AbstractVector{<:Vec3}, Ls::NTuple{3, Real}) -> Bool

Fold a set of coordinates onto the main unit cell.

The idea is that the filament coordinates are (mainly) in the main unit cell, given by $[0, L₁] × [0, L₂] × [0, L₃]$, where $Lᵢ$ is the period in each direction.

This is nice for visualisations, and might also improve performance of the folding required by long-range computations.

To avoid creating discontinuities in the filament, what this function actually does is to make sure that the average among all filament nodes is in the main unit cell.

Returns true if coordinates were modified, false otherwise. In the first case, one may want to call update_coefficients! to update the curve representation.

fold_periodic!(f::AbstractFilament, Ls) -> Bool

Fold filament nodes onto the main unit cell.

Returns true if coordinates were modified, false otherwise.

If true, coefficients may need to be updated afterwards using update_coefficients!.

redistribute_nodes!(f::AbstractFilament) -> f

Redistribute nodes of the filament so that they are (approximately) equally spaced.

More precisely, this function repositions the filament nodes such that the knot spacing $t_{i + 1} - t_{i}$ is constant. In other words, the new locations satisfy f[i] = f((i - 1) * Δt) where $Δt = t_{N + 1} / N$ is the knot spacing, $N$ is the number of nodes, and the index $N + 1$ refers to the filament endpoint (which is equal to the starting point for a closed filament).

One does not need to call update_coefficients! after calling this function.

split!(f::ClosedFilament, i::Int, j::Int; p⃗ = Vec3(0, 0, 0)) -> (f₁, f₂)

Split closed filament into two filaments.

Assuming j > i, the resulting filaments are respectively composed of nodes f[i + 1:j] and f[(j + 1:end) ∪ (begin:i)].

In practice, a split makes sense when the nodes f[i] and f[j] - p⃗ are spatially "close". Here p⃗ is an optional offset which usually takes into account domain periodicity (see also find_min_distance).

Note that this function modifies the filament f, which should then be discarded.

One should generally call update_coefficients! on both filaments after a split.

merge!(f::ClosedFilament, g::ClosedFilament, i::Int, j::Int; p⃗ = Vec3(0, 0, 0)) -> ClosedFilament

Merge two closed filaments into one.

The resulting filament is composed of nodes:

f[begin:i] ∪ {g[(j + 1:end) ∪ (begin:j)] - p⃗} ∪ f[i + 1:end]

Here p⃗ is an optional offset which usually takes into account domain periodicity (see also find_min_distance).

This function returns a new filament h which may share memory with f. The filament g is not modified.

One should generally call update_coefficients! on the returned filament after merging.

Merging a filament with its shifted image

This function supports merging a filament f with a shifted version of itself, g = f + p⃗. For this, one should simply pass g = f and a non-zero offset p⃗.

AbstractParametrisation

Abstract type representing a curve parametrisation.

For a parametric curve $\bm{s}(t)$, the choice of the discrete knots $t_i$ at given discrete nodes $\bm{s}_i$ is not unique. The parametrisation determines the way the knots $t_i$ are chosen.

Some predefined parametrisations include:

• ChordalParametrisation,
• CentripetalParametrisation,
• FourierParametrisation.

Moreover, CustomParametrisation allows to define custom parametrisation functions.

ChordalParametrisation <: AbstractParametrisation
ChordalParametrisation()

Represents the chordal parametrisation.

In this case, the increment between two knots is roughly equal to the length of the segment between the knots:

\[t_{i + 1} = t_{i} + |\bm{s}_{i + 1} - \bm{s}_i|\]

This is known as the chordal parametrisation in the context of Catmull–Rom splines.

This is the recommended (and default) parametrisation for general cases.

CentripetalParametrisation <: AbstractParametrisation
CentripetalParametrisation()

Represents the centripetal parametrisation.

In this case, the increment between two knots is given by:

\[t_{i + 1} = t_{i} + |\bm{s}_{i + 1} - \bm{s}_i|^{1/2}\]

This is known as the centripetal parametrisation in the context of Catmull–Rom splines.

FourierParametrisation <: AbstractParametrisation
FourierParametrisation()

Curve parametrisation adapted for Fourier series representation of curves.

In this case, the increment between two knots is constant and given by:

\[t_{i + 1} = t_{i} + \frac{2π}{N}\]

where $N$ is the number of nodes.

CustomParametrisation{F <: Function} <: AbstractParametrisation
CustomParametrisation(f::Function)

Allows to define a custom curve parametrisation.

Here the function f should have the signature f(Xs, i) → Δt where Xs is the vector of discretisation nodes. It must return the knot increment Δt = t_{i + 1} - t_{i}.

For instance, the ChordalParametrisation corresponds to defining

f(Xs, i) = norm(Xs[i + 1] - Xs[i])

Base.getindex(f::AbstractFilament{T}, i::Int) -> Vec3{T}
Base.getindex(f::AbstractFilament{T}, i::Int, ::Derivative{n}) -> Vec3{T}
Base.getindex(f::AbstractFilament{T}, i::Int, ::GeometricQuantity)

Return coordinates of discretisation point $\bm{s}_i$.

One may also obtain derivatives and other geometric quantities at point $\bm{s}_i$ by passing an optional Derivative or GeometricQuantity.

Base.setindex!(f::AbstractFilament{T}, v, i::Int) -> Vec3{T}

Set coordinates of discretisation point $\bm{s}_i$.

normalise_derivatives(ṡ::Vec3, s̈::Vec3) -> (s⃗′, s⃗″)
normalise_derivatives((ṡ, s̈)::NTuple)   -> (s⃗′, s⃗″)

Return derivatives with respect to the arc length $ξ$, from derivatives with respect to the parameter $t$.

The returned derivatives satisfy:

• $\bm{s}' ≡ t̂$ is the unit tangent vector;

• $\bm{s}'' ≡ ρ n̂$ is the curvature vector, where $n̂$ is the normal unit vector (with $t̂ ⋅ n̂ = 0$) and $ρ = R^{-1}$ is the curvature (and R the curvature radius).

normalise_derivatives!(Ẋ::AbstractVector, Ẍ::AbstractVector)

Normalise vectors containing derivatives at filament locations.

If possible, prefer using normalise_derivatives, which works on a single filament location at a time.

See normalise_derivatives for more details.

integrate(integrand::Function, f::AbstractFilament, i::Int, quad::AbstractQuadrature; limits = nothing)
integrate(integrand::Function, s::Segment, quad::AbstractQuadrature; limits = nothing)

Estimate integral along a filament segment using the chosen quadrature.

Integration is performed in the segment (f[i], f[i + 1]).

The function integrand(ζ) takes the arguments f, i and $ζ ∈ [0, 1]$. The first two are a bit redundant since they're the same filament and node index passed to this function, while ζ corresponds to the position of a point within the segment. See further below for some examples.

By default, the integration is performed along the whole segment, that is, for the range $ζ ∈ [0, 1]$. It is possible to integrate over a subset of the segment, i.e. for $ζ ∈ [a, b]$ with $0 ≤ a ≤ b ≤ 1$. For this, one should pass the keyword argument limits = (a, b).

Examples

Estimate arc length of segment $[i, i + 1]$, given by $ℓ = ∫_{t_{i}}^{t_{i + 1}} |∂_t \bm{X}(t)| \, \mathrm{d}t$:

quad = GaussLegendre(4)  # quadrature rule
ℓ = integrate(f, i, quad) do f, i, ζ
    norm(f(i, ζ, Derivative(1)))  # = |∂ₜ𝐗|
end

Alternatively:

quad = GaussLegendre(4)  # quadrature rule
s = Segment(f, i)
ℓ = integrate(s, quad) do s, ζ
    norm(s(ζ, Derivative(1)))  # = |∂ₜ𝐗|
end

integrate(integrand::Function, f::AbstractFilament, quad::AbstractQuadrature; limits = nothing)

Estimate integral over a whole filament.

The integral is computed as the sum of the integrals over each filament segment. The given quadrature rule is applied over each segment, meaning that the total number of evaluations is given by the number of segments multiplied by the length of the quadrature rule.

The signature of the integrated function must be integrand(f::AbstractFilament, i::Int, ζ::Real), where i is the index of the segment of interest. See below for some examples.

Examples

Estimate the total length of a closed filament, $L = ∮ |∂_t \bm{X}(t)| \, \mathrm{d}t$:

quad = GaussLegendre(4)  # quadrature rule

integrand(f, i, ζ) = norm(f(i, ζ, Derivative(1)))  # = |∂ₜ𝐗|

# Here `f` is an existent filament.
L = integrate(integrand, f, quad)

# Or, more conveniently, using a `do` block to define an anonymous function.
L = integrate(f, quad) do ff, i, ζ
    norm(ff(i, ζ, Derivative(1)))  # = |∂ₜ𝐗|
end

find_min_distance(
    fx::AbstractFilament, fy::AbstractFilament, i::Int, j::Int;
    periods::NTuple{3, Real} = (Infinity(), Infinity(), Infinity()),
    maxiter = 4, rtol = 1e-2,
)

Determine the minimum distance between filament segments.

This function estimates the minimum distance between filament segments fx[i:i+1] and fy[j:j+1] via an iterative (Newton–Raphson) method.

The filaments fx and fy may be the same filament.

Returns

Returns a NamedTuple with the following fields:

• ζx, ζy: optimal locations in $[0, 1]$ within each segment;

• x⃗, y⃗: optimal locations within each segment;

• p⃗: periodic offset (each component is a multiple of the domain period along that direction);

• d⃗: minimum distance vector, d⃗ = x⃗ - y⃗ + p⃗.

Optional keyword arguments

• periods: the period of the spatial domain. This should be given if one wants to take into account periodic images of the filaments;

• maxiter = 4: maximum number of iterations of the Newton method;

• rtol = 1e-2: relative tolerance for stopping the iterations. By default it's not very small since in general we don't need to be very precise.

check_nodes(Bool, f::AbstractFilament) -> Bool
check_nodes(f::AbstractFilament)

Check whether current filament nodes are compatible with the filament discretisation method.

In its first form, this function returns false in case of incompatibility, while it throws an error in its second form.

For now, the only requirement is that the number of nodes must be larger than some small value. In particular, one can't have a closed filament with less than 3 nodes (but the specific discretisation method might impose some other small value).

number_type(x) -> Type{<:Number}

Obtain the number type associated to a container x.

This is expected to return a concrete type T <: Number.

This function can be useful when the actual number type is hidden behind many nested array types.

Some examples:

• if x is a single AbstractFilament{T}, this returns T;
• if x is a vector or tuple of AbstractFilament{T}, this also returns T;
• if x is a vector of vectors of SVector{3, T}, this returns T.

Note that, in the last two cases, this corresponds to eltype(eltype(eltype(x))), which is less readable and prone to errors.

Filaments.curl(vf::Function) -> Function
Filaments.curl(vf::Function, x⃗::Vec3) -> Vec3

Return curl of a vector-valued field, $\bm{ω}(\bm{x}) = \bm{∇} × \bm{v}(\bm{x})$.

Here vf is a function $\bm{v}(\bm{x})$ which takes a 3D vector x⃗ and returns a vector value v⃗.

The curl of vf is obtained via automatic differentiation using ForwardDiff.jl.

The first variant returns a function which can be then evaluated at any x⃗. The second variant directly evaluates the curl at some x⃗.

RefinementCriterion

Abstract type describing a curve refinement criterion.

Implemented refinement criteria are:

• NoRefinement: disables refinement;

• RefineBasedOnSegmentLength: enforces a minimum and maximum distance between neighbouring filament nodes;

• RefineBasedOnCurvature: inserts more nodes on highly-curved filament segments.

NoRefinement <: RefinementCriterion
NoRefinement()

Used to disable filament refinement.

RefineBasedOnSegmentLength <: RefinementCriterion
RefineBasedOnSegmentLength(ℓ_min, ℓ_max = 2 * ℓ_min)

Refinement criterion imposing a minimum segment length.

This refinement criterion imposes neighbouring filament nodes to be at a distance $ℓ ∈ [ℓ_{\min}, ℓ_{\max}]$. This means that:

• nodes are inserted if the distance between two nodes is $ℓ > ℓ_{\max}$. The insertion is done at an intermediate position using the functional representation of the filament (e.g. splines or Hermite interpolation);

• nodes are removed if the distance between two nodes is $ℓ < ℓ_{\min}$. For a filament which is strongly curved at that point, this means that local information is lost and that the filament is smoothed.

RefineBasedOnCurvature <: RefinementCriterion
RefineBasedOnCurvature(ρℓ_max::Real, ρℓ_min = ρℓ_max / 2.5; ℓ_min = 0.0, ℓ_max = Inf)

Curvature-based refinement criterion.

Node insertion

According to this criterion, a filament is locally refined if:

\[ρ \left|\bm{X}_{i + 1} - \bm{X}_{i}\right| > (ρℓ)_{\max}\]

where $ρ$ is some estimation of the curvature of segment $[i, i + 1]$.

For splines, this uses a classical knot insertion algorithm which preserves the shape of the curve.

Node removal

Similarly, filaments nodes are removed based on the value of ρℓ_min. This value should be less than ρℓ_max / 2 to avoid alternatively adding and removing nodes when repeatedly calling refine!.

For safety, two adjacent nodes will never be removed in a single call to refine!.

Note that, when filaments are nearly straight, this may lead to the removal of most nodes. To limit this, set the keyword argument ℓ_max to some finite value determining the maximum length of a segment. Similarly, the keyword argument ℓ_min sets a lower limit for the distance between neighbouring nodes.

refine!(f::AbstractFilament, crit::RefinementCriterion) -> (Int, Int)

Refine the filament according to a given criterion.

More precisely, this function can add and remove discretisation points according to the chosen refinement criterion.

Returns the number of added and removed nodes.

Example usage:

crit = RefineBasedOnCurvature(0.5)
refine!(f, crit)

insert_node!(f::AbstractFilament, i::Integer, [ζ = 0.5]) -> Vec3

Insert node in-between locations f[i] and f[i + 1].

The optional argument ζ ∈ [0, 1] corresponds to the relative location of the new node within the segment. By default it is set to ζ = 0.5, which corresponds to an estimation of the middle of the segment.

Note that update_after_changing_nodes! must be called after inserting one or more nodes.

See also remove_node!.

Returns the inserted node.

remove_node!(f::AbstractFilament, i::Integer) -> Vec3

Remove node at location f[i].

Note that update_after_changing_nodes! must be called after removing one or more nodes.

See also insert_node!.

Returns the removed node.

update_after_changing_nodes!(f::AbstractFilament)

Update filament fields after changing nodes.

Depending on the filament discretisation method, this can recompute derivatives, knots or discretisation coefficients.

Should be called after inserting or removing filament nodes. See insert_node! and remove_node!.

GeometricQuantity

Abstract type defining a geometric quantity defined for a filament.

Some available geometric quantities include:

• UnitTangent,

• CurvatureVector,

• CurvatureScalar,

• CurvatureBinormal,

• TorsionScalar.

Evaluating geometric quantities works in the same way as evaluating derivatives.

UnitTangent <: GeometricQuantity

Represents the unit tangent vector $\bm{T}$ at a filament location.

In terms of an arbitrary parametrisation $\bm{s}(t)$ (where $t$ is in general different from the arc length $ξ$), the unit tangent vector is

\[\bm{T} = \frac{\bm{s}'}{|\bm{s}'|}\]

where derivatives are with respect to the parameter $t$.

CurvatureVector <: GeometricQuantity

Represents the curvature vector associated to a filament.

In terms of an arbitrary parametrisation $\bm{s}(t)$ (where $t$ is in general different from the arc length $ξ$), the curvature vector is

\[\bm{ρ} = \frac{{|\bm{s}'|}^2 \bm{s}'' - (\bm{s}' ⋅ \bm{s}'') \, \bm{s}'}{|\bm{s}'|^4} = \frac{\bm{s}' × (\bm{s}'' × \bm{s}')}{|\bm{s}'|^4}\]

where derivatives are with respect to the parameter $t$.

The curvature vector can be written as $\bm{ρ} = ρ \bm{N}$ where $\bm{N}$ is the unit normal vector and $ρ$ the scalar curvature (the inverse of the curvature radius).

CurvatureScalar <: GeometricQuantity

Represents the scalar curvature associated to a filament.

This is simply the norm of CurvatureVector. It is explicitly given by

\[ρ = \frac{\bm{s}' × \bm{s}''}{|\bm{s}'|^3}\]

where derivatives are with respect to the arbitrary parametrisation $\bm{s}(t)$.

CurvatureBinormal <: GeometricQuantity

Represents the scaled binormal vector associated to a filament.

The scaled binormal vector is defined as $\bm{b} = \bm{T} × ρ⃗ = ρ \, (\bm{T} × \bm{N}) = ρ \bm{B}$, where $\bm{B}$ is the (unit) binormal vector and $ρ$ is the scalar curvature.

TorsionScalar <: GeometricQuantity

Torsion of a filament.

The torsion $τ$ describes the variation of the binormal vector along the curve.

It can be obtained as

\[τ = \frac{(\bm{s}' × \bm{s}'') ⋅ \bm{s}'''}{|\bm{s}' × \bm{s}''|^2}\]

where derivatives are with respect to an arbitrary curve parametrisation.

SegmentIterator{Filament <: AbstractFilament}

Convenience type allowing to iterate over the segments of a filament.

segments(f::AbstractFilament) -> SegmentIterator(f)

Create a SegmentIterator object for iterating over the segments of a filament.

Base.length(s::SegmentIterator{<:AbstractFilament}) -> Int

Return the number of segments in a filament.

Base.eachindex(s::SegmentIterator{<:AbstractFilament}) -> AbstractRange

Return the indices associated to the segments of a filament.

Segment{<:Filament}
Segment(f::AbstractFilament, i::Integer)

Represents a single filament segment.

The segment goes from nodes f[i] to f[i + 1].

midpoint(s::Segment) -> Vec3

Return an estimation of the segment midpoint (prioritising performance over accuracy).

Filaments.segment_length(s::Segment; quad = nothing)

Estimate length of a filament segment.

One may pass a quadrature rule as quad for better accuracy. Otherwise, if quad = nothing, this simply returns the straight distance between the two segment extremities.

filamentplot!([ax,] f::AbstractFilament, [velocities]; kws...)
MakieCore.plot!(ax, f::AbstractFilament, [velocities]; kws...)

Plot filament onto existent 3D axis.

The first argument should typically be an Axis3 or an LScene.

Example usage:

using GLMakie
fig = Figure()
ax = Axis3(fig[1, 1])
plot!(ax, f)  # f is a filament

Optional arguments and their defaults

• refinement::Int = 1: level of refinement of the curves (must be ≥ 1)

• periods = (nothing, nothing, nothing). This can be a tuple of values representing the domain period, for instance (2π, 2π, 2π). In that case, filaments will be "broken" when they exit the main unit cell, so that all vortex elements are within the cell.

• color = :black

• linewidth = 1.5f0

• linestyle = :solid

• markercolor = nothing (nothing → same as color)

• marker = :circle

• markersize = 10.0f0

• colormap = :viridis

• colorrange = MakieCore.Automatic()

Arrow arguments

The following are used when plotting arrows (tangents, curvatures, velocities, …):

• arrowscale = 1.0f0 allows to scale vectors (controls their length). Corresponds to lengthscale in Makie.arrows;

• arrowsize = MakieCore.Automatic() controls the head size. It has the same name in Makie.arrows. It should be a tuple (sx, sy, sz), where the first 2 set the width of the cone, and sz its height;

• arrowwidth = MakieCore.Automatic() sets the linewidth of the arrow. Corresponds to linewidth in Makie.arrows.

See also the Makie docs on arrows.

Plotting tangent and curvature vectors

Tangent and curvature vectors can be optionally plotted via the tangents and curvatures arguments. A single vector will be plotted for each filament segments. By default, vectors are evaluated at filament nodes, but one can also evaluate them in-between nodes using the vectorpos argument.

• tangents::Bool = false: plot unit tangent vectors.

• curvatures::Bool = false: plot curvature vectors. Note that the magnitude is the local curvature $ρ = 1 / R$, where $R$ is the curvature radius.

• tangentcolor = nothing

• curvaturecolor = nothing

• vectorpos = 0.0: relative vector positions within each segment. Must be in $[0, 1]$.

Plotting velocities of filament nodes

Similarly, it is possible to plot vector quantities attached to filament nodes, such as filament velocities. For this pass a vector of velocities as a positional argument after the filament f.

Associated keyword arguments:

• velocitycolor = nothing colour of velocity vectors

filamentplot(f::AbstractFilament, [velocities]; kws...)
MakieCore.plot(f::AbstractFilament, [velocities]; kws...)

Plot a filament using Makie.jl.

Example usage:

using GLMakie
plot(f; refinement = 4)  # f is a filament

See filamentplot! for details and for optional keyword arguments.