% \chapter{Numerical simulations of implicit curve dynamics}
\chapter{Hidden degrees of freedom in implicit vortex filaments}\label{ch:hidden_degrees}

\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{images_hidden_degrees/torus2/TorusKnotComp.png}
    \caption{The time evolution of  highly knotted vortex filaments. Our implicit description represents the dynamics of filaments as the zero levelsets of space-time complex-valued functions, which automatically handles topological changes of curves without yielding singularities.} 
    \label{fig:torus}
\end{figure}

This chapter is largely a reprint of the article:

\bibentry{iwc2022implicit_filaments}.

This chapter develops a computational method for the dynamics of space curves, with emphasis on its applications to vortex filaments in fluid dynamics. Instead of representing these filaments with explicit curve geometry and Lagrangian equations of motion, we describe them using implicit representations via a class of complex functions. We leverage their redundant degrees of freedom in both the configuration and the dynamics, which can be tailored specifically to mitigate numerical instability that native approaches suffer from. 

% Furthermore, we note how these hidden degrees of freedom perfectly map to a {\em Clebsch variables} from symplectic geometry.
% Motivated by these observations, 
For handling chaotic systems specifically,
we introduce {\em untwisted} Clebsch variables and {\em non-swirling} dynamics, which regularize sources of numerical instability, particularly in the twisting modes around curve filaments.
Built on these, our simulation method robustly produces the dynamics of large numbers of interacting vortex filaments and effortlessly handles topological changes and re-connection events.


\paragraph{Changes from the original publication}
This chapter contains a few updates from the article \cite{iwc2022implicit_filaments}. We provide a more precise definition of the equivalence class of complex functions that represent the same space curve. The discussion and results concerning the twist of complex functions have also been revised to reflect this refined definition.

\paragraph{Dependency of the chapter}
This chapter is self-contained, as it reintroduces the necessary concepts from previous chapters. Readers interested in the theoretical foundations of implicit representations may find \autoref{ch:implicit_representations} useful.

\section{Introduction}

The deformation of space curves is an interesting topic in many subjects such as differential geometry, low-dimensional topology, classical and quantum fluid mechanics, and electromagnetism. One example from fluid mechanics is the dynamics of vortex filaments.
In a nearly inviscid fluid, vorticity originates from {\em codimension-1} interfaces or obstacle surfaces. The vortex sheets subsequently roll up into {\em codimension-2} vortex filaments, due to the Kelvin--Helmholtz instability.  Hence, most physical inviscid fluids have their vorticity concentrated into a sparse set of space curves, rather than distributed evenly throughout space.
Based on this observation, certain physical equations model fluids only with dynamically deforming space curves. Many fluid simulation methods take advantage of this sparsity structure.

One major challenge for an explicit (Lagrangian) filament-based fluid solver is to handle reconnection events when filaments collide.
Without any reconnection, the total length of filaments can grow exponentially, exploding the computational cost and halting the solver. Hence, existing explicit filament simulators include a tedious process of collision detection followed by non-differentiable heuristic geometry surgeries.

To that end, implicit (Eulerian) curve representations are more appealing.
The recently emerging \emph{Clebsch representation} expresses vortex lines as level sets of a 2-dimensional-valued function called \emph{Clebsch variables} 
\cite{Clebsch:1859:IHG}.
Like any level set method, topological changes of level set geometries occur gracefully.  
The difficulty, however, in a Clebsch-based fluid solver is in the dynamics of the Clebsch variables.  The Clebsch variables satisfy the transport equation advected with the fluid velocity, which unfortunately behaves in a swirling motion with a high-spatial frequency and singularities near the vortex filaments.
Such a rough transporting vector field is hard to resolve accurately in a computational grid.  Even if the transport equation is computed accurately, the level set function will quickly evolve into a twisted and distorted function that is difficult to deal with.

This chapter develops a new approach for describing the geometry and dynamics of filaments with implicit curve functions. Our main insight is that the problem has a huge number of redundant degrees of freedom:  both the velocity field and the level set function (i.e., the Clebsch representation) can be varied in ways that do not change the solution. We exploit these additional degrees of freedom to ensure stable numerical simulation and automatic handling of topological changes, without sacrificing  accuracy. 
In particular, we choose an {\em untwisted} Clebsch representation for the level set geometry, and {\em non-swirling} dynamics for advecting  vortex filaments. We regularize these functions by identifying and constraining hidden degrees of freedom in their representations, allowing us to greatly improve numerical robustness compared to naive implementations. 

Our algorithm is the first method for animating implicit vortex filament geometry with automatic topological changes. Our mathematical formulation offers new tools for future research on fluid simulation and curve geometry processing, and our results show greatly improved stability compared to a standard level set implementation, with fewer user parameters than for explicit Lagrangian filament techniques.

%!TEX root = VortexUntwisted_paper.tex
\section{Related work}

We review previous work, highlighting numerical methods for simulating the time evolution of vorticity in incompressible fluids.
These vortex-capturing schemes seek a useful representation of the vorticity and solve its governing equations of motion.

\paragraph{Explicit vortex methods}
In previous vortex methods, vortices are represented either as
particles \cite{gamito1995two,park2005vortex,Selle:2005:VPM,zhang2014pppm,Angelidis:2017:MVF},
filaments \cite{Cottet:2000:VM,Angelidis:2005:SSB,Weissmann:2009:FBS,Weissmann:2010:filaments,Padilla2019bubbleRings},
segments \cite{Chorin:1990:HRV,Xiong:2021:VSC}, 
sheets \cite{brochu2012linear,Pfaff:2012:LVS,da2015double} or 
volumes \cite{Elcott:2007:SCP,Zhang:2015:RMV}.
Vortex particle methods
represent vortices as a disconnected point cloud. However, the strength of vortex per particle or per unit volume \cite{Zhang:2015:RMV} undergoes a numerically unstable \emph{vortex stretching}, requiring an artificial clamping or diffusion that sacrifices accuracy.
The stretching problem is avoided by representing vorticity per filament, segment, or sheet, or per unit area using differential 2-forms \cite{Elcott:2007:SCP}.  However, describing vortex explicitly (Lagrangian method) with filaments, segments and sheets comes with a cost of sophisticated and heuristic treatment for changes of vortex topology.  Volumetric (Eulerian) methods \cite{Elcott:2007:SCP} do not require managing topological changes, but they do not have a handle on codimensional structures. 
Our codimension-2 level set method is an Eulerian method that can represent sharp filament structures without any additional difficulty from vortex reconnection.

\paragraph{Clebsch representations}
Another Eulerian representation of vorticity is to describe vortex lines as the level sets of a 2D-valued function known as the Clebsch variable \cite{Clebsch:1859:IHG,Lamb:1895:HD}.  The representation was not widely adopted since all \(\RR^2\)-valued Clebsch variables can only describe fluids with zero \emph{helicity} \cite{Chern:2017:IF}.  The helicity problem is solved by using a sphere-valued Clebsch variable \cite{Kuznetsov:1980:TMC} which can represent nonzero (though quantized) values of helicity. Since the recognition of these ``spherical'' Clebsch maps, they have become an established method for vortex representation in computer graphics \cite{chernSchroedinger,Chern:2017:IF,Yang:2021:CGF}.  However, Clebsch variables represent a {\em smooth vorticity field} by the continuum of level sets of all values. It remains a challenge to represent a {\em sharp codimension-2 filament}, especially with a limited grid resolution. By contrast, our method represents a sharp filament just as the \emph{zero set} of a %\(\CC\)-valued 
complex-valued
function.   Moreover, we show in \autoref{sec:UntwistedClebschVariables} 
that our representation can actually resolve continuous values of helicity (as opposed to the quantized values mentioned earlier).

\paragraph{Dynamics of Clebsch variables}
In addition to the implicit representation of vorticity, Clebsch variables also play significant roles in a variational and Hamiltonian formulation for the incompressible Euler equation \cite{Clebsch:1859:IHG,Lamb:1895:HD,Morrison:1998:HDI,Chern:2017:FD}.
In short, one (of many possible) governing equation(s) for the Clebsch variable is exceedingly simple: the Clebsch variable is advected by the fluid velocity.  This equation of motion is recently adopted by ~\cite{Yang:2021:CGF}.  However, referred to as the \emph{Lagrangian chaos}, a direct transportation by the fluid velocity quickly stirs and twists any variable to a distorted one unresolved by the finite computational grid \cite{Qu:2019:ECF}. The Clebsch variable is no exception under such dynamics.
The method of \emph{Schr\"odinger's Smoke} \cite{chernSchroedinger} bypassed the Lagrangian chaos: its total energy (Hamiltonian) includes the Dirichlet energy of the Clebsch variable, which is therefore bounded for all time.  
However, while the dynamics of \emph{Schr\"odinger's Smoke} appear to be similar to that of Euler's equation, it is only an approximation to the Euler fluid. 
There is still a large degree of freedom in the Clebsch representation and its dynamical system.  Finding an equation of motion for the Clebsch variable that both describes the correct Euler fluid without Lagrangian chaos is an unexplored research topic.  This chapter describes an instance of a Lagrangian-chaos-free dynamical system for an implicit representation of vortex filaments. 


\paragraph{Implicit filament representations}
We represent filaments as the zero set of a complex-valued function. These zero sets of a complex phase field are widely studied in condensed matter physics as \emph{topological defects} appearing in superfluids and superconductors \cite{Bethuel:1994:GLV,Pismen:1999:VNF}. These topological defect models also facilitate singularity placements in flow analysis and geometry processing \cite{weissmann2014smokerings,Solomon:2017:BEO,Palmer:2020:ARV}. 
Complex phase field models are taken more generally as high-codimensional level set representations by \cite{Ambrosio:1996:LSA,MBO2001diffusion,Burchard:2001:MOC,Min:2004:LLS}.
%However, in the physics or geometry processing literature, the level set functions have other physical meaning, leaving little rooms for a smoother representative.
However, in the physics and geometry processing literature, the level set functions have specific physical meanings such as the phase of a wave, leaving little room for a smoother representative. %\cCW{$\leftarrow$ I didn't understand this sentence, so I tried to make it clearer. Please tell me (or fix it) if I got the meaning wrong.} 
In most cases, the phase field has norm 1 except for a sudden dip to 0 near the filaments, creating a configuration that is difficult to resolve efficiently on a computational grid. 
In the level set method literature,
the norm-1 condition is often adopted for (re)initializing the level set functions despite the discontinuity \cite{MBO2001diffusion}. The complex phase is constructed locally with little discussion about global topological obstruction. \cite{Burchard:2001:MOC} addressed the challenges in reinitializing the multi-component level set functions; by mimicking the codimension-1 signed distance functions, they propose a sophisticated reinitialization by solving a ``manifold eikonal equation'' along the isosurface of each function component. Unfortunately the process will not resolve the twists of the framed curve.
To our knowledge, there has not been a thorough discussion about the degrees of freedom in the implicit filament representations or in their dynamics until now.

In this chapter, we explore the degrees of freedom of both the codimension-2 level set functions and their equations of motion.  We further provide a simple reinitialization method comparable to the codimension-1 signed distance function.  Comparisons show that exploiting these degrees of freedom are essential to a robust simulation.


\section{Representations for evolving curves}
\label{sec:representation}
We now begin our description of implicit filament dynamics.
The main mathematical object is a union of closed space curves.  This section describes an implicit representation for these curves and their dynamics, as well as the degrees of freedom in the representation. 

\subsection{Representations for curve configuration}
Let the physical domain be an open region \(M\subset\RR^3\).  
We use $\gamma$ to represent a collection of $m$ closed curves
\begin{align}
    \label{eq:gamma}
     \gamma\colon\textstyle\bigsqcup_{i=1}^m\SS^1\rightarrow M,
\end{align}
where $\bigsqcup$ denotes a disjoint union, \(\SS^1\) is the topological circle, $m$ is the number of filaments, and $\gamma$ is the mapping from the 1D curves into 3D.
The configuration space \(\cX\) of these filaments is the space of all possible placements of these curves:\footnote{Precisely, reparametrizations mean the quotient by the orientation preserving diffeomorphisms on $\bigsqcup_{i=1}^m\SS^1$ for each fixed $m$ (See \autoref{sec:explicit_shapes} for detail).}
\begin{align}
    \label{eq:StateSpace}
    \cX = 
    \left.
    \bigsqcup_{m=1}^\infty\left\{
    \gamma\colon\textstyle\bigsqcup_{i=1}^m\SS^1\rightarrow M
    \right\}
    \right/\text{reparametrizations}.
\end{align}
Although it is straightforward to represent curves {\em explicitly} via objects \(\gamma\) in \eqref{eq:StateSpace} as parameterized curves 
or their discrete counterparts, topological changes of curves such as splitting or merging are difficult to describe mathematically and algorithmically.  For instance, the number \(m\) of components can change when curves reconnect or split apart.

Instead of relying on an explicit curve representation, our work adopts an {\em  implicit} representation for the elements in \(\cX\).
We model every collection of closed curves in \(M\) as the zero set of a complex-valued level set function \(\psi\colon M\rightarrow\CC\):
\begin{align}
    \gamma = \{p\in M\,\vert\,\psi(p) = 0\} = \{\Re\psi = 0\}\cap\{\Im \psi=0\},
\end{align}
In other words, an alternative definition for $\gamma$ is the set of all points $p$ where {\em both} the real and imaginary components of a level set function $\psi$ evaluate to zero as in \autoref{fig:level_surfaces}. We find it useful to draw an analogy to the scalar-valued level sets commonly used in computer graphics applications: the zero level-set of a scalar-valued function represents shapes of codimension 1 ({\em a.k.a.} surfaces), while our complex-valued level set has twice as many variables and thus represents shapes of codimension 2 (curves).

We next note that different functions $\psi$ can represent the same collection of curves if they share the same zeros and  orientation. 
To formulate this redundancy precisely, we define the following equivalence relation \(\sim\) on the function space \((M\rightarrow\CC)\).
We say \(\psi_1\sim\psi_2\) if their zeros agree \ie, $\psi_1^{-1}(0)=\psi_2^{-1}(0)$, and share the same orientation, that is, 
\begin{align}
    \int_{\partial \Sigma} \Im \frac{d\psi_1}{\psi_1} = \int_{\partial \Sigma} \Im \frac{d\psi_2}{\psi_2}
\end{align}
for any oriented topological disc $\Sigma$ intersecting transversely with $\gamma$. For details, see \autoref{sec:implicit_representations}.

% and only if \(\psi_1 = \varphi\psi_2\) for some nowhere-vanishing function \(\varphi\colon M\to\CC\setminus\{0\}\).

% This means \(\psi_1\sim\psi_2\) if and only if they share the same zero set.

\begin{figure}[tbp]
         \centering
\includegraphics[width=0.85\linewidth, 
,trim={0cm 0cm 0cm 0cm}
, clip]
{{images_hidden_degrees/illustrations/level_surfaces_and_ribbons/level_surface_re_im_3}.png}
    \caption{An example of level surfaces of two linked curves. The blue and red surfaces are cross sections of \(\{\textup{Im} \psi =0 \}\) and  \(\{\textup{Re} \psi =0 \}\) respectively. 
    % \cSI{If you think we don't need all of them. Feel free to use only the one with both real and imaginary level surfaces.}
    }\label{fig:level_surfaces}
\end{figure}


% \cCW{I removed the mathematical definition of congruence between configuration spaces since it is fairly obvious from the writing and only adds potential confusion if you don't take the effor to define and explain your symbols. I also removed the sentence following it about the free degree of freedom, since it is redundant with the text above.}
 
% \cCW{Update after meeting Sadashige: The equation said more than I realized, so need to find a way to restore the lost information while giving the reader complete confidence that he/she knows what is the important information they should take away from it.}
 

 
% \cAC{I moved the regularity assumption up here}
 
%\subsubsection{Regularity}
%\label{sec:Regularity}
 For the implicit representation to work properly, we make a regularity assumption that the level set functions \(\psi\) are sufficiently ``nice'' around the zero set.  Precisely, we assume \(\psi\) is smooth, and
 the differential \(d\psi\vert_p\colon T_pM\rightarrow\CC\cong\RR^2\) is surjective, i.e. its matrix form has rank 2 at any point \(p\)  on the zero-set.
%  it's first spatial derivative has rank 2 at all points on the zero-set.
%\cSI{The differential being a linear map between vector spaces and its surjectivity are important, so I retrieved this.}
% 
% \footnote{The full-rank requirement does not modify the degrees of freedom in the rescaling \(\varphi\) since \(d(\varphi\psi)\vert_p = d\varphi\vert_p\underbrace{\psi_p}_{=0} + {\varphi_p} d\psi\vert_p = \varphi_p d\psi\vert_p\) and a non-zero complex multiplication by \(\varphi_p\) always preserves the rank of \(d\psi\vert_p\).}
%\cCW{simplified the explanation in the text and demoted the subsubsection to a paragraph. The footnote and remaining text still preserve the mathematical info needed for later proofs.}

%  \cAC{How about the following:}
To summarize, our configuration space of filaments \(\cX\) from \eqref{eq:StateSpace} is replaced by the space of complex-valued functions \(\psi\) modulo the above equivalence class:
 \begin{align}
     \label{eq:GDefinition}
      \cF \coloneqq \{\psi\colon M\rightarrow\CC\}/\sim.
  \end{align}
While \(\cX\) and \(\cF\) describe the same configuration space of filaments, the objects in \(\cF\) are much more continuous compared to the disjoint spaces of \(\cX\).  

\paragraph{Relation to Clebsch representations}
The representation of codimension-2 curve geometries in 3D is known in fluid dynamics as \emph{Clebsch representations} \cite{Clebsch:1859:IHG,Lamb:1895:HD,Chern:2017:IF,Yang:2021:CGF}. 
For a fluid flow with a smooth vorticity field, the vortices are geometrically depicted as fibrous vortex lines diffusely distributed over the fluid domain.  A Clebsch representation aims at an implicit representation for such fibrous structure.  The representation uses a map \(s\colon M\rightarrow\Sigma\) from the 3D fluid domain \(M\) to a 2-dimensional manifold \(\Sigma\) with a measure \(\sigma\) 
%\cAC{just \(\sigma\) and remove \(\Omega^2(\Sigma)\)}
 to describe the \emph{vortex lines} as preimages \(s^{-1}\{p\}\) of points \(p\in\Sigma\), and the density of the vortex lines 
%  \sout{(\ie\@ the vorticity 2-form)} \cAC{\(\leftarrow\) remove} 
 as the pullback \(s^*(\sigma)\) of the measure \(\sigma\).
A smooth Clebsch map \(s\) and a smooth measure \(\sigma\) yields a smooth distribution of vortex lines.  To achieve a more singular and concentrated vorticity field such as vortex filaments, one would consider \(s\) with larger derivatives \cite{MarsdenWeinstein1983coadjoint} (\(s\) sweeps out more measure \(\sigma\) over a small area in \(M\)). 

In our setup, we want to represent singular curves with a Dirac-\(\delta\) density, instead of diffused distribution of vortex lines.  Previous considerations in Clebsch representations would set the Clebsch map \(s\) with enormous derivative.  By contrast, we obtain such concentrated filaments by setting \(\sigma\) singular while keeping the Clebsch map smooth.   Our complex level function \(\psi\colon M\rightarrow\CC\) is a Clebsch map with the target space \(\Sigma = \CC\) equipped with a \(\delta\)-measure \(\sigma = \delta_0\) at the origin.

An important discussion about Clebsch representations \cite{Chern:2017:IF} is whether or not a fluid configuration can be represented with the choice of \(\Sigma\) and \(\sigma\).
Previous Clebsch representations \cite{Chern:2017:IF,Yang:2021:CGF} adopt \(\Sigma=\SS^2\), since the more straightforward choice of \(\Sigma = \RR^2\cong\CC\) with the standard area measure \(\sigma\) can only represent fluid flows with zero \emph{helicity}.  The helicity obstruction is reduced for \(\Sigma=\SS^2\) as it can admit a discrete set of nonzero helicity \cite{Chern:2017:IF}.  Our Clebsch representation can represent any space curve without any obstruction.  In particular, the helicity of a vortex filament is proportional to its \emph{writhe} \cite{ArnoldKhesin2021} which can take any real value.

\subsection{Representations for curve dynamics}
\label{sec:RepresentationsForFilamentDynamics}

In the explicit representation, a first-order time evolution of curves  \(\gamma\colon(\bigsqcup^{m}\SS^1)\times\RR\rightarrow M\) can be described by an  equation of the form:
\begin{align}
    \label{eq:ExplicitDynamics}
    {\partial\gamma\over\partial t}(s,t) =    V_{\gamma_t}(s),\quad s\in\textstyle\bigsqcup\SS^1, t\in\RR.
\end{align}
Here, $s$ is the parameterization of the curve, $t$ is time, and the \emph{velocity} \(   V_{(\cdot)}\colon\cX\rightarrow (\bigsqcup\SS^1\rightarrow\RR^3)\) is a dynamical model that tells the filament how to move based on the current filament shape and position. 

\begin{example}
    \label{ex:BiotSavart}
In the context of fluid dynamics, important examples for the velocity model \(   V\) are the ones that govern the motion of vortex filaments.  When \(M=\RR^3\), \ie\@ there are no obstacles or boundaries, the velocity models are the Biot--Savart model
\begin{align}\label{eq:biot_savart}
        V_{\gamma}^{\on{BS}}(s) \coloneqq \frac{\Gamma}{4 \pi}  \oint \gamma^{\prime}(\tilde s) \times \frac{\gamma(s)-\gamma(\tilde s)}{\left|\gamma(s)-\gamma(\tilde s)\right|^{3}} d\tilde s
\end{align}
and the more regular Rosenhead--Moore model \cite[pp.~213]{Saffman1993vortex}
\begin{align}\label{eq:Rosenhead-Moor}
        V_\gamma^{\on{RM}}(s)\coloneqq \frac{\Gamma}{4 \pi}  \oint \gamma^{\prime}(\tilde s) \times \frac{\gamma(s)-\gamma(\tilde s)}{\sqrt{e^{-\nicefrac{3}{2}}a^{2}+\left| \gamma(s)-\gamma(\tilde s)\right|^{2}}^{3}} d\tilde s
 \end{align}
where the constants $\Gamma$ and $a$ are the vortex strength and vortex thickness respectively, and the integrating measure \(d\tilde s\) is the arclength element (set \(a=0\) for the Biot--Savart model). 
Note that  \eqref{eq:biot_savart} and \eqref{eq:Rosenhead-Moor} are the restrictions at the curve of the entire fluid velocity field over the 3D domain
\begin{align}
    \label{eq:RMAll}
       U^{\on{RM}}_{\gamma}(  x)
    \coloneqq \frac{\Gamma}{4 \pi}  \oint  \frac{\gamma^{\prime}(\tilde s) \times(  x-\gamma(\tilde s))}{\sqrt{e^{-\nicefrac{3}{2}}a^{2}+\left|   x-\gamma(\tilde s)\right|^{2}}^{3}} d\tilde s,\quad  x\in\RR^3
\end{align}
That is, \(   V^{\on{RM}}_\gamma(s) =    U^{\on{RM}}_\gamma(\gamma(s)).\)
\end{example}

Now, we translate the dynamical system \eqref{eq:ExplicitDynamics} into an evolution equation for a time-dependent complex level function \(\psi\).
First, we note that the evolution of \(\psi\) around the zeros is given as the \emph{transport equation} along some vector field, which is formally the following lemma.
\begin{lemma}\label{lem:transport_equation}
    For any time-dependent complex level set function \(\psi\) with 
    the regularity assumptions in \autoref{sec:RepresentationsForFilamentDynamics},
    there exists a neighborhood \(U\subset M\) of the zero set of \(\psi\) and a vector field \(   v\colon U\rightarrow\RR^3\) such that
    \begin{align}
        \label{eq:AdvectionInU}
        {\partial\psi\over\partial t} +    v\cdot\nabla\psi = 0\quad\text{in \(U\)}.
    \end{align}
\end{lemma}
\begin{proof}
    %Using \autoref{sec:Regularity} 
    By the regularity assumptions in \autoref{sec:RepresentationsForFilamentDynamics},
    there exists a neighborhood \(U\) of the zero set of \(\psi\) where the 3D-to-2D linear map \(d\psi\) is full-rank and thus surjective.  
    Hence at every point \(  x\in U\) and for any value of \(\nicefrac{\partial\psi}{\partial t}(  x)\), there exists a vector \(   v_{  x}\in\RR^3\) such that \(d\psi\vert_  x(   v_{  x}) = -\nicefrac{\partial\psi}{\partial t}(  x)\).
    With this construction, we obtain a vector field \(   v\colon U\rightarrow\RR^3\) satisfying \(\nicefrac{\partial\psi}{\partial t} +    v\cdot\nabla\psi=0\) in \(U\).
\end{proof}

Hence the representation of the dynamics for \(\psi\) can be encapsulated into a \emph{vector field}.
Observe that the zero level curve \(\gamma\) of \(\psi\) evolves by \(\nicefrac{\partial\gamma}{\partial t} =    v\vert_\gamma\) under the transport equation \eqref{eq:AdvectionInU}, and that \(\operatorname{ker} (d\psi)\) at each point on \(\gamma\) is spanned by the tangent \(\gamma'\) as \(\gamma'\) lies on the tangent spaces of both \(\{\textup{Im} \psi =0 \}\) and  \(\{\textup{Re} \psi =0 \}\).  By matching these induced curve dynamics in Lemma \ref{lem:transport_equation} with \eqref{eq:ExplicitDynamics} we conclude:
\begin{theorem}
    \label{thm:VelocityCondition}
    The zero level curve \(\gamma\) of \(\psi\) evolves according to \eqref{eq:ExplicitDynamics} if and only if \(\psi\) satisfies \eqref{eq:AdvectionInU} for some vector field \(   v\) that agrees with the curve velocity at the curve:
    \begin{align}
        \label{eq:VelocityCondition}
           v(\gamma(s)) =    V_\gamma(s) + f(s)\gamma'(s)
    \end{align}
    where \(f\gamma'\) is the tangent vector \(\gamma'\) multiplied by an arbitrary scalar function \(f\).
    The degrees of freedom of the dynamics for \(\psi\) are the degrees of freedom for choosing \(   v\) with the condition \eqref{eq:VelocityCondition}.
\end{theorem}

Essentially, the only velocities that really matter for the evolution of the curves are the velocities located {\em on the zero level set}. Furthermore, the locations of the curves will not change if we slide them around their tangent direction (like spinning a circle around its axis of symmetry), so we only need to pin down their {\em normal and bi-normal components}. So we have a huge number of velocity variables (3 for each point in the 3D domain) with very few constraints (2 for each point on the 1D curves). This under-determined system gives us a {\em redundancy} in possible velocity fields for curve dynamics, which is largely unexplored by previous work.

Let us apply \autoref{thm:VelocityCondition} to vortex filament dynamics (\cf\@ \autoref{ex:BiotSavart}): Plugging in \eqref{eq:RMAll} for $   v$ gives us the following dynamics for $\psi$:
\begin{align}
    \label{eq:NaiveAdvection}
    {\partial\over\partial t}\psi +   v\cdot\nabla\psi = 0,\quad   v(  x) =    U^{\on{RM}}_\gamma(  x).
\end{align}
This is the most straightforward way to do it: simply advect the level set $\psi$ in the exact same way as the rest of the fluid. However, we know that \(   v=   U^{\on{RM}}_\gamma\) is an extremely sensitive function to deal with numerically --- it tends to infinity near $\gamma$, has unbounded derivatives, and rapidly changes direction in very tight swirls. Small errors in $\gamma$ inevitably create huge errors in velocity, making simulations unstable, as demonstrated in our accompanying video. Fortunately, according to \autoref{thm:VelocityCondition} we now know that there are infinitely many velocity fields that will all theoretically give us the same filament motions; our mission in the next section is to swap out this unstable Biot-Savart velocity field for one that is much more numerically robust.

\begin{remark}
When a curve has a self-intersection, that is, there is a point $p=\gamma(s)=\gamma(s')$ for some $s\neq s'$, the differential $d\psi_p$ fails to be surjective. Consequently, \autoref{lem:transport_equation} and hence \autoref{thm:VelocityCondition} do not apply. Intuitively, the time evolution of a self-intersecting curve is not described by the advection of $\psi$ along a single vector field as the overlapping points $\gamma(s)$ and $\gamma(s')$ may have different velocities.

In the numerical scheme we present in \autoref{sec:algorithms}, we handle this issue by introducing a mechanism that automatically resolves self-intersections as soon as they appear.
\end{remark}

\section{Untwisted Clebsch variables and non-swirling dynamics}
\label{sec:untwisted}


\subsection{Untwisted Clebsch variables}
\label{sec:UntwistedClebschVariables}

Like the common codimension-1 real-valued level set methods, the implicit representation benefits from the regularity of the level set function.  There, a level set function is well-conditioned if the magnitude of the gradient is close to one.  For that reason, the level set function is typically initialized as the \emph{signed distance function}, and this property is typically maintained as the level set is evolved (called \emph{re-distancing}). 

For our codimension-2 complex level set representation, we shall also characterize a set of desirable qualities of the complex level set function \(\psi\colon M\rightarrow\CC\). Due to the higher codimension, the discussion involves the notion of \emph{twist} from the mathematical  \emph{ribbon theory}. Finally, we describe a concrete construction of \(\psi\) that will be used for initialization and re-distancing. 

\subsubsection{Conditioning of a complex level set function}
We want \(\psi\) to be continuous everywhere and non-zero outside the curves \(\gamma=\{\psi=0\}\). 

Near the curves \(\gamma\), we want the differential \(d\psi\) to be well-conditioned:
if we only consider the function restricted to the plane spanned by the curve normal and bi-normal, \(d\psi\vert_{\gamma\bot}\) is close to an isometry; that is \(d\psi\vert_{\gamma^\bot}\colon\gamma^\bot\to\CC\cong\RR^2\) has singular values \(\approx(1,1)\).
Notice, however, that even when \(d\psi\vert_{\gamma^\bot}\) is well-behaved on each normal plane, the level set function \(\psi\) can still exhibit significant {\em shearing} if the complex phase varies significantly along the tangent direction.

This variation of complex phase along the curve's tangent can be seen more intuitively via the geometry of the surface \(S_{\psi}\) formed by the level sets of the positive real part of \(\psi\) (\(S_{\psi} = \{\Re(\psi)< 0, \Im(\psi)=0\} = \{\arg(\psi)= \pi\}\)).  Note that the boundary of the surface \(S_{\psi}\) is \(\gamma\), as illustrated in \autoref{fig:ribbon}. The tangential variation of the phase \(\arg(\psi)\) is embodied by the twist of \(S_{\psi}\) at its boundary \(\gamma\).

For the specific method presented in this chapter, we take a simple construction of \(\psi\) that is known to have a smooth \(S_\psi\) with little twisting at \(\gamma\). We first lay out a few geometric and topological properties about the twist. Readers who are only concerned about the simulation method can skip ahead to the last part of this subsection, \emph{Solid-angle distance function}.

\subsubsection{The twist of a complex level set function}
To discuss the twist precisely, we first review the standard definition of twist for a ribbon of a curve. Then, we apply this definition to the ribbon induced by $\psi$, which is tangent to the surface \(S_{\psi}\).

\begin{definition}[Ribbon and twist]\label{def:ribbon_and_twist}
   Let $U$ be a unit normal vector field over a collection of non self-intersecting space curves $\gamma$. A \emph{ribbon} $R$ of $\gamma$ is defined as another collection of curves $R=\gamma+\epsilon U$ where $\epsilon>0$ is taken sufficiently small so that $R$ does not intersect with itself and $\gamma$.

   The \emph{twist} of the pair \((\gamma,R)\) at each point on \(\gamma\) is given by
    \begin{align}
        \omega \coloneqq U'\cdot(\gamma'\times U)
    \end{align}
    where \((\cdot)'\) is the derivative along curve's tangent relative to arclength.

    Then the \emph{total twist} of  \((\gamma,R)\) is defined as
\begin{align}
    \operatorname{Tw}(\gamma,R)\coloneqq {1\over 2\pi}\oint_{\gamma}\omega.
\end{align}
\end{definition} 

We now consider the ribbon induced by $\psi$.
\begin{definition}[\(\psi\)-induced framing]
    Each complex level set function \(\psi\colon M\rightarrow\CC\) for \(\gamma\) gives rise to a normal vector field \(U_\psi\colon\gamma\to\RR^3\), \(|U_\psi|=1\), \(\langle \gamma',U_\psi\rangle = 0\) that points to the direction where \(\psi\) is real and positive.  Explicitly, at each point on \(\gamma\),
    \begin{align}
        U_\psi \coloneqq \frac{(d\psi\vert_{\gamma^\bot})^{-1}(-1)}
        {|(d\psi\vert_{\gamma^\bot})^{-1}(-1)|}.
    \end{align}
    This defines the ribbon $R_\psi\coloneqq \gamma+\epsilon U_\psi$ with a small enough $\epsilon>0$ as in \autoref{def:ribbon_and_twist}. 
    Using this, we define the twist $\omega_\psi$ and the total twist $\on{Tw}(\psi)$ as those of the ribbon $(\gamma,R_\psi)$.
\end{definition} 
% This normal vector field \(U_\psi\) makes \(\gamma\) a \emph{framed curve}.

% \begin{definition}[Ribbon and twist]


%     The \emph{twist} of \(\psi\) at each point on \(\gamma\) is given by
%     \begin{align}
%         \omega \coloneqq U_\psi'\cdot(\gamma'\times U_\psi).
%     \end{align}
%     Here, \((\cdot)'\) is the derivative along curve's tangent relative to arclength.
% \end{definition} 



When \(d\psi\vert_{\gamma^\bot}\colon \gamma^\bot\to\CC\) is conformal (which is true, for the solid-angle distance function we will introduce in this subsection, since it is isometric), the twist directly relates to the derivative of \(\psi\) along the curve:
 \begin{align}
     \label{eq:TwistDPsi}
     \omega_\psi = -(d\theta)(\gamma')= -\Re\left({(d\psi)(\gamma')\over i\psi}\right),\quad\theta = \arg\psi,
 \end{align}
in a small neighborhood of \(\gamma\).

\begin{figure}[htbp]
 \centering
     \begin{minipage}[b]{.8\linewidth}
         \centering
         \includegraphics
         [width=0.98\linewidth,trim={6cm 7cm 6cm 6cm},clip]         {{images_hidden_degrees/illustrations/level_surfaces_and_ribbons/Active_Render.0002.48_color_adjusted}.png}
     \end{minipage}
    \caption{The surface \(S_\psi\) (translucent blue) and the 
    %\cCW{twisted} 
    ribbon \(U_\psi\) (opaque cyan).
    \label{fig:ribbon}
    }
\end{figure}

Now, a natural consideration for designing \(\psi\) is to minimize the twist.  For example, one may try to construct \emph{Bishop's parallel frame} \cite{bergou2008elastic}, which has no twist.  However, one would find the construction \emph{impossible} for a closed curve and a union of closed curves in general.  This is due to the following theorem that the total twist $\on{Tw}(\psi)$ is fixed by the curve geometry \(\gamma\in\cX\) or the equivalent class of level functions \([\psi]\in\cF\), and hence non-zero unless $\on{Tw}(\gamma)=0$.

\begin{theorem}[Invariance of total twist]\label{th:total_twist_invariance}
    Suppose \(M=\RR^3\).  Let \([\psi]\in\cF\). Then any representative \(\psi\in[\psi]\) has the same total twist \(\operatorname{Tw}(\psi)\).
\end{theorem}
To prove the theorem, we use an auxiliary setting.

\begin{definition}[Seifert surface and Seifert framing]\label{def:Seifert_surface_framing}
    A Seifert surface $\Sigma$ of $\gamma$ is an oriented and compact surface such that $\partial\Sigma=\gamma$. \footnote{Note that a standard definition of Seifert surface additionally requires connectedness \cite{Seifert1935seifertsurface, musasugi2008knottheory}, but we drop this condition for our purposes.}
    A ribbon $R$ of $\gamma$ is called a Seifert framing if there is a Seifert surface $\Sigma$ such that $R$ is tangent to $\Sigma$ on $\gamma=\partial \Sigma$.
\end{definition}


\begin{definition}[Ribbon helicity]\label{def:ribbon_helicity}
    Let $\gamma_i$ be each connected component of $\gamma$ and $R_i$ be a ribbon of $\gamma_i$ Then the \emph{ribbon helicity} of $\gamma_i$ is 
    \begin{align}
        H(\gamma_i,R_i)= \on{Link}(\gamma_i,R_i)+\sum_{j\neq i} \on{Link}(\gamma_i,\gamma_j)
    \end{align}
    where $\on{Link}$ is the linking number between non-intersecting curves.
\end{definition}
The linking number in \autoref{def:ribbon_helicity} measures how much two curves are linked. For details we refer to \cite{ArnoldKhesin2021}.

When the ribbon is a Seifert framing, the helicity of the ribbon vanishes (\cite[Theorem 3]{deRicca2021zeroHelicity}):
\begin{lemma}[Zero ribbon helicity of Seifert framing]\label{lem:zero_helicity_Seifert}
    If the ribbon $R$ of a collection of curves $\gamma$ is a Seifert framing, we have
    \begin{align}
        H(\gamma_i,R_i)=0, \quad \forall i
    \end{align}
    where $(\gamma_i,R_i)$ are each connected component of $\gamma$ and its associated component of $R$.
\end{lemma}

Another key ingredient is the C\u{a}lug\u{a}reanu theorem:
\begin{lemma}[C\u{a}lug\u{a}reanu]\label{lem:Calugareanu}
For a single curve $\gamma$ and a ribbon $R$, we have 
\begin{align}
    \on{Link}(\gamma, R)= \on{Tw}(\gamma,R)+\on{Wr}(\gamma)
\end{align}
where $\on{Wr}(\gamma)$ is the writh of $\gamma$ (for details, see \cite{ArnoldKhesin2021}).
\end{lemma}

We are now ready to prove the theorem.
\begin{proof}[Proof of \autoref{th:total_twist_invariance}]
    Without loss of generality, we can assume a constant phase-shift to \(\psi\) so that \(S_{\psi}\) is a Seifert surface of \(\gamma\).\footnote{The connectedness of $S_\psi$ is not necessarily guaranteed, but \autoref{def:Seifert_surface_framing} and \autoref{lem:zero_helicity_Seifert} do not require connectedness.}
    Hence for the ribbon $R_\psi$, we have   $H(\gamma_i, R_{\psi,i})=0$ and namely \( \on{Link}(\gamma_i,R_{\psi,i})=-\sum_{j\neq i} \on{Link}(\gamma_i,\gamma_j)\) for each $i$, due to \autoref{lem:zero_helicity_Seifert}. 
    Applying the C\u{a}lug\u{a}reanu theorem (\autoref{lem:Calugareanu}), we have
    \begin{align}
        \on{Tw}(\psi)
        = \on{Tw}(\gamma, R_\psi)
        &= \sum_i \on{Tw}(\gamma_i,R_{\psi,i})\\
        &= \sum_i - \on{Wr}(\gamma_i) + \on{Link}(\gamma_i, R_{\psi,i} )\\
        & = - \sum_i\left( \on{Wr}(\gamma_i) +\sum_{j\neq i} \on{Link}(\gamma_i,\gamma_j)\right).
    \end{align}
    This show that $\on{Tw}(\psi)$ depends only on $\gamma$.
    % We finally apply the C\u{a}lug\u{a}reanu theorem, stating that \(\on{Link}(\gamma_i, R_i)= \on{Tw}(\gamma_i,R_i)+\on{Wr}(\gamma_i)\) for any ribbon $R_i$, where $\on{Wr}$ is the writh of $\gamma_i$ (for details, see \cite{ArnoldKhesin2021}). 
\end{proof}

% \begin{proof}
%     In the special case that \(\gamma\) is a single connected curve in \(\RR^3\), the statement is the result of the C\u{a}lug\u{a}reanu Theorem: The total twist of the frame \(\operatorname{Tw}(\psi)\) and the total writhe of the curve \(\operatorname{Wr}(\gamma)\) must add up to the linking number between the two boundaries of the ribbon swept out by the frame \(U_\psi\).  Here, we note that any closed curve admits a Seifert surface i.e. a compact, connected, and oriented surface spanning the curve 
%     \cite{
%         Seifert1935seifertsurface,
%         musasugi2008knottheory}.
%     Without loss of generality, we can assume a constant phase-shift to \(\psi\) so that \(S_{\psi}\) is a Seifert surface of \(\gamma\).
%      Since the ribbon \(U_\psi\) lies in a Seifert surface, the linking must be zero.  Therefore, \(\operatorname{Tw}(\psi) = -\operatorname{Wr}(\gamma)\), which only depends on \(\gamma\). See \autoref{fig:ribbon} for a visualization of \(S_\psi\) and \(U_\psi\). For the definitions of quantities for curves such as writhe, twist, or linking numbers, we refer to 
%     \cite{ArnoldKhesin2021}.
    
%     In the general case, observe that under any change of representative \(\tilde\psi = \varphi\psi\) in \([\psi]\) with a nowhere-vanishing \(\varphi\colon M\rightarrow\CC\setminus \{0\}\), the twist transforms according to \(\tilde\omega = \omega - d_{\gamma'}(\arg(\varphi))\).  Since \(\varphi\) has no zero in a simply connected \(M\), the angle \(\arg(\varphi)\) can be defined as a real-valued function globally over \(M\) rather than a \(2\pi\)-ambiguous angle-valued function. In particular, Stokes' theorem applies and the integral along the closed curve vanishes: \(\oint_{\gamma}d(\arg(\varphi)) = 0\). So \(\oint_\gamma\tilde\omega = \oint_\gamma\omega\).
% \end{proof}

Note that even though the total twist is fixed among possible choices of $\psi$, strong {\em local} twist can still be present.  
 
\subsubsection{Solid-angle distance function}
% \label{sec:SolidAngleDistanceFunction}

As explained above, we desire a complex level set function $\psi$ that is numerically well-conditioned: it should be nearly isometric near the curve $\gamma$ and non-zero outside $\gamma$. In order to isolate these properties and enforce them explicitly, we model $\psi$ with a complex wave function:
\begin{align}
    \psi(  x) = r(  x) e^{i\theta(  x)}.
\end{align}
Similar to \emph{signed distance functions} for the codimension-1 level set functions, we set
\begin{align}
    |\psi(  x)| = r(  x) \coloneqq \operatorname{dist}(\gamma,  x),\quad  x\in M\subset\RR^3.
\end{align}
This setup ensures that the value \(\psi(  x)\) is non-zero for \(  x\notin\gamma\).
What remains is a choice for the complex phase \(\theta\colon M\setminus\gamma\to\SS^1 = \RR_{\text{mod }2\pi}\),
which we set to
the \emph{half solid-angle subtended by \(\gamma\)}:
\begin{align}
    \theta(  x) \coloneqq {1\over 2}\operatorname{SolidAngle}(\gamma;  x)\mod 2\pi
\end{align}
where \(\operatorname{SolidAngle}(\gamma;  x)\in\RR_{\text{mod $4\pi$}}\) is a dimensionless quantity given by the signed spherical area \footnote{Our definition of the solid angle is modulated by \(4\pi\) as the spherical polygon of projected curves may cover the sphere multiple times.} enclosed by the projection of \(\gamma\) on the unit sphere centered at \(  x\) i.e. \(\operatorname{Proj}\gamma_{  x}(s):=(\gamma(s)-  x)/|\gamma(s)-  x| \).
Combining these two, we now have the \emph{solid-angle distance function},
\begin{align}\label{eq:solid_angle_distance}
    \psi(x)=\on{dist}(\gamma,x) e^{i \on{SolidAngle}(\gamma,x)/2}
\end{align}
illustrated in \autoref{fig:projection_onto_sphere}.


% \begin{figure}[htb]
%     \centering
%     \hspace{20pt}
%     \setlength\unitlength{1pt}
%     \begin{picture}(110,110)
%         \put(0,0){\includegraphics[width=0.6\columnwidth]{images_hidden_degrees/illustrations/projected_curve/SolidAngle.png}}
%         \put(105,65){\small\(\gamma\)}
%         \put(15,100){\small\(  x\)}
%         \put(-10,80){\small\(\operatorname{SolidAngle}(\gamma;  x)\)}
%         \put(25,105){\small\(\operatorname{dist(\gamma,  x)}\)}
%         \put(-50,40){\small\(\psi(  x) = \operatorname{dist}(\gamma,  x)e^{{i\over 2}\operatorname{SolidAngle}(\gamma;  x)}\)}
%     \end{picture}

    
%     \caption{The solid-angle distance function \(\psi\) for a space curve \(\gamma\) is constructed by the distance and the angle subtended by the curve.}
%     \label{fig:projection_onto_sphere}
% \end{figure}

\begin{figure}
\centering
\begin{overpic}[width=0.7\textwidth]{images_hidden_degrees/illustrations/projected_curve/SolidAngle.png}
% \put(x,y). 0,0 is left bottom. 
        \put(75,40){\(\gamma\)}
        \put(10,70){\(  x\)}
        \put(10,55){\(\operatorname{SolidAngle}(\gamma;  x)\)}
        \put(25,70){\(\operatorname{dist(\gamma,  x)}\)}
        \put(-10,40){\(\psi(  x) = \operatorname{dist}(\gamma,  x)e^{i\operatorname{SolidAngle}(\gamma;  x)/2}\)}
        %\graphpaper(0,0)(110,110)
\end{overpic}
        \caption{The solid-angle distance function \(\psi\) for a space curve \(\gamma\) is constructed by the distance and the angle subtended by the curve.}
    \label{fig:projection_onto_sphere}
\end{figure}

The solid-angle distance function meets our desired conditions for the codimension-2 level set representation.  On each normal plane \(\theta\vert_{\gamma^\bot}\) is asymptotically the 2D angle function about the zero \(\gamma\), so that \(d\psi\vert_{\gamma^\bot}\colon \gamma^\bot\to\CC\) is close to an isometry.  

Moreover, as studied by \cite{Binysh_2018solidangle}, 
 the surface \(S_\psi=\{\theta=0\}\) features little twist at \(\gamma\). \autoref{fig:solid_angle} illustrates this concept.
 
 \begin{figure}[htbp]
    \centering
    \includegraphics[width=1.0\columnwidth]{images_hidden_degrees/illustrations/solid_angle/solidAngle_with_bar.png}
    \caption{Plotting $\theta$ on a 2D plane which intersects a circular vortex ring at two points (left). The color indicates the value of the $\theta$, and the black lines are its level curves. The curves meet where the filament intersects the plane. Zooming into the white box (right) shows evenly-spaced curves closer to the filament, where $e^{i\theta}$ resembles the complex plane $\CC$.}
    \label{fig:solid_angle}
\end{figure}

\begin{remark}[Properties of untwisted Clebsch variable]
    We note that our choice of the Clebsch variables \(\psi\) is an instance of many reasonable ones that exploit degrees of freedom in curve representations rather than the optimal one for a specific dynamical system such as vortex filaments.
    
Nevertheless, our \(\psi\) has a number of desirable properties. For example $\psi$ is uniquely computed from given curves with an explicit formula, and is locally little twisted with  \(d\psi\vert_{\gamma^\bot}\) being locally nearly-isometry, as we saw in this section. In addition, $\arg \psi$ is harmonic and the zeros of the real and the imaginary parts intersect orthogonally, which ensures the zeros to be always codimension-2 and makes numerically extracting zeros robust. 
\end{remark}


\begin{remark}[Open curves]
Our construction of \(\psi\) is not limited to closed curves but is also valid for open curves 
with end points located on the boundary or obstacles subject to the following integrability condition.  The orientation of each curve assigns a positive or negative signature for its two end points.  The collection of curves is said to be integrable if all ends of the curves lie on the boundary and each connected component of the boundary contains an equal number of positive and negative ends.  Note that vortex filaments must be integrable since this integrability condition is precisely the integrability condition for curl: The vorticity 2-form is the exterior derivative of a velocity 1-form if and only if it is closed (divergence-free) and its restriction to every boundary component has zero total flux. In this case, the integrability of curl ensures that the endpoints of open filaments landing on boundaries must give equal positive and negative ends per boundary component. Therefore, it is possible to pair the endpoints along the boundary and complete the filaments as closed curves, from which we know how to construct \(\psi\). 
\end{remark}

\subsection{Non-swirling dynamics}
\label{sec:NonSwirlingDynamics}
In \autoref{sec:UntwistedClebschVariables}, we leveraged the degrees of freedom in the complex level set function \(\psi\) to design a sufficiently regular implicit representation.  
Here, we exploit similar degrees of freedom to construct evolution equations that produce theoretically equivalent dynamics but are more numerically robust. 

Many dynamical systems for curves already come with a known physical evolution equation.  For example, the vortex filament dynamics can be simulated with \eqref{eq:NaiveAdvection}, \ie\@ by advecting \(\psi\) using the Biot--Savart or the Rosenhead--Moore flow (\autoref{ex:BiotSavart}).   Hence, redesigning the equation may seem unnecessary.  However, when it comes to numerically advancing the variables, the highly oscillatory or discontinuous nature of the Biot--Savart and the Rosenhead--Moore flows near the vortex core (as illustrated in \autoref{fig:different_velocity_fields}) can cause significant interpolation error.  We point out that these errors are avoidable by redesigning the flow of the advection.

\begin{figure}[htbp]
\centering
     \begin{minipage}[b]{1.0\linewidth}
          \centering
          \includegraphics
          [width=0.95\linewidth,trim={0cm 4cm 0cm 4cm},clip]
          {{images_hidden_degrees/illustrations/vector_fields/bs_t2}.png}
      \end{minipage}
      \begin{minipage}[b]{1.0\linewidth}
          \centering
          \includegraphics
          [width=0.95\linewidth,trim={0cm 4cm 0cm 4cm},clip]
          {{images_hidden_degrees/illustrations/vector_fields/ne_t2}.png}
      \end{minipage}
      \begin{minipage}[b]{1.0\linewidth}
          \centering
          \includegraphics
          [width=0.95\linewidth,trim={0cm 4cm 0cm 4cm},clip]
          {{images_hidden_degrees/illustrations/vector_fields/good_t2}.png}
      \end{minipage}
     \caption{Different vector fields for two linked rings. The original Rosenhead--Moore model (top), the nearest point velocity field (middle), and a smooth  weighted average field (bottom). While these three velocity fields coincide on the filaments, they differ significantly outside the filaments.
     }
     \label{fig:different_velocity_fields}
 \end{figure}

We consider dynamical systems as discussed in \autoref{sec:RepresentationsForFilamentDynamics}.  
Suppose the evolution of the curve is given by \({\partial\gamma\over\partial t} =    V_{\gamma}\) where \(   V_{\gamma}\) is the velocity field defined on the curve.  The evolution for \(\psi\) must be an advection by an extension \(   v\) of the velocity field \(   V_\gamma\) in a neighborhood of the curve  (\cf\@ \eqref{eq:AdvectionInU}).   A straightforward construction  of \(   v\) is a constant extrapolation.  That is, \(   v(  x)\) is set to \(   V_\gamma\) at the closest point on \(\gamma\) from \(  x\). This extrapolation is, however, singular where closest points are not unique as in the middle row of \autoref{fig:different_velocity_fields}.

To gain continuity without changing the velocity on the filaments, we 
smooth away these singularities by taking the weighted average of the filament velocity as 
\begin{align}\label{eq:non-swirling_velocity}
        v(  x)\coloneqq \frac{1}{\cN(  x)}\oint_\gamma    V_\gamma(\gamma(s)) w(  x,\gamma(s))\,  ds.
\end{align}
Here, \(\cN(  x)\) is the normalization factor
\begin{align} \label{eq:normalization_factor}
     \cN(  x) = \oint_\gamma w(  x,\gamma(s)) ds,
\end{align}
and \(w\) is some weight function 
that applies less smoothing as it gets closer to the filament, i.e. 
{\(w(  x, \gamma(s)) / \cN(  x) \rightarrow \delta(\gamma^{-1}(  x)-s)\)} as  {\( \textup{dist}(  x,\gamma) \rightarrow 0\)}.
For example, we observed that a Gaussian function with distance-dependent variance works stably:
\begin{align} \label{eq:weight_function}
    w(  x,s)=\exp\left(-\frac{|  x-\gamma(s)|^2}{ \sigma^2 \textup{dist}(  x,\gamma)}\right),\quad\sigma\text{ some constant}.
\end{align}

To accelerate computation for \(   v\), we can further multiply the integrand of \autoref{eq:non-swirling_velocity} by a smooth cutoff function which equals to \(1\) near \(\gamma\) and \(0\) far away from \(\gamma\). Then \(   v\) is non-vanishing only near \(\gamma\).  By applying this smooth cutoff, we only need to evaluate \(   v\) in a narrow band close to the filament.

Note that \(   v\) is in general not divergence-free. The velocitiy near the curve is determined according to the curve velocities so that the motion of the zeros of \(\psi\) emulates the motion of the curves. Imposing an additional constraint like incompressibility to the velocitity field may trade off the fidelity to the original curve motion or the numerical smoothness of the surrounding vector field.


%!TEX root = VortexUntwisted_paper.tex
\section{Algorithms}
\label{sec:algorithms}
In this section, we describe an algorithm for simulating filament dynamics.

Throughout the simulation, we maintain a complex level set function \(\psi\). The main algorithm computes the transport equation of \(\psi\) along a velocity field \(   v\) in a neighborhood of the zeros of \(\psi\).  This main algorithm is accompanied by a few subroutines for evaluating the velocity and redistancing: one subroutine extracts the zero set \(\gamma\) of \(\psi\); another subroutine constructs the solid-angle distance function \(\psi\) from \(\gamma\) \eqref{eq:solid_angle_distance}; a third subroutine evaluates the filament motion \(   V_\gamma\) using \(\gamma\); and the last subroutine extends \(   V_\gamma\) to a velocity field \(   v\) in a neighborhood of \(\gamma\) (\autoref{sec:NonSwirlingDynamics}).

We store the level set function \(\psi\) and the velocity \(   v\) on a 3D lattice and discretize the curves \(\gamma\) as oriented collections of line segments. 

\begin{algorithm}[htbp]
% \begin{algorithmic}[1]
%      \Require{Initial filament \(\gamma\in\cX\); %timestep \(\Deltait t>0\);
%      }
%      \State{\(\psi\gets\) construct \(\psi\) from \(\gamma\);}
%      \Comment{\autoref{sec:ComputePhase}} 
% \end{algorithmic}
  \SetKw{BreakFor}{break for}
  \KwIn{Initial filament \(\gamma\in\cX\);}
  \(\psi\gets\) construct \(\psi\) from \(\gamma\); \Comment{Subroutine \ref{subroutine:ComputePhase}} 
  
  \While{simulating}
  {
    \(   v\gets\) evaluate \(   V_\gamma\) and extend it on grids near \(\gamma\);\Comment{Subroutine \ref{subroutine:extendVelocity}}

    \(\psi\gets\) advect \(\psi\) along \(   v\);\Comment{Subroutine \ref{subroutine:Advection}}

    \(\gamma\gets\) extract the zero set of \(\psi\); 
    \Comment{Algorithm \ref{alg:construct_curve}}
    
    \(\psi\gets\) construct \(\psi\) from \(\gamma\);
    \Comment{Subroutine \ref{subroutine:ComputePhase}}
  }
  \caption{The main time integration}
  \label{alg:time_integration} % \label has to be placed AFTER \caption to produce correct cross-references.
\end{algorithm}

\begin{algorithm}[htbp]
  \SetKw{BreakFor}{break for}
  \ForEach{ cell $c$ }
  {
     \ForEach{face $f$ in $c$}
     {

     Compute incidence \(n_f\in\{-1,0,1\}\); \Comment{Eq.~\eqref{eq:ArgumentPrinciple}}

     \If{$\lvert\vec{x}^{(k)}-\vec{x}^{(k-1)}\rvert < \epsilon$}
     {Find \(   p_f^\pm = \psi^{-1}(0)\in\RR^3\)  in \(f\); \Comment{bilinear interp.}}
    }

     Connect \(   p_f^-\) and \(    p_g^+ \) of some faces \(f, g\) in \(c\);
     
  }
   \caption{Extract the zero curve \(\gamma\) from \(\psi\)}
  \label{alg:construct_curve} % \label has to be placed AFTER \caption to produce correct cross-references.
\end{algorithm}

\subsection{Details of subroutines}

\begin{enumerate}
\item \textbf{Advection}
    \label{subroutine:Advection}
    
To advect \(\psi\) with a given flow \(   v\), one can adopt any Eulerian advection scheme.  In our implementation, we use the modified MacCormack method \cite{Selle:2008:MCM} with 4th order Runge--Kutta back-tracing.

\item \textbf{Construction of  \(\gamma\) from \(\psi\)}
\label{subroutine:gammaFromPsi}

After updating \(\psi\), we need to %construct new
update the
filaments \(\gamma\) by extracting \(\{\psi=0\}\). 
We summarize this subroutine in \autoref{alg:construct_curve}, which is adopted from \cite{weissmann2014smokerings}.

In our setting, 
each  vertex of \(\gamma\) lives on a face \(f\) of 
the volumetric grid.
We first evaluate for each face \(f\) the \(\{-1,0,+1\}\)-valued \emph{signed intersection} \(n_f\) with the zero curve of \(\psi\) using the \emph{argument principle}: If the vertices of a face \(f\) are \(i,j,k,\ell\) in an oriented order, then
\begin{align}
    \label{eq:ArgumentPrinciple}
    n_f = {1\over 2\pi}\left(\arg(\tfrac{\psi_j}{\psi_i}) + \arg(\tfrac{\psi_k}{\psi_j}) + \arg(\tfrac{\psi_\ell}{\psi_k}) + \arg(\tfrac{\psi_i}{\psi_\ell})\right)
\end{align}
using the principal branch \(-\pi<\arg(\cdot)\leq\pi\).
% \cAC{[the original writing \(\omega_f=\sum_{i=\{0,1,2,3\}}^{\textup{mod }4} \theta_f^{i+1} - \theta_f^i\) can be problematic.  The branching is subtle.]}
Geometrically, \eqref{eq:ArgumentPrinciple} describes how many times the quadrilateral \(\psi_i,\psi_j,\psi_k,\psi_\ell\in\CC\) winds around the origin.
For each face where \(n_f\neq 0\) we evaluate the more precise location of the zero using a bilinear interpolation. That is, we regard \(f=[0,1]^2\) by scaling and \(\psi\) is bilinearly interpolated as
\begin{align}
     \psi_f(x,y)
     := & (1-x)(1-y) \psi(0,0) + x(1-y)\psi(1,0)\nonumber\\ 
     &+(1-x)y\psi(0,1)+xy \psi(1,1),
\end{align}
and  \(\psi_f:f\rightarrow \CC\) has the inverse when \(n_f=\pm 1\).
The location \(\psi_f ^{-1}(0)\)  is a vertex \(   p_f\) of the curve \(\gamma\).
Finally, we build the edges of \(\gamma\) by running over the grid cells where we connect the pairs of zeros on the face with a consistent orientation. 
%\cCW{We are not discussing potentially ambiguous cases with 4 or 6 intersections per cell. Do we agree that this is too much low-level detail to be included here? Or is it not a problem at all?} \cSI{I think it's good to discuss it here. And actually I originally had an argument about it, so I added it back.}
%By a geometric consideration, it is easy to see that 
Each cube \(c\) may have up to two pairs of zeros with positive and negative  \(n_f\). A cube with two pairs of vertices has an ambiguity similarly to the marching cube algorithm \cite{Lorensen1987marchingcubes}. We resolve this ambiguity by connecting vertices arbitrarily in a way that preserves curve orientation.
%Within this work, w
We have not investigated higher order algorithms to connect curves more accurately at the sub-grid scale.

\item \textbf{Construction of \(\psi\) from \(\gamma\)}
\label{subroutine:ComputePhase}

Given \(\gamma\), we construct \(\psi\) as the solid-angle distance function \eqref{eq:solid_angle_distance}.  For each grid point \(  x\) near \(\gamma\), we evaluate \(|\psi(  x)|=\on{ dist}(\gamma,  x)\) as the distance to the closest polygon edge. 
We evaluate \(\arg(\psi(  x)) = \theta(  x) = {1\over 2}\on{SolidAngle}(\gamma;  x)\) by computing the signed area of a spherical polygon with vertices \(\{p_i\}_{i=1}^N\) that are the discrete points of \(\gamma\) projected onto the unit sphere centered at \(x\). To ensure differentiability and the local isometry properties of \(d\psi\vert_{\gamma\bot}\), the sign of $\theta$ must depend on curve orientation: we compute the {\em signed} area by introducing a pole \(Z\colon=(1,0,0)\):
\begin{align}
     \on{SignedArea}\left( \{p_i\}_{i=1}^N\right)=\sum_{i=1}^{N} \on{sign}\left((p_i \times p_{i+1})\cdot Z \right) \on{Area}(p_i,p_{i+1},Z)
\end{align}
where the {\em unsigned} area of each spherical triangle \(\{q_1,q_2,q_3\}\) is computed using a standard area formula \cite{Bevis1987ComputingTA,arvo1995applications}, 
\begin{align}
     \on{Area}(q_1,q_2,q_3)= -\pi+\sum_{i=1}^3 
     \arccos\left(
          \frac{
          \left(q_{i-1} \times q_{i}\right) \cdot\left(q_{i} \times q_{i+1}\right)} {\left\|q_{i-1} \times q_{i}\right\|\left\|q_{i} \times q_{i+1}\right\|}\right).
\end{align} 
These equations assume cyclic vertex 
indexing, so ${p_{N+1}=p_1}$, ${q_{3+1}=q_1}$, and ${q_{1-1}=q_3}$.

Alternatively, we can use the area formula for spherical polygons introduced in \autoref{ch:spherical_polygon}, which is robust even when the polygon $\{p_i\}_{i=0}^N$ has degenerate points.

\item\textbf{Extending velocity to grid}
\label{subroutine:extendVelocity}

In order to advect \(\psi\), we need to extend the velocity field defined on \(\gamma\) to the grid points near the curves. To produce non-swirling dynamics, we used a smooth average field (Equation \ref{eq:non-swirling_velocity}, \ref{eq:normalization_factor}, and \ref{eq:weight_function}) for the examples in this chapter, unless otherwise explained. In \autoref{eq:weight_function}, we observed that \(\sigma\) ranging from \(0.1\) to \(10\) times the grid size works stably without smoothing out detailed dynamics.
\end{enumerate}

% \subsubsection{Advection}
% \label{alg:Advection}
% To advect \(\psi\) with a given flow \(   v\), one can adopt any Eulerian advection scheme.  In our implementation, we use the modified MacCormack method \cite{Selle:2008:MCM} with 4th order Runge--Kutta back-tracing.

% \subsubsection{Construction of  \(\gamma\) from \(\psi\)}
% After updating \(\psi\), we need to %construct new
% update the
% filaments \(\gamma\) by extracting \(\{\psi=0\}\). 
% We summarize this subroutine in \autoref{alg:construct_curve}, which is adopted from \cite{weissmann2014smokerings}.

% In our setting, 
% each  vertex of \(\gamma\) lives on a face \(f\) of 
% the volumetric grid.
% We first evaluate for each face \(f\) the \(\{-1,0,+1\}\)-valued \emph{signed intersection} \(n_f\) with the zero curve of \(\psi\) using the \emph{argument principle}: If the vertices of a face \(f\) are \(i,j,k,\ell\) in an oriented order, then
% \begin{align}
%     \label{eq:ArgumentPrinciple}
%     n_f = {1\over 2\pi}\left(\arg(\tfrac{\psi_j}{\psi_i}) + \arg(\tfrac{\psi_k}{\psi_j}) + \arg(\tfrac{\psi_\ell}{\psi_k}) + \arg(\tfrac{\psi_i}{\psi_\ell})\right)
% \end{align}
% using the principal branch \(-\pi<\arg(\cdot)\leq\pi\).
% % \cAC{[the original writing \(\omega_f=\sum_{i=\{0,1,2,3\}}^{\textup{mod }4} \theta_f^{i+1} - \theta_f^i\) can be problematic.  The branching is subtle.]}
% Geometrically, \eqref{eq:ArgumentPrinciple} describes how many times the quadrilateral \(\psi_i,\psi_j,\psi_k,\psi_\ell\in\CC\) winds around the origin.
% For each face where \(n_f\neq 0\) we evaluate the more precise location of the zero using a bilinear interpolation. That is, we regard \(f=[0,1]^2\) by scaling and \(\psi\) is bilinearly interpolated as
% \begin{align}
%      \psi_f(x,y)
%      := & (1-x)(1-y) \psi(0,0) + x(1-y)\psi(1,0)\nonumber\\ 
%      &+(1-x)y\psi(0,1)+xy \psi(1,1),
% \end{align}
% and  \(\psi_f:f\rightarrow \CC\) has the inverse when \(n_f=\pm 1\).
% The location \(\psi_f ^{-1}(0)\)  is a vertex \(   p_f\) of the curve \(\gamma\).
% Finally, we build the edges of \(\gamma\) by running over the grid cells where we connect the pairs of zeros on the face with a consistent orientation. 
% %\cCW{We are not discussing potentially ambiguous cases with 4 or 6 intersections per cell. Do we agree that this is too much low-level detail to be included here? Or is it not a problem at all?} \cSI{I think it's good to discuss it here. And actually I originally had an argument about it, so I added it back.}
% %By a geometric consideration, it is easy to see that 
% Each cube \(c\) may have up to two pairs of zeros with positive and negative  \(n_f\). A cube with two pairs of vertices has an ambiguity similarly to the marching cube algorithm \cite{Lorensen1987marchingcubes}. We resolve this ambiguity by connecting vertices arbitrarily in a way that preserves curve orientation.
% %Within this work, w
% We have not investigated higher order algorithms to connect curves more accurately at the sub-grid scale.


% %To accelerate computation, we first evaluate local vorticity \(w_f\) on each face \(f\). Note that \(w_f\) is zero or \(\pm 2\pi\). When it is \(\pm 2\pi\), we further evaluate the precise location of \(p_f^{\pm}:=\{\psi=0\}\cap f\). We finally make an edge by connecting vertices \(p_f^+\) and \(p_f^-\) within cube \(c\). 


% % \begin{algorithm}[htb]
% %      \caption{Construct \(\gamma\) as \(\{ \psi=0\}\)}
% %      \label{alg:construct_curve}
% %      \begin{algorithmic}[1]
% %      \For{each cube \(c\)}
% %           \For{each face \(f\) in \(c\)}
% %           \State{Compute local vorticity \(\omega_f=\sum_{i=\{0,1,2,3\}}^{\textup{mod }4} \theta_f^{i+1} - \theta_f^i \)}  %\Comment{this is a comment}
% %                \If {\(\omega_f=\pm 2\pi\)}
% %                     \State{Find \(p_f^\pm = \psi^{-1}(0)\)  in \(f\)}
% %                \EndIf
% %
% %           \EndFor
% %           \State{Connect \(p_f^-\) and \( p_g^+ \) of some faces \(f, g\) in \(c\)}
% %      \EndFor
% %
% %
% %      \end{algorithmic}
% % \end{algorithm}

% \subsubsection{Construction of \(\psi\) from \(\gamma\)}
% \label{sec:ComputePhase}
% Given \(\gamma\), we construct \(\psi\) as described in \autoref{sec:SolidAngleDistanceFunction}.  For each grid point \(  x\) near \(\gamma\), we evaluate \(|\psi(  x)|=\on{ dist}(\gamma,  x)\) as the distance to the closest polygon edge. 
% We evaluate \(\arg(\psi(  x)) = \theta(  x) = {1\over 2}\on{SolidAngle}(\gamma;  x)\) by computing the signed area of a spherical polygon with vertices \(\{p_i\}_{i=1}^N\) that are the discrete points of \(\gamma\) projected onto the unit sphere centered at \(x\). To ensure differentiability and the local isometry properties of \(d\psi\vert_{\gamma\bot}\), the sign of $\theta$ must depend on curve orientation: we compute the {\em signed} area by introducing a pole \(Z\colon=(1,0,0)\):
% \begin{align}
%      \on{SignedArea}\left( \{p_i\}_{i=1}^N\right)=\sum_{i=1}^{N} \on{sign}\left((p_i \times p_{i+1})\cdot Z \right) \on{Area}(p_i,p_{i+1},Z)
% \end{align}
% where the {\em unsigned} area of each spherical triangle \(\{q_1,q_2,q_3\}\) is computed using a standard area formula \cite{Bevis1987ComputingTA,arvo1995applications}, 
% \begin{align}
%      \on{Area}(q_1,q_2,q_3)= -\pi+\sum_{i=1}^3 
%      \arccos\left(
%           \frac{
%           \left(q_{i-1} \times q_{i}\right) \cdot\left(q_{i} \times q_{i+1}\right)} {\left\|q_{i-1} \times q_{i}\right\|\left\|q_{i} \times q_{i+1}\right\|}\right).
% \end{align} 
% These equations assume cyclic vertex 
% indexing, so ${p_{N+1}=p_1}$, ${q_{3+1}=q_1}$, and ${q_{1-1}=q_3}$.

% \subsubsection{Extending velocity to grid}
% In order to advect \(\psi\), we need to extend the velocity field defined on \(\gamma\) to the grid points near the curves. To produce non-swirling dynamics, we used a smooth average field (Equation \ref{eq:non-swirling_velocity}, \ref{eq:normalization_factor}, and \ref{eq:weight_function}) for the examples in this chapter, unless otherwise explained. In \autoref{eq:weight_function}, we observed that \(\sigma\) ranging from \(0.1\) to \(10\) times the grid size works stably without smoothing out detailed dynamics.


\section{Applications}
\label{sec:applications}

This section discusses applications of our approach, specifically applied to vortex filament dynamics.
We implemented our algorithms on Houdini 18.5.759 and ran all simulations on a MacBook Pro (13-inch, 2020) with a 2.3 GHz Quad-Core Intel Core i7 processor. For an example implementation, see
\href{https://github.com/sdsgisd/ImplicitVortexFilaments}{https:\slash\slash{}github.com/sdsgisd\slash{}ImplicitVortexFilaments}


We use our algorithm to animate two ``leapfrogging'' vortex rings in \autoref{fig:leapfrog} and our accompanying video. We note that the system remains stable and highly symmetric even at the end of a long simulation with several high-speed ring interactions. \autoref{fig:oblique_rings} visualizes two vortex rings colliding with one another at right angles, reconnecting, and detaching into two new rings. We note that the final rings retain plenty of energy after the collision event, in contrast to Eulerian simulations of this phenomenon which tend to damp out over time. The visual detail in our simulations is also practically independent of grid resolution, as the motion of marker particles are described analytically by velocities (\autoref{eq:RMAll}), instead of a vector field stored on a coarse grid.

\begin{figure}[t]
    \centering{
       \includegraphics[width=\linewidth, trim={10cm 1cm 3cm 5cm},clip]{{images_hidden_degrees/leapfrog/leapfrog_particles_black_yellowblue.0170}.jpg}
       \caption{
           %Leapfrogging filaments and marker particles. From left to right are frames 1, 70, and 170. The initial setting is two circles located on the same plane whereas one circle has the half radius as the other.
           Two filaments leapfrog through one another, dragging marker particles into the shape of a mushroom cloud. The initial filament geometry consists of two co-planar vortex rings, one with half the radius of the other.
       }
       \label{fig:leapfrog}
    }
\end{figure}

\begin{figure}[t]
    \centering{
    \begin{subfigure}[b]{0.3\linewidth}
        \centering
        \scalebox{-1}[1]{
        \includegraphics[width=\linewidth,  trim={14cm 0cm 12cm 0cm},clip]
        {{images_hidden_degrees/obelique/obelique_curves.0001}.png}
        % {{images_hidden_degrees/obelique/obelique_curves.0001}.jpg}
        }
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.3\linewidth}
        \centering
        \scalebox{-1}[1]{
        \includegraphics[width=\linewidth, trim={14cm 0cm 13cm 0cm},clip]{{images_hidden_degrees/obelique/obelique_curves.0050}.png}
        % \includegraphics[width=\linewidth, trim={14cm 0cm 12cm 0cm},clip]{{images_hidden_degrees/obelique/obelique_curves.0050}.jpg}
        }
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.3\linewidth}
        \centering
        \scalebox{-1}[1]{
        \includegraphics[width=\linewidth, trim={14cm 0cm 14cm 0cm},clip]{{images_hidden_degrees/obelique/obelique_curves.0205}.png}
        % \includegraphics[width=\linewidth, trim={14cm 0cm 12cm 0cm},clip]{{images_hidden_degrees/obelique/obelique_curves.0205}.jpg}
        }
    \end{subfigure}
    }

    \centering{
    \begin{subfigure}[b]{0.3\linewidth}
        \centering
        \scalebox{-1}[1]{
        \includegraphics[width=\linewidth, trim={14cm 0cm 12cm 0cm},clip]
        % {{images_hidden_degrees/obelique/obelique_particles.0001}.png}
        {{images_hidden_degrees/obelique/obelique_particles.0001}.jpg}
        }
        % \caption*{Initial state}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.3\linewidth}
        \centering
        \scalebox{-1}[1]{
        % \includegraphics[width=\linewidth, trim={14cm 0cm 13cm 0cm},clip]{{images_hidden_degrees/obelique/obelique_particles.0050}.png}
        \includegraphics[width=\linewidth, trim={14cm 0cm 12cm 0cm},clip]{{images_hidden_degrees/obelique/obelique_particles.0050}.jpg}
        }
        % \caption*{Frame = 50}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.3\linewidth}
        \centering
        \scalebox{-1}[1]{
        % \includegraphics[width=\linewidth, trim={14cm 0cm 14cm 0cm},clip]{{images_hidden_degrees/obelique/obelique_particles.0205}.png}
        \includegraphics[width=\linewidth, trim={14cm 0cm 12cm 0cm},clip]{{images_hidden_degrees/obelique/obelique_particles.0205}.jpg}
        }
        % \caption*{Frame = 205}
    \end{subfigure}
    }

    \caption{Two smoke rings (left) colliding at orthogonal angles (middle)  and re-connecting (right). Rendered as filaments (top) and marker particles (bottom).
    The colliding rings leave swirly trails of smoke particles after their collision and reconnection.
    }
    \label{fig:oblique_rings}
\end{figure}

\autoref{fig:jet} illustrates a jet of smoke created by generating a new smoke ring at the left side of the domain every three time steps. 
For transporting smoke as a scalar field stored on grid points, we again used \autoref{eq:RMAll}.
This simulation shows the robustness of our topology changes: each re-connection event is the result of a curve extraction from a level set function, so there is no possibility of any unexpected edge cases, and no need for any geometric intersection code. The simple set-up creates a large variety of chaotic motions resulting from fast leap-frogging rings squeezing in between others and reconnecting filaments causing sudden changes in direction. When rings shrink smaller than the grid resolution, our algorithm deletes them (similar to codimension-1 level set methods). 

Lastly, \autoref{fig:torus} illustrates how our method can evolve intricate filament geometry, specifically the (5,8)-torus knot defined by
\begin{align}
    &\gamma(s)=\left(
        \left(\cos(q s)+2\right)\cos(ps),
        \left(\cos(q s)+2\right) \sin(ps),-sin(qs)
        \right) \notag
\end{align}
with \((p,q)=(5,8)\) and \(s\in [0,2\pi)\).


\begin{figure}[h]
\centering
    \begin{minipage}[b]{.95\linewidth}
        \centering
        \includegraphics[width=0.99\linewidth, trim={0cm 2cm 0cm 2cm},clip]{{images_hidden_degrees/jet/jet_curves.0189}.png}
        % \includegraphics[width=0.99\linewidth, trim={0cm 2cm 0cm 2cm},clip]{{images_hidden_degrees/jet/jet_curves.0189}.jpg}
    %   \caption{jet curves}
    \end{minipage}
    \begin{minipage}[b]{.95\linewidth}
        \centering
        \includegraphics[width=0.99\linewidth, trim={0cm 2cm 0cm 2cm},clip]{{images_hidden_degrees/jet/jet_smoke_intensity5.0189}.png}
        % \includegraphics[width=0.99\linewidth, trim={0cm 2cm 0cm 2cm},clip]{{images_hidden_degrees/jet/jet_smoke_intensity5.0189}.jpg}
    %   \caption{jet smoke}
    \end{minipage}
    \caption{A jet of smoke, rendered as raw filament geometry (top) and an advected smoke density function (bottom).
    }
    \label{fig:jet}
\end{figure}

\subsection{Influence of numerical parameters}

\autoref{fig:different_settings} demonstrates the importance of each step in our approach by selectively removing or modifying different algorithmic components and illustrating the consequences. First we illustrate what happens when we vary the free degrees of freedom in the velocity field $   v$ used to advect the filaments. In agreement with the discussion in \autoref{sec:RepresentationsForFilamentDynamics}, we see that setting $   v$ to the fluid velocity (based on Biot-Savart kernels) field causes the level set function $\psi$ to rapidly twist up and become unstable. Setting $   v$ to equal the velocity at the nearest point on the filament creates similar noise, presumably due to spatial discontinuities in the field. Compare these results to the smooth geometry generated by our velocity field at the bottom of \autoref{fig:different_settings}.

Lower down in the same figure, we illustrate the effect of varying the free degrees of freedom in our level set function $\psi$. Starting with an initially smooth $\psi$ and advecting it without any re-distancing or regularization works well at the beginning, but it eventually accumulates topological noise. To illustrate the impact of $\psi$'s free parameters on numerical stability and accuracy, the fourth row in \autoref{fig:different_settings} replaces our smooth choice of $\psi$ with one that is intentionally twisted by a phase shift of \(\Delta \theta(x):=0.05\  \operatorname{dist}(x,\gamma)\); the twisted $\psi$ causes high-frequency geometric noise and artificially shrinks the filaments. Again, we can compare these results to the smooth geometry generated by our un-twisted $\psi$ at the bottom of \autoref{fig:different_settings}.

\begin{figure}[htbp] 
\centering
    \begin{overpic}[width=0.35\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_BS.0020}.png}
        % {{images_hidden_degrees/twoRings/2rings_BS.0020}.jpg}
        % \put(20,0){ \small BS Frame = 20}
        % \put(40,0){\small Advection along the Biot-Savart velocity field}
    \end{overpic}
    % \hspace*{\fill}
    \begin{overpic}[width=0.35\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_BS.0045}.png}
        % {{images_hidden_degrees/twoRings/2rings_BS.0045}.jpg}
        % \put(20,0){ \small BS Frame = 45}
    \end{overpic}
    \vspace{-20pt}
    \caption*{\small Advection along the Biot-Savart velocity field}
    % \hspace*{\fill}
    % \vspace{-20pt}

    % \medskip
    \begin{overpic}[width=0.35\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_normal_extension.0020}.png}
        % {{images_hidden_degrees/twoRings/2rings_normal_extension.0020}.jpg}
        % \put(20,0){ \small NV Frame = 20}
        % \put(40,0){\small Advection along the nearest point velocity field}
    \end{overpic}
    % \hspace*{\fill}
    \begin{overpic}[width=0.35\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_normal_extension.0045}.png}
        % {{images_hidden_degrees/twoRings/2rings_normal_extension.0045}.jpg}
        % \put(20,0){ \small NV Frame = 45}
    \end{overpic}
    % \hspace*{\fill}
    \vspace{-30pt}
    \caption*{\small Advection along the nearest point velocity field}
    % \vspace{-20pt}

    \medskip
    \begin{overpic}[width=0.35\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_wo_redistance.0020}.png}
        % {{images_hidden_degrees/twoRings/2rings_wo_redistance.0020}.jpg}
        % \put(20,0){ \small NoRedist Frame = 20}
    \end{overpic}
    % \hspace*{\fill}
    \begin{overpic}[width=0.35\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_wo_redistance.0045}.png}
        % {{images_hidden_degrees/twoRings/2rings_wo_redistance.0045}.jpg}
        % \put(20,0){ \small NoRedist Frame = 45}
    \end{overpic}
    % \hspace*{\fill}
    \vspace{-30pt}
    \caption*{\small Without redistancing \(\psi\)}
    % \vspace{-20pt}

    \medskip
    \begin{overpic}[width=0.35\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_twisted.0020}.png}
        % {{images_hidden_degrees/twoRings/2rings_twisted.0020}.jpg}
        % \put(20,0){ \small Twisted Frame = 20}
    \end{overpic}
    % \hspace*{\fill}
    \begin{overpic}[width=0.35\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_twisted.0045}.png}
        % {{images_hidden_degrees/twoRings/2rings_twisted.0045}.jpg}
        % \put(20,0){ \small Twisted Frame = 45}
    \end{overpic}
    % \hspace*{\fill}
    \vspace{-30pt}
    \caption*{\small Artificially twisted \(\psi\)}
    % \vspace{-20pt}

    \medskip
    \begin{overpic}[width=0.4\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_ours.0020}.png}
        % {{images_hidden_degrees/twoRings/2rings_ours.0020}.jpg}
        % \put(20,0){ \small Ours Frame = 20}
        
    \end{overpic}
    % \hspace*{\fill}
    \begin{overpic}[width=0.4\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_ours.0045}.png}
        % {{images_hidden_degrees/twoRings/2rings_ours.0045}.jpg}
        % \put(20,0){ \small Ours Frame = 45}
   
    \end{overpic}
    % \hspace*{\fill}
        \vspace{-30pt}
    \caption*{\small Our proposed setting}
    % \vspace{-10pt}
    % \caption*{\small Frame = 20 \hspace{80pt} Frame = 45}
    \caption{Comparisons with different settings. Initial state is the 1st row of \autoref{fig:two_rings}. The left and right columns show states at frames 20 and 45.  
    } \label{fig:different_settings}
\end{figure}

\begin{figure}
    \centering{
        \begin{subfigure}[b]{0.48\linewidth}
            \centering
            \includegraphics[width=\linewidth,trim={6cm 2cm 6cm 2cm},clip]{{images_hidden_degrees/trefoil/trefoil_low_res_vert.0125}.png}
            % \includegraphics[width=\linewidth,trim={6cm 2cm 6cm 2cm},clip]{{images_hidden_degrees/trefoil/trefoil_low_res_vert.0125}.jpg}
        \end{subfigure}
        \hfill
        \begin{subfigure}[b]{0.4\linewidth}
            \centering
            \includegraphics[width=\linewidth,trim={6cm 2cm 6cm 2cm},clip]{{images_hidden_degrees/trefoil/trefoil_low_res.0125}.png}
            % \includegraphics[width=\linewidth,trim={6cm 2cm 6cm 2cm},clip]{{images_hidden_degrees/trefoil/trefoil_low_res.0125}.jpg}
        \end{subfigure}
        \vspace{-10pt}
        \caption*{Resolution \(50\times50\times60\)}
    }

    \centering{
    \begin{subfigure}[b]{0.4\linewidth}
        \centering
        \includegraphics[width=\linewidth,trim={6cm 2cm 6cm 2cm},clip]{{images_hidden_degrees/trefoil/trefoil_normal_res_vert.0125}.png}
        % \includegraphics[width=\linewidth,trim={6cm 2cm 6cm 2cm},clip]{{images_hidden_degrees/trefoil/trefoil_normal_res_vert.0125}.jpg}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.4\linewidth}
        \centering
        \includegraphics[width=\linewidth,trim={6cm 2cm 6cm 2cm},clip]{{images_hidden_degrees/trefoil/trefoil_normal_res.0125}.png}
        % \includegraphics[width=\linewidth,trim={6cm 2cm 6cm 2cm},clip]{{images_hidden_degrees/trefoil/trefoil_normal_res.0125}.jpg}
    \end{subfigure}
    \vspace{-10pt}
    \caption*{Resolution \(100\times100\times120\)}
    }

    \caption{
        %Intermediate states of the trefoil knot in different grid resolutions. The left and right columns show a side and a front view. Filaments are under-resolved in the half resolution setting and tend to shrink and lose details.
        Simulation of a trefoil knot with different grid resolutions, as viewed from the side (left column) and front (right column).
    }
    \label{fig:different_resolutions}
\end{figure}

Another important numerical parameter is the spatial resolution of the filament. Lagrangian methods constrain the curve resolution by subdividing and collapsing edges when they become too long or short. In contrast, our method controls the curve details via the resolution of the grid used for the level set $\psi$. \autoref{fig:different_resolutions} illustrates a simulation of an evolving trefoil knot on a \(50\times50\times60\) grid and one twice as detailed at \(100\times100\times120\). As expected, higher resolution grids create filaments with sharper details. Note, however, that the common way to visualize fluids is with marker particles or smoke densities, not by visualizing the filaments themselves. Thus, even very low resolution vortex filaments can still produce high resolution visual details. 
More detailed filaments make themselves evident via there more detailed velocity fields and complex smoke dynamics.

\begin{figure}[htbp!] 
 
     \begin{overpic}[width=0.32\linewidth,tics=0,trim={12cm -2cm 12cm 8cm},clip]
        {{images_hidden_degrees/trefoil/trefoil.0001}.png}
        % {{images_hidden_degrees/trefoil/trefoil.0001}.jpg}
    \end{overpic}
    \begin{overpic}[width=0.32\linewidth,tics=0,trim={12cm 4cm 12cm 2cm},clip]
        {{images_hidden_degrees/trefoil/trefoil_WP.0125}.png}
        % {{images_hidden_degrees/trefoil/trefoil_WP.0125}.jpg}
    \end{overpic}
    \begin{overpic}[width=0.32\linewidth,tics=0,trim={12cm 4cm 12cm 2cm},clip]
        {{images_hidden_degrees/trefoil/trefoil_WP.0150}.png}
        % {{images_hidden_degrees/trefoil/trefoil_WP.0150}.jpg}
    \end{overpic}
    \vspace{-20pt}
    \caption*{Wei\ss{}mann and Pinkall~\cite{Weissmann:2010:filaments}}

    \begin{overpic}[width=0.32\linewidth,tics=0,trim={12cm -2cm 12cm 8cm},clip]
        {{images_hidden_degrees/trefoil/trefoil.0001}.png}
        % {{images_hidden_degrees/trefoil/trefoil.0001}.jpg}
    \end{overpic}
    \begin{overpic}[width=0.32\linewidth,tics=0,trim={12cm 4cm 12cm 2cm},clip]
        {{images_hidden_degrees/trefoil/trefoil_ours.0125}.png}
        % {{images_hidden_degrees/trefoil/trefoil_ours.0125}.jpg}
    \end{overpic}
    \begin{overpic}[width=0.32\linewidth,tics=0,trim={12cm 4cm 12cm 2cm},clip]
        {{images_hidden_degrees/trefoil/trefoil_ours.0150}.png}
        % {{images_hidden_degrees/trefoil/trefoil_ours.0150}.jpg}
    \end{overpic}
    \vspace{-20pt}
    \caption*{Our method}

    \caption{Comparing a buoyant trefoil knot simulation by Wei\ss{}mann and Pinkall~\cite{Weissmann:2010:filaments} to ours. The simulations evolve from left to right. %The left and right columns show frames  at 125 and 150.
    } 
    \label{fig:trefoil_hidden}
\end{figure}

\begin{figure}[htbp!] 
\centering
    \begin{overpic}[width=0.4\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
            {{images_hidden_degrees/twoRings/2rings_0001}.png}
            % {{images_hidden_degrees/twoRings/2rings_0001}.jpg}
            % \put(30,0){ \small Initial state}
    \end{overpic}
    \vspace{-20pt}
    \caption*{Initial state}
    % \medskip
    
        % \medskip
    \begin{overpic}[width=0.4\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_ours.0070}.png}
        % {{images_hidden_degrees/twoRings/2rings_ours.0070}.jpg}
        % \put(20,0){ \small Ours Frame = 70}
    \end{overpic}
    % \hspace*{\fill}
    \begin{overpic}[width=0.4\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_ours.0100}.png}
        % {{images_hidden_degrees/twoRings/2rings_ours.0100}.jpg}
        % \put(20,0){ \small Ours Frame = 100}
    \end{overpic}
    % \hspace*{\fill}
    \vspace{-20pt}
    \caption*{Our method}

    \begin{overpic}[width=0.4\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_WP.0070}.png}
        % {{images_hidden_degrees/twoRings/2rings_WP.0070}.jpg}
        % \put(20,0){ \small W\&P Frame = 70}
    \end{overpic}
    % \hspace*{\fill}
    \begin{overpic}[width=0.4\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings_WP.0100}.png}
        % {{images_hidden_degrees/twoRings/2rings_WP.0100}.jpg}
        % \put(20,0){ \small W\&P Frame = 100}
    \end{overpic}
    % \hspace*{\fill}
    \vspace{-20pt}
    \caption*{Wei\ss{}mann and Pinkall~\cite{Weissmann:2010:filaments}, {\em with aggressive re-connection} }

 % \medskip
    \begin{overpic}[width=0.4\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings.WP_sameDOF.0078}.png}
        % {{images_hidden_degrees/twoRings/2rings.WP_sameDOF.0078}.jpg}
        % \put(20,0){ \small Ours Frame = 70}
    \end{overpic}
    % \hspace*{\fill}
    \begin{overpic}[width=0.4\linewidth,tics=0,trim={10cm 3cm 10cm 2cm},clip]
        {{images_hidden_degrees/twoRings/2rings.WP_sameDOF.0112}.png}
        % {{images_hidden_degrees/twoRings/2rings.WP_sameDOF.0112}.jpg}
        % \put(20,0){ \small Ours Frame = 100}
    \end{overpic}
    % \hspace*{\fill}
    \vspace{-20pt}
    \caption*{Wei\ss{}mann and Pinkall~\cite{Weissmann:2010:filaments}, {\em with conservative re-connection} }
    % \caption*{}
    
    % \vspace{-10pt}
    % \caption*{\small Frame = 70 \hspace{80pt} Frame = 100}

    %\caption{Comparison with Wei\ss{}mann and Pinkall~\cite{Weissmann:2010:filaments} by the linked rings. The left and right columns show frames  at 70 and 100. } \label{fig:two_rings}
    \caption{Simulating two linked vortex rings (top) with our method ($2^{nd}$ panel), and with the explicit filament approach of Wei\ss{}mann and Pinkall~\cite{Weissmann:2010:filaments} (bottom two panels). Lagrangian methods can be sensitive to numerical parameters for topological changes.} \label{fig:two_rings}
\end{figure}
    
%\subsection{Comparisons with \cite{Weissmann:2010:filaments}}
\subsection{Comparisons with Lagrangian filaments}
Next, we qualitatively compare implicit and explicit representations of curve dynamics using our algorithm and the Lagrangian vortex filament technique of Wei\ss{}mann and Pinkall~\cite{Weissmann:2010:filaments}, as implemented in Houdini software by SideFX.
 \autoref{fig:trefoil_hidden} shows the evolution of a knotted vortex filament with both methods. The filament is initialized as
 \begin{align}
    \gamma(s)=\left(\sin(s)+2\sin(2s),\cos(s)-2\cos(2s),-\sin(3s)\right) \notag
 \end{align}
%${\gamma(s)=\left(\sin(s)+2\sin(2s),\cos(t)-2\cos(2s),-\sin(3s)\right)}$
with \(s\in [0,2\pi)\).
 Aside from some small differences arising from the particulars of how filaments break apart and reconnect, the two methods produce roughly the same dynamics. 
 
On the other hand, the two methods have significantly different mechanisms for handling topological changes, which can produce divergent results. The Lagrangian method depends on user parameters like the thresholds for distances and angles between curves; the only relevant user parameter for changing topology in our method is the grid spacing, which prescribes the resolution of the level set $\psi$. \autoref{fig:two_rings} shows a simulation of two linked rings: our approach both preserves long thin tendrils and filters out topological noise. The Lagrangian simulation is sensitive to re-connection parameters: setting these parameters too aggressively leads to smooth geometry but loses thin features, while setting parameters too conservatively preserves thin features but creates noisy, persistent, high-speed ``ringlets'' that dominate the fluid velocity field. 

We stress here that the purpose of this comparison is to show how explicit and implicit descriptions handle topological changes of curves result in qualitatively different ways; the accuracy of these methods are not easily comparable as they have different mechanisms, and the accuracy depends on the type of curve dynamics. 

\subsection{Obstacles}

Like many methods for vortex dynamics, our method can also make filaments circumvent obstacles. 
A typical approach is to find a smooth harmonic potential \(\phi\) such that for a given obstacle \(B\subset M\), it solves 
\begin{align}\label{eq:pressure_projection}
     \langle    v_\gamma  -   v_{\partial B}-\nabla \phi,  n\rangle_{\RR^3}=0 \text{ on } \partial B \\
     \phi\rightarrow 0 \text{ at infinity} 
\end{align}
where \(   v_{\partial B}\) is the boundary velocity at each point of \(\partial B\).
Wei\ss{}mann and Pinkall \cite{Weissmann:2012:UnderWaterRigidBody} construct such a potential by regarding points inside \(B\) as sources of localized potentials so their weighted sum solves \autoref{eq:pressure_projection}, and Brochu et al. \cite{brochu2012linear} and Zhang et al. \cite{zhang2014pppm} solve a similar system using boundary element techniques.
 Nabizadeh et al. \cite{nabizadeh:2021:kelvin} address an equivalent problem by solving linear PDEs on infinite domains using the Kelvin transform.

Another approach by Wei\ss{}mann and Pinkall \cite{Weissmann:2010:filaments} is more specialized to vortex filaments and phenomena like vortex shedding; it regards the obstacle as a collection of artificial filaments \(\gamma_M\) such that the normal component of the velocity is zero.
Other approaches, like that of Park and Kim \cite{park2005vortex} and Da et al. \cite{da2015double} add point constraints to the boundary which zero out both the normal and tangential velocity components.
 
Our vortex filament algorithm is compatible with any of these obstacle-handling methods; our particular implementation uses \cite{nabizadeh:2021:kelvin}, as seen in \autoref{fig:obstacle}. We can observe that filaments near the obstacle are accelerated due to the induced mirrored image of themselves in the obstacle (or vortex sheet on the surface). 


Note that for the evolution of curves, only the evaluation of the velocity field needs the treatment of the boundary. Other components of the algorithm including the construction of levelset function ignore the boundary and do not require any additional treatments. 

\begin{figure}[t]
    \begin{minipage}[b]{1.0\linewidth}
        \centering
        \includegraphics[width=1.0\linewidth]{{images_hidden_degrees/obstacle/jet_obstacle_curves.0248.B}.jpg}
    %   \caption{obstacle curves}
    \end{minipage}
    \begin{minipage}[b]{1.0\linewidth}
        \centering
        \includegraphics[width=1.0\linewidth]{{images_hidden_degrees/obstacle/jet_obstacle.0248.B}.jpg}
    %   \caption{obstacle}
    \end{minipage}
    \caption{Turbulence caused by a spherical obstacle. We generate a new vortex filament ring every 10 frames.}
    \label{fig:obstacle}
\end{figure}


\begin{figure}
    \centering{
    \begin{subfigure}[b]{0.24\linewidth}
        \centering
        \includegraphics[width=\linewidth,trim={16cm 8cm 14cm 8cm},clip]{{images_hidden_degrees/twoRings/2rings_curvature_flow_big.0090}.png}
        % \includegraphics[width=\linewidth,trim={16cm 8cm 14cm 8cm},clip]{{images_hidden_degrees/twoRings/2rings_curvature_flow_big.0090}.jpg}
        % \caption*{Frame = 90}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.24\linewidth}
        \centering
        \includegraphics[width=\linewidth,trim={16cm 8cm 14cm 8cm},clip]{{images_hidden_degrees/twoRings/2rings_curvature_flow_big.0100}.png}
        % \includegraphics[width=\linewidth,trim={16cm 8cm 14cm 8cm},clip]{{images_hidden_degrees/twoRings/2rings_curvature_flow_big.0100}.jpg}
        % \caption*{Frame = 100}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.24\linewidth}
        \centering
        \includegraphics[width=\linewidth,trim={16cm 8cm 14cm 8cm},clip]{{images_hidden_degrees/twoRings/2rings_curvature_flow_big.0130}.png}
        % \includegraphics[width=\linewidth,trim={16cm 8cm 14cm 8cm},clip]{{images_hidden_degrees/twoRings/2rings_curvature_flow_big.0130}.jpg}
        % \caption*{Frame = 130}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.24\linewidth}
        \centering
        \includegraphics[width=\linewidth,trim={16cm 8cm 14cm 8cm},clip]{{images_hidden_degrees/twoRings/2rings_curvature_flow_big.0155}.png}
        % \includegraphics[width=\linewidth,trim={16cm 8cm 14cm 8cm},clip]{{images_hidden_degrees/twoRings/2rings_curvature_flow_big.0155}.jpg}
        % \caption*{Frame = 130}
    \end{subfigure}
   
    }

    \caption{
        Curves evolving under the curve shortening flow, from left to right.
    }
    \label{fig:curve_shortening_flow}
\end{figure}

\subsection{Other types of filament dynamics}
Our idea of implicit representation of filament dynamics is not limited to vortex filaments. In theory, this should be applicable to any first-order time evolution of curves. \autoref{fig:curve_shortening_flow} shows an example with the curve-shortening flow,
\begin{align}
    \frac{\partial\gamma}{\partial t}(s,t)=-\gamma^{\prime\prime}(s,t).
\end{align}
As expected, this velocity causes the filament to shrink over time, similar to mean-curvature flow for surfaces~\cite{osher1988fronts}.


\arrayrulecolor{black}
\begin{table*}[htbp!]
    \caption{Parameters and timing breakdown per frame for all simulations in this chapter and the accompanying video. Symbols \(\Gamma\) and \(a\) are the intensity and the thickness of filaments. All simulations are 24fps and and the time step size \(1.0/24\)s except for \autoref{fig:curve_shortening_flow}, which we used 120fps. Average timings are taken over the entire simulation. The ``Other'' column includes the remaining operations, including obstacle handling and generation of new curves at the sources in \autoref{fig:jet} and \autoref{fig:obstacle}. }

\centering
    %Shrink table size
\scalebox{0.75}{
    \begin{tabular}{|c|| c | c |c ||c|| c|c|c|c|c|c|} 
     \hline
     Scene & Resolution & \(\Gamma\) & \(a\)  & Total time & Evaluate \(   v\) &  Advect \(\psi\) & Construct \(\gamma\) & Construct \(\psi\) & Other\\ 
    %  [0.5ex] 
     \hline
     \hline
     \autoref{fig:two_rings}  & \(100 \times 100 \times 100\) & 1.0 &0.05 & 0.159s & 0.041s  & 0.017s &0.041s& 0.056s & 0.002s \\ 
     \hline
     \autoref{fig:trefoil_hidden}   & \(100 \times 100 \times 120\) & 1.0 &0.05 & 0.274s  & 0.056s  & 0.026s &0.049s& 0.102s & 0.001s\\ 
     \hline
     \autoref{fig:leapfrog}   & \(60 \times 60 \times 200\) & 2.0 &0.08 & 0.090s  & 0.013s &0.021s & 0.023s & 0.031s & 0.001s\\ 
     \hline
     \autoref{fig:oblique_rings}   & \(80 \times 120 \times 140\) & 1.0 &0.05&  0.152s  & 0.031s  & 0.030s & 0.046s & 0.044s & 0.001s\\ 
     \hline
     \autoref{fig:torus}   & \(100 \times 100 \times 140\) & 1.0 &0.05 & 0.655s & 0.144s  & 0.085s & 0.091s & 0.332s & 0.002s \\ 
     \hline
     \autoref{fig:jet}   & \(100 \times 100 \times 180\) & 2.0 &0.08 & 1.423s  & 0.177s  & 0.099s & 0.812s & 0.324s & 0.011s\\ 
     \hline
     \autoref{fig:obstacle}   & \(100 \times 100 \times 160\) & 2.0 &0.08 & 0.969s  & 0.047s  & 0.087s & 0.473s & 0.226s & 0.135s\\ 
     \hline
     \autoref{fig:curve_shortening_flow}  & \(100 \times 100 \times 100\) & N/A &N/A & 0.152s  & 0.028s  & 0.012s & 0.068s & 0.043s & 0.001s\\ 
     \hline
    \end{tabular}
}% end scalebox
    \label{table:parameters}
\end{table*}

\begin{table}[hbtp!]
    \caption{Computational timings compared with Houdini's  built-in implementation of Weissmann \& Pinkall \cite{Weissmann:2010:filaments}. ``Same DOF'' refers to simulations with approximately the same number of computational degrees of freedom as our method: we set the relevant parameters (re-connection distance, minimum and maximum edge lengths ) so that the number of explicit curve vertices are similar to ours.
    % For the same DOF comparisons, we set #Vertices = 320 for the linked rings and #Vertices = 640 for the trefoil.
    }
%Shrink table size
\centering
\scalebox{0.75}{
    \begin{tabular}{|c|| c |c|c|} 
     \hline
     Scene description &  W\&P &  W\&P (Same DOF) &  Ours  \\ 
    %  [0.5ex] 
     \hline
     \hline
     Linked rings (\autoref{fig:two_rings}) & 0.052s &  0.065s& 0.159s  \\ 
     \hline
     Trefoil knot (\autoref{fig:trefoil_hidden} & 0.045s & 0.101s&0.274s  \\ 
     \hline
    \end{tabular}
}% End scalebox
    
    \label{table:comparison_with_WP}
\end{table}

\subsection{Discussion}

Our method represents a fundamental new way to animate fluids and vortex filament dynamics. Its unique features give it some strengths and weaknesses relative to existing approaches. 

First of all, we represent our filaments with implicit functions $\psi$. While there are many ways to encode an implicit function, our implementation uses a regular grid, which implies a finite bounding box. While common for Eulerian fluid simulations, bounded domains are a constraint not shared by Lagrangian methods. This constraint could be mitigated if we use sparse grids or trees~\cite{museth2019openvdb}.

The main parameter in our method is the grid resolution. As discussed earlier, this parameter affects the geometric detail and topology of our filaments, influencing the velocity field directly, but influencing the final visual results only indirectly. 
The grid resolution also directly influences topological changes; the only mechanism for topological changes in our algorithm is to merge curves when they intersect the same grid cell. We believe this automatic and robust method for handling topological changes is a strength of our method; it minimizes the need to fine-tune parameters (especially the relationship between minimum/maximum edge length and topological change interaction lengths) and seems immune to the types of numerical blow-ups that we have seen in Lagrangian methods when filaments get close or exhibit near-degenerate geometry. 

On the other hand, this grid-based method for re-connecting curves will also delete small features when they shrink below the grid cell size (similar to level set methods for surfaces). This behavior is most evident in our jet example (\autoref{fig:jet}). We believe that the deletion of small features can be reduced in the future in a number of ways. More accurate advection schemes will probably preserve higher frequency features of $\psi$ better without deleting them. Also, although we took care to introduce a $   v$ and $\psi$ which alleviates egregious numerical stability problems, we have not quite optimized for accuracy or geometric durability for the specific dynamics of vortex filaments. Finding a pair of $   v$ and $\psi$ that is more suitable for each curve dynamics in the degrees of freedom discussed in \autoref{sec:representation} would lead to even better numerical performance.

Although our method relies on Eulerian advection to evolve the geometry, it {\em does  not} suffer from the artificial viscosity typically associated with Eulerian fluid simulations. Our fluid velocity is reconstructed from filament dynamics, so the velocity field is not recursively re-sampled and does not accumulate damping errors. 
Consequently, our method produces swirly and energetic fluid flows even at low grid resolutions. 

\autoref{table:parameters} lists the simulation parameters and timing breakdown for each of the simulations in this chapter and accompanying video. We stress that our prototype implementation is meant as a proof of concept, and it has plenty of room for optimization. Our current implementation employs regular grids for ease of implementation; future implementations can make great use of sparse grids, since our filaments only use a 1-dimensional path through the 3D grid. Our implementation of ``Evaluate \(   v\)'' and ``Construct \(\psi\)'' iterates over the entire curve geometry for each point in the narrow band where \(\gamma\) and \(\psi\) are required; future implementations could use a sparse data structure and approximated fast summation via tree-codes like the fast multipole method to reduce total evaluation time. Additionally, our current implementation of ``Construct \(\gamma\)'' redundantly doubles the work per curve vertex, so that step can be sped up by at least a factor of 2. 

With these inefficiencies in mind, our implementation appears to run modestly slower than Houdini's optimized implementation of  \cite{Weissmann:2010:filaments} for the scenarios we tested. The two methods have completely different numerical degrees of freedom, so we find it difficult to compare them directly.
\autoref{table:comparison_with_WP} compares the computational cost of the two methods for figures in this chapter. 
We aimed to keep the number of curve vertices roughly the same as our method's in the ``Same DOF'' simulations, so we believe these are the most relevant for comparing timings. 

In our examples of vortex filaments,
we assumed inviscid filaments, but we can handle viscous motions too once the time-derivative of viscous curves is speficied.
We also assumed a uniform vortex strength. Just as in the explicit description, our implicit description requires a well-defined equation of motion for curves, and it outputs only new curve configurations without processing additional quantities like vortex strengths. Simulating vortex filaments with different strengths is an independent challenge. 
Merging filaments with different strengths would create a graph of filaments rather than just disjoint closed curves, which would require a new implicit representation for filament graphs. 


\section{Conclusion and future work}
\label{sec:conclusion}
We have shown that implicit representations of geometric curves exhibit large degrees of freedom in both their mathematical representation, as well as their dynamics. We then took advantage of these redundant degrees of freedom to improve the stability of vortex filament simulations. 

We see a number of avenues that can be explored in future work. 
Our current strategy exploits redundancy in the 3D velocity field and level set outside the filaments, but we can also incorporate free degrees of freedom in the tangential components of the velocity on the filament itself.
More generally, one can attempt to formalize the regularization of $   v$ and $\psi$ for various purposes. A possible instance is an optimization problem for certain energies aiming provable guarantees on numerical accuracy and stability. Besides, the use of redundancy in dynamics should also be possible for other codimensional cases such as level surfaces in 3D or level sets of an arbitrary codimension in a higher dimensional space.

Our numerical scheme can be made more sophisticated as well: higher order advection schemes and geometric curve representations, as well as sparse and adaptive grids can make our method both more efficient and more numerically accurate. 

Finally, this chapter explores vortex filaments and curve shortening flows, but our ideas are not limited to these specific dynamical systems. We expect that the idea of exploiting hidden degrees of freedom in implicitly represented curve dynamics will generalize to many more types of dynamics appearing in both scientific fields and engineering applications. 


