\documentclass[lettersize,journal]{IEEEtran}
\usepackage[utf8]{inputenc}
\DeclareUnicodeCharacter{3002}{\textvisiblespace}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{tabu}
\usepackage{booktabs}
\usepackage{lipsum}
\usepackage{mwe}
\usepackage{mathptmx} 
\usepackage{multirow}
\usepackage{multicol}
\usepackage{algorithm2e}
\usepackage{float}
\usepackage{xcolor}
\usepackage{enumitem}
\usepackage{cleveref}
\usepackage{thmtools,thm-restate}

\newtheorem{definition}{Definition}
\newtheorem{lem}{Lemma}
\newtheorem{assumption}{Assumption}
\newtheorem{col}{Corollary}
\newtheorem{pro}{Proclaim}

\SetKwComment{Comment}{/* }{ */}
\RestyleAlgo{ruled}
\newcommand{\tabincell}[2][c]{\begin{tabular}[#1]{@{}c@{}}#2\end{tabular}}

\newcommand{\HG}[1]{\textcolor{red}{HG: [[#1]]}}
\newcommand{\CR}[1]{\textcolor{blue}{CR: #1}}
\newcommand{\removed}[1]{\textcolor{red}{\sout{#1}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\title{Meshing Deforming Spacetime for\\Visualization and Analysis}
\author{Congrong~Ren and~Hanqi~Guo
\IEEEcompsocitemizethanks{
\IEEEcompsocthanksitem C. Ren and H. Guo are with the Department of Computer Science and Engineering at The Ohio State University, Columbus, OH 43210, USA.\protect\\
E-mails: \{ren.452, guo.2154\}@osu.edu%
\IEEEcompsocthanksitem H. Guo is also with the Mathematics and Computer Science Division, Argonne National Laboratory, Lemont, IL 60439, USA.\protect\\
E-mail: hguo@anl.gov}% <-this % stops an unwanted space
% \thanks{Manuscript received April 19, 2005; revised August 26, 2015.}
}

\markboth{IEEE~Transactions~on~Visualization~and~Computer~Graphics,~Vol.~X, No.~Y, August~2023}%
{Ren \MakeLowercase{\textit{et al.}}: Meshing Deforming Spacetime for Visualization and Analysis}

\IEEEtitleabstractindextext{
\begin{abstract}
We introduce a novel paradigm that simplifies the visualization and analysis of data that have a spatially/temporally varying frame of reference. The primary application driver is tokamak fusion plasma, where science variables (e.g., density and temperature) are interpolated in a complex magnetic field-line-following coordinate system.  We also see a similar challenge in rotational fluid mechanics, cosmology, and Lagrangian ocean analysis; such physics implies a deforming spacetime and induces high complexity in volume rendering, isosurfacing, and feature tracking, among various visualization tasks.  
% 
Without loss of generality, this paper proposes an algorithm to build a simplicial complex---a tetrahedral mesh, for the deforming 3D spacetime derived from two 2D triangular meshes representing consecutive timesteps. Without introducing new nodes, the resulting mesh fills the gap between 2D meshes with tetrahedral cells while satisfying given constraints on how nodes connect between the two input meshes. In the algorithm we first divide the spacetime into smaller partitions to reduce complexity based on the input geometries and constraints. We then independently search for a feasible tessellation of each partition taking nonconvexity into consideration.  We demonstrate multiple use cases for a simplified visualization analysis scheme with a synthetic case and fusion plasma applications.
% and fluid dynamics applications.
% 
% We propose a framework for constructing simplicial mesh for deforming 3D spacetime. Given a 2D triangular mesh and its varying dependency between mesh nodes in successive timesteps, our framework build a tetrahedral mesh in 3D spacetime domain. Our tetrahedralization method is based on \textit{divide-and-conquer} paradigm, which splits the 3D deforming spacetime into smaller pieces first and eliminates nodes one by one by a decision-tree search algorithm. The derived simplicial spacetime meshing provides us a discretization of continuous time-varying scientific data as well as a common ground to fill the gaps between discrete timesteps. We present how our framework can be used on time-varying meshes including XGC mesh, adaptive multi-resolution (AMR) meshes, and mesh built on flow field around turbine cascade, and show the feature tracking results on them.
\end{abstract}

\begin{IEEEkeywords}
Spacetime meshing, triangulation, time-varying mesh, isosurfacing, volume rendering.
\end{IEEEkeywords}}

\maketitle
\IEEEdisplaynontitleabstractindextext
\IEEEpeerreviewmaketitle

% \IEEEraisesectionheading{\section{Introduction}
\section{Introduction}
\label{sec:intro}
\IEEEPARstart{S}{cience} applications often involve spacetime that is constantly changing.  For example, the cosmos evolved from a dense matter and expanded to today's universe with sparsely distributed galaxies~\cite{takle2012case, woodring2011situ}.  In the general theory of relativity, distortions of spacetime are created by mass and energy~\cite{dirac1996general}.  In fluid dynamics, fluid parcels may be rotated, translated, compressed, and expanded, driven by flow transport~\cite{dutta2016situ, koehler2011vortex}.  Likewise, in electrodynamics, electromagnetic fields deform infinitesimal volumes of particle parcels~\cite{altenbach2016inelastic}.

% Deforming spacetime are derived from scientific scenarios involving time-varying physical behaviors, such as fluid dynamics of flows induced by rotating machinery~\cite{dutta2016situ, koehler2011vortex}, evolution of matters in cosmological simulations~\cite{takle2012case, woodring2011situ}, and eddy analysis in ocean models~\cite{hollt2014ovis, woodring2015situ}. 

In our observation, to date, \emph{deforming spacetime} is either overlooked or induces high complexity in scientific visualization.  In the former case, physical insights between discrete timesteps are ignored because usually the time-varying data are visualized only at discrete timesteps.  In the latter case, deforming spacetime challenges visualization and analysis tasks. For example, nonlinear physics-driven interpolation incurs prohibitively high computational costs in visualization tasks such as volume rendering: interpolating science variables (e.g., density and temperature) in tokamak fusion plasma\footnotemark~requires tracing streamlines in a magnetic field~\cite{white2017theory}. Furthermore, such a nonlinear interpolation basis makes it even more difficult for root-finding tasks, such as extracting isosurfaces where the interpolated value equals a given value and finding critical points where interpolated vector field values vanish.  Later sections will refer to the interpolation and root-finding processes as \emph{forward} and \emph{inverse approximations}, respectively, as we further motivate visualization with deforming spacetime.
\footnotetext{Technically, this paper is not concerned with physical time in tokamaks but presents 3D deformation by treating toroidal direction ($\phi$) as time without losing generality, as explained in case studies.}

% \HG{I would avoid wordings like huge, enormous, and impressive; they are too marketing and oftentimes non-scientific.  Use plain, specific, and accurate languages will be appropriate for research papers.}

% However, deforming property of spacetime introduces huge complexity to visualization and analysis tasks including upsampling, volume rendering, and feature extraction and tracking due to the difficulty on discretization: high temporal resolution rises both storage and computational cost, while low temporal resolution weakens the accuracy of interpolation.

% Scientists meet with analyses of data behavior over space and time in many scenarios, such as evolution of matters in cosmological simulations~\cite{takle2012case, woodring2011situ}, eddy analysis in ocean models~\cite{hollt2014ovis, woodring2015situ}, and fluid dynamics of flows induced by rotating machinery~\cite{dutta2016situ, koehler2011vortex}. One computational way is formulating a discretization for the space (i.e., meshing the space) while treating time dimension separately. Another possibility is directly meshing the \textit{spacetime}, a coupled domain of space and time dimensions which is consist of time ``slabs'' generated by extruding either fixed or time-varying spatial mesh along temporal direction while building correspondence between nodes in two meshes of consecutive timesteps~\cite{diosady2018linear}. Discretization of spacetime not only corresponds with fundamental physical laws of nature (e.g., conservation law)~\cite{lowrie1998space} but also drastically reduces the complexity compared to time-marching methods when tracking position of propagating features~\cite{ishii2019solving, yano2012optimization}.

In this work we explore a novel paradigm to represent deforming spacetime by (simplicial) meshes to simplify the visualization and analysis of data with spatially/temporally varying frames of reference.  Our spacetime mesh offers a continuous (and also invertible for root-finding problems) data representation of time-varying data, making forward/inverse approximation possible with less computational cost.  For example, volume rendering in the fusion plasma applications relies on a complex and nonlinear magnetic-driven interpolation; with the deformed space represented as an unstructured grid, the data can be directly visualized with off-the-shelf visualization tools.  Likewise, isosurfacing in the deformed space is also made possible with the spacetime meshing scheme.
% We introduce a novel paradigm that simplifies the visualization and analysis of data that have spatially/temporally varying frame-of-references. 
Although our primary application driver is tokamak fusion plasma, %, where science variables (e.g., density and temperature) interpolate with a complex magnetic-line-following coordinate system.  
 our methodology is applicable also in rotational fluid mechanics~\cite{dutta2016situ, koehler2011vortex}, cosmology~\cite{takle2012case, woodring2011situ}, and Lagrangian ocean analysis~\cite{friederici2018lagrangian, van2018lagrangian}; such physics implies a deforming spacetime and induces high complexity in volume rendering
% , \poincare plots\HG{I would entirely remove Poincare plots in this paper.},
and feature extraction and tracking, among various visualization tasks.

% \HG{You did not mention what are the inputs and constraints of our algorithm but directly started talking about the challenges.  I hereby added the following two paragraphs.}

This paper focuses on meshing 3D spacetime (2D space and 1D time) and assumes that the inputs are a series of 2D simplicial meshes, each representing the spatial discretization in a timestep.  The assumption reflects real application scenarios: most simulations output data as a spatial grid. The deformation of spacetime is explicitly defined as temporal connectivity of mesh nodes in adjacent steps (i.e., the ``next node'' in tokamak applications, as explained in Section~\ref{sec:spacetime_formulation}).  In other words, the connectivities are a discrete representation of how spatial locations are transformed in adjacent steps. The temporal connectivities may be available from data % (e.g., fusion plasma simulations, as further explained in Section~\ref{sec:spacetime_formulation}) 
or derived by tracing particle trajectories based on the underlying physics.  % domain knowledge.  
% \footnotetext{Depending on applications, the 2D mesh structure may be identical over timesteps; still, one can define deformation with temporal connectivity.}


Spacetime meshing from spatial meshes and temporal connectivities imposes multiple geometrical constraints, leading to a convoluted problem that requires a computationally feasible method.  The objectives and constraints for the resulting spacetime mesh include the following:
\begin{enumerate}[topsep=1ex,itemsep=-1ex,partopsep=1ex,parsep=1ex]
  \item The spacetime mesh includes all simplices (e.g., triangles, edges, and nodes) from the input spatial meshes without adding new nodes (i.e., \textit{Steiner points}).
  \item The spacetime mesh contains all edges that connect nodes in adjacent steps, as specified by temporal connectivities.
  \item The spacetime mesh is simplicial (containing only tetrahedral cells) for visualization and analysis tasks.
\end{enumerate}

\begin{figure*}
    \centering
    \includegraphics[width=0.9\textwidth]{figs/spacetime.pdf}
    \caption{
    % \HG{Avoid using subfloat but just make the figure as one single file; refer to subfigures with (a), (b), and (c).  This is just a personal preference but can make it easier to adjust sizes.} 
    Two examples of spacetime meshes (a and c) and their simplicial subdivisions (b and d). 
    % Two examples of spacetimes and their simplicial subdivisions. 
    The prismatic mesh in (a) is obtained by extruding a 2D spatial mesh. One possible triangulation is given by staircase triangulation~\cite{de2010triangulations} in (b).
    % \HG{cite the staircase paper here.}
    Time-varying behavior introduces deforming spacetime in (c), which has more complex correspondence between the lower and upper meshes: nodes $n_3$ and $n_4$ have the same next node; node $n_x$ has no previous nodes; triangulations in the lower and upper meshes are different; the lower mesh needs translating, rotating, or hybrid transformation to be aligned with the upper mesh. All these differences make it difficult building a simplicial meshing (d) for deforming spacetime.}
    \label{fig:spacetime}
\end{figure*}


% % \HG{I would avoid using the term topology because it may be interpreted as other things.}
% Meshing the deforming spacetime can be prohibitively complex, and we observe several technical challenges. 
% % This work provides a framework for constructing simplicial mesh for deforming 3D spacetime extruded by 2D simplicial meshes, with consideration of three main challenges.
% % \HG{Do we have to assume 2D meshes being fixed?}
% % The challenges of subdividing deforming 3D spacetime into 3-simplices are threefold. 
% First, temporal dependencies between nodes often make the spacetime nonconvex (both locally and globally), making it challenging to leverage existing algorithms for convex polytopes.  
% % Some connections of mesh nodes might cross the surface or go outside of the volume covered by spacetime and are not allowed to create.
% Second, even without considering constraints, the time complexity of triangulating a 3D nonconvex polyhedron is given by $O(n_v^3)$ ($n_v$ stands for the number of mesh nodes in a spacetime)~\cite{si2019simple}, which is considerably costly for meshes in real applications; incorporating constraints from next nodes will further complicate the triangulation.
% our methodology is the first attempt to search a triangulation spacetime with the constraints defined by spatial meshes and next nodes.
% % Third, given connectivity of the spatial mesh and temporal connections between mesh nodes given by physical behaviors introduce constraints on triangulation of spacetime. Although spacetime forms a polytope, the surface of it is not composed of polygons but triangular mesh instead, and the temporal connections inside of it highly narrow down the scope of feasible tetrahedralizations as well.

Considering the constraints induced by the spatial meshes and the temporal connectivities, meshing the deforming spacetime can be extremely complex;
% To the best of our knowledge, 
a new methodology is required to search for a triangulation of the spacetime subject to the above-described constraints with a reasonable time cost.  First, most existing triangulation algorithms are for convex polyhedra; however, the temporal dependencies between nodes often make the spacetime nonconvex (both locally and globally).  Second, even without considering constraints, the time complexity of triangulating a 3D nonconvex polyhedron is $O(|\mathcal{V}|^3)$~\cite{si2019simple} ($|\mathcal{V}|$ stands for the number of mesh nodes in a spacetime), which is considerably costly for meshes in real applications.

% To cope with all the challenges above, 
In this paper we propose a divide-and-conquer algorithm that divides the spacetime into small and independent partitions (the \textit{divide} stage) and then triangulates each partition (the \textit{conquer} stage). 
%GAIL - "searches triangulation" sounds odd -- do you mean "triangulates" or searches for possible triangulations in?
The divide stage aims to reduce computational costs by decomposing the problem into independent and manageable subproblems of triangulating partitions.  For shared boundaries of neighboring partitions, we define a rule to triangulate the boundaries so that each partition can be processed independently.  
% speculate how the boundaries are triangulated
% Besides cutting the spacetime into as small as possible blocks, if tetrahedralization of each block is independent from those of other blocks and if the triangulation on the cross sections of neighboring partitions is identical are also taken care of. 
The conquer stage focuses on each individual partition, which is triangulated by iterating over nodes and eliminating all tetrahedra that connect to a node. 
% we eliminate nodes one by one by removing tetrahedra incident to them. 
A decision-tree-based search is applied to decide which node to eliminate and how to triangulate its surrounding with all geometry (e.g., nonconvexity) and connectivity constraints considered.  Once we triangulate all partitions, the resulting spacetime mesh can enable and simplify many visualization tasks, and we demonstrate use cases of our spacetime meshes with tokamak fusion plasma data and a synthetic dataset. 
% This algorithm is tested on tokamak fusion plasma and a synthetic dataset.
In summary, our framework makes the following contributions:
\begin{itemize}
    \item A novel paradigm to simplify visualization of scientific data that have spatially/temporally varying frames of reference
    \item A method to partition 3D deforming spacetime into smaller independent components for triangulation
    % time complexity reduction method by segmenting a 3D deforming spacetime into smaller independent partitions.
    \item A decision-tree-based algorithm for searching for feasible triangulation for a 3D deforming spacetime with geometry and connectivity constraints.
    %GAIL - searches for feasible? 
    % \item A generalization scheme for spatiotemporal
    % \HG{spatiotemporal, no hyphen}
    % data visualization and analysis.
    % of algorithm. \HG{I don't understand what you mean by ``data analysis of algorithm.''}
\end{itemize}

% \HG{Put the boilerplate here: ``The remainder of this paper is organized as follows...''}
The remainder of this paper is organized as follows.  Section~\ref{sec:related} reviews related work, and Section~\ref{sec:preliminaries} presents fundamentals and the driver application of our methodology.  We describe the divide-and-conquer algorithm in Section~\ref{sec:method}.  Section~\ref{sec:case-studies} demonstrates the uses cases and effectiveness of our method both quantitatively and qualitatively. Section~\ref{sec:discussion} discusses limitations and future work. We conclude with a brief summary in Section~\ref{sec:conclusion}. 

% Cases studies in Section 6 verifies the effectiveness of our method quantitatively and qualitatively. At last, in section 7, we conclude the paper with a summary, limitations, and future work.
% In Section 2, how feature tracking tasks motivate this work, what has been achieved by existing feature tracking algorithms, and how our tetrahedralization methodology is different than triangulation in computational geometry are discussed. In Section 3, we describe an example of deformed spacetime formed by magnetic field in fusion plasma, which drives this work. Immediately in Section 4, we formulate the problem of meshing 3D spacetime into a geometric problem. Section 5 describes our methodology in detail. 


\section{Related Work}
\label{sec:related}

% This work is closely related to research topics including feature tracking, computational science, and triangulation.
We review literature on spacetime meshing and triangulation. 
% for visualization and computational science.  We also review possible triangulation schemes for spacetime meshing. 

 % feature tracking, spacetime meshing, and triangulation. 

% We review literature related to feature tracking, spacetime mesh, and triangulation in geometry to illustrate motivation and challenges of this work. The work is motivated by feature tracking tasks, which are to identify trajectories of features (i.e., objects of interest including critical points, vortices, vortex cores, and isosurfaces) over a series of timesteps. There have been some well-known feature tracking algorithms,; different than these works, we focus on challenges caused by deforming spacetime. Besides feature tracking, we also review other tasks that can be solved when spacetime mesh is given. Our formulated problem is similar to 3D triangulation problem in computational geometry domain, but has more mesh connectivity and geometry constraints.

% We briefly review literature related to feature tracking, spacetime mesh, and triangulation in geometry to illustrate motivation and challenges of this work. The work is motivated by feature tracking tasks, which are to identify trajectories of features (i.e., objects of interest including critical points, vortices, vortex cores, and isosurfaces) over a series of timesteps. There have been some well-known feature tracking algorithms,; different than these works, we focus on challenges caused by deforming spacetime. Besides feature tracking, we also review other tasks that can be solved when spacetime mesh is given. Our formulated problem is similar to 3D triangulation problem in computational geometry domain, but has more mesh connectivity and geometry constraints.

% \HG{Technically, this is Sec 2.0 here, and you have to write a couple of lead sentences here.}


\subsection{Spacetime Meshing for Scientific Visualization}
\label{sec:feature_tracking}
% \HG{Missing the rationale on how/why our research is related to feature tracking.  It may be worth reviewing more on the spacetime approach.  Possible to merge with 2.1?}

The visualization research community has been investigating spacetime meshing since the early 2000s, primarily motivated by extracting and tracking features that evolve over the spacetime continuum.  Spacetime meshing provides a mathematically sound basis for feature tracking; however, to our knowledge, no methodology has considered meshing deforming spacetime for feature tracking.  As exemplified below, the fundamental rationale of meshing spacetime for feature tracking is establishing a continuous data representation to generalize spatial feature descriptors to spacetime directly.  

% feature tracking problems, and thus enables us to describe how features evolve over time. Besides spacetime approaches, there are other two classes of feature tracking algorithms according to how they treat time dimensions, which also leads to different methods for feature extraction and matching problem (i.e., to build correspondence between features in consecutive timesteps that depict the same feature at different times)\cite{post2003state}.
% It is involved with both feature extraction and matching problem, i.e., to build correspondence between features in consecutive timesteps that depict the same feature at different times. Different tracking approaches extract and match features by different methods. There are three categories of them according to how they treat time dimension, which is also related to how they solve matching problem~\cite{post2003state}.

\textbf{Critical point tracking in vector fields.}  In the early work by Tricoche et al.~\cite{tricoche2002topology} for tracking critical points in time-varying 2D vector fields, the authors constructed a 3D mesh by first placing and copying the spatial triangulation for different timesteps and then building a prismatic cell for each corresponding triangle pair in adjacent timesteps, as illustrated in Figure~\ref{fig:spacetime}(a).  Garth et al.~\cite{garth2004tracking} further enabled 3D critical point tracking based on 4D prisms.  Assuming a piecewise linear interpolation applies to the spatial triangulation and a linear interpolation over time, one can derive the exact location and time of critical points in the prism cells' boundaries.  One can further define and reconstruct critical point trajectories and identify events such as birth/death and split/merge. 

\textbf{Singularity tracking in complex-valued scalar fields.} A similar success of generalizing feature descriptors to 4D is the tracking of vortices---singularity curves in 3D complex fields (represented as magnitudes and phases)---from superfluidity, superconductivity, and Bose--Einstein condensates.  In 3D, a vortex is a locus of points where the contour integral of the phase field over an infinitesimal loop is nonzero~\cite{phillips2015detecting}.  The 4D definition of vortices~\cite{guo2017situ, guo2015extracting, phillips2016tracking} generalizes the contour integral to spacetime so that one can reconstruct the trajectory surfaces of a vortex with the 4D mesh complex.  

\textbf{Isosurface tracking in scalar fields.} We view higher-dimensional isosurfacing~\cite{BhaniramkaWC04}---a generalization of marching cubes~\cite{LorensenC87} in 4D and beyond---as a tracking technique based on spacetime meshes. Because there are a finite number of possible ways for a hypercube to intersect an isosurface, one can establish a lookup table to extract time-varying isosurfaces in 4D. Inspired by marching tetrahedra~\cite{DoiK91} that reduces ambiguity cases in marching cubes, Guo et al.~\cite{guo2021ftk} developed a higher-dimensional simplicial meshing to help eliminate ambiguous cases in isosurface tracking.

% \textit{Spacetime approaches} regard time as an additional dimension to spatial dimensions. Features are extracted directly in spacetime domain, and tracked by temporal interpolation. Tracking by linear (including barycentric, bilinear, and trilinear) interpolation assumes the field is linear in every cell and thus piecewise-linear (PL) over spacetime domain~\cite{garth2004tracking, guo2021ftk, tricoche2002topology}, which puts forward higher requirements on reasonable spacetime meshing: the values of cell nodes need to be highly related to the values of other points in the cell.

% Spacetime meshing has been used to track features like singularities in vector~\cite{garth2004tracking, guo2021ftk, tricoche2002topology} or phase~\cite{guo2015extracting, guo2017situ, phillips2015detecting, phillips2016tracking} fields in scientific visualization domain. A spacetime is generated by connection of spatial meshes for successive timesteps. However, most of the works only consider prismatic meshes or subdivision of prismatic meshes, and thus are not applicable for deforming spacetime where ``prisms'' are squeezed or twisted into other shapes, which motivates this work.

We refer readers to the literature~\cite{post2003state} for a comprehensive review of feature tracking in scientific visualization.  
% 
For example, \textit{isolated time approaches} first extract features in every timestep and then determine their correspondence using distance metrics as matching criteria.  Features can be extracted by numerical~\cite{liang2020toward} or topological~\cite{carr2003computing} methods.  Some link two features if their distance is lower than some threshold~\cite{laney2006understanding}.  Some apply topological descriptors based on Morse theory and search for the pairing of nodes (i.e., critical points) in two graphs (e.g., merge tree~\cite{yan2023multilevel}, Reeb graph/contour tree~\cite{tu2019propagate}) that minimizes the distance between the two graphs.
% 
Some approaches \textit{implicitly incorporate temporal dimension}, such as the feature flow field (FFF) method~\cite{theisel2003feature} and its variants~\cite{reininghaus2011efficient, weinkauf2010stable}. Theisel and Seidel proved that critical lines in a $n$-D scalar field are equivalent to streamlines of a derived $(n+1)$-D vector field.
%so-called FFF~\cite{theisel2003feature}.  
This equivalence transfers a critical point tracking problem to a streamline integration problem.  Stable FFF is then introduced to reduce accumulated errors incurred by numerical integration~\cite{weinkauf2010stable}.  Combinatorial FFF further filters out less important critical lines by their integrated persistence to make FFF more robust for noisy data~\cite{reininghaus2011efficient}.

Besides feature tracking, some other works are related to spacetime. For example, in flow visualization, spacetime is used to define the domain of pathlines in time-dependent vector fields. Wilde et al.~\cite{wilde2020flow} presented a technique to modify \textit{flow maps} in spacetime according to space deformation, which maps a particle seeded at some position and time to its destination after a given time interval and thus explicitly encodes pathlines. G{\"u}nther~\cite{gunther2020visibility} considered flow visualization in applications related to finite-sized particles, where the trajectories of particles are described in a spacetime domain.

% \subsection{Spacetime Triangulation for Spacetime Meshing} % in Geometry}
\subsection{(Spacetime) Triangulation} % in Geometry}

We review several basic concepts in triangulation and the most relevant literature on it. While few directly tackle triangulation in spacetime, triangulation (primarily in 2D and 3D spaces for real-world applications) is widely studied in geometric modeling, computer graphics, and computational sciences; readers are referred to the work of De Loera et al.~\cite{de2010triangulations} for a comprehensive understanding.

\textbf{Simplices and simplicial prisms}. An \textit{$n$-simplex} is a convex hull with $n+1$ nodes that are affinely independent in $\mathbb{R}^n$. For example, a 0-simplex is a point, a 1-simplex is a line segment, a 2-simplex is a triangle, and a 3-simplex is a tetrahedron. An \textit{$(n+1)$-D simplicial prism} (also referred to as $(n+1)$-prism or prism) is derived by extruding an $n$-simplex to one dimension higher.  Note that prisms are usually nonsimplicial; for example, a 3D triangular prism consists of two congruent triangles and three quadrilaterals. 
 % faces.  
% ~\cite{guo2021ftk}.

\textbf{Polytopes and their triangulation}. An \textit{$n$-polytope} is a geometry object in $\mathbb{R}^n$ whose faces are all \textit{flats} that can be described by a system of linear equations~\cite{sommerville2020introduction}. A triangulation of an $n$-polytope is a subdivision of the polytope into a finite collection of $n$-simplices such that the union of all simplices is equal to the polytope (\textit{union property}) and intersection of any pair of simplices is a common $(n-1)$-facet or empty (\textit{intersection property})~\cite{de2010triangulations}. We focus on 3-polytopes, namely,  \textit{polyhedra}, that may or may not be convex in this study. Triangulation of a polyhedron is also referred to as \textit{tetrahedralization} in the following paragraphs.

\textbf{Triangulation problems related to this research}.  We relate our research to two notoriously challenging triangulation problems: (1) triangulation without additional (Steiner) points~\cite{si2019simple} and (2) triangulation of nonconvex polytopes~\cite{chazelle1990triangulating, si2019simple}.  First, we intend to avoid Steiner points in spacetime because data are usually not immediately available beyond spatial mesh vertices.  Triangulation without Steiner points adds complexity and is not always achievable; for example, there exist nontriangulable polytopes (e.g., \textit{Sch\"onhardt polyhedron} shown in~\Cref{fig:staircase} (a))~\cite{de2010triangulations}, which can be created by rotating one of two parallel equilateral faces of a 3D triangular prism and inserting opposite diagonals in previous rectangles.  Second, we must handle nonconvex polytopes introduced by a deforming spacetime, whereas most triangulation focus on convex polytopes (i.e., convex hulls of sets of finite 3D points) only~\cite{below2004complexity, de2010triangulations}.  Nonconvexity causes severe problems because not all newly created connections remain in the polytope.  Also, the time complexity of triangulating a nonconvex polytope is high ($O(|\mathcal{V}|^3)$~\cite{chazelle1990triangulating, si2019simple}); we must incorporate even more complex constraints in this research, as described in the next section. 

\textbf{Spacetime triangulation for computational sciences}.  Scientists have recently explored using 4D spacetime meshes to simulate partial differential equations directly without using traditional 3D meshes and timestepping methods~\cite{erickson2005building, ishii2019solving}.  Although related, the goal of this paper is fundamentally different from 4D spacetime meshing in computational sciences; we attempt to represent science variables available on spatial grids, which are still the prevailing way to represent field data in today's computational sciences. While native spacetime meshes may interpolate and represent time-varying variables as conventional meshes, spacetime simulations are subject to high complexity, increased memory footprint, and numerical instabilities.  
%GAIL - you had "variables than" which did not make sense -- I substituted "as" - please check
% 
% Spacetime meshing in computational science serves for finite element method (FEM), a class of numerical methods based on discretization of the domain that are frequently used for building partial differential equation (PDE) solvers in mathematically modeling~\cite{erickson2005building, ishii2019solving}. 
For example, in spacetime finite element methods, one can create and update an $(n+1)$-D spacetime mesh by adding $(n+1)$-simplices on an $n$-D spatial mesh along the time dimension~\cite{thite2009adaptive}. %, just like building a wall by adding bricks on foundation along height dimension.
% while considering adaptivity of the algorithm, i.e., the ability to update mesh for varying phenomena and guarantee that the updated mesh is still suitable for PDE solvers in the meanwhile~\cite{abedi2004spacetime, thite2009adaptive}.
This method forms a simplicial complex with a terrain-like surface as spacetime. Researchers have expressed interest in different phenomena, such as wave propagation~\cite{thite2009adaptive} and rotation~\cite{von2021four}, that form spacetimes with different shapes (e.g., cone, prism) or time-variant topology of the spatial mesh.

\textbf{Spacetime triangulation in visualization} has focused on nondeforming spacetime so far, primarily for feature tracking problems.  Given a time-invariant spatial discretization, for example, a triangular mesh, one can establish a prismatic mesh connecting all corresponding nodes in spacetime.  One can further tesselate the prismatic mesh with the staircase triangulation scheme~\cite{de2010triangulations, guo2021ftk}, which we  review in the next section.

% Spacetime triangulation in visualization has been studied for non-deforming spacetime based on a staircase triangulation scheme~\cite{guo2021ftk, de2010triangulations}, which is further reviewed in the next section.


% Although there have been some methods to triangulate polyhedra, 3D deforming spacetime meshing still needs more effort. 

% First, 
% most of triangulation methods focus on convex polyhedra (i.e., convex hulls of sets of finite 3D points) only~\cite{below2004complexity, de2010triangulations}. Nonconvexity causes serious problems because not all newly created connections still stay in the polyhedron. Second, some works consider nonconvex polydehra but meet with time complexity issue (e.g., $O(n^3)$)~\cite{chazelle1990triangulating, si2019simple}. More importantly, unlike the spacetime generated by connecting nodes in two spatial meshes, a polyhedron in geometry do not have any inner edges on its polygonal faces or inside of it as constraints of triangulation. It is possible to exhaust all valid triangulations to find those satisfying all the constraints in spacetime, but it even increases time complexity.


% \subsection{Spacetime Meshing for Computational Sciences}



% Besides feature tracking, spacetime meshing are also used for other computational and visualization tasks. This section summarizes several problems that spacetime meshing can deal with in computational science and scientific visualization domains.

% \HG{Again, need lead sentence(s) to mention what aspects are related and why.}

% % \textit{Computational science}. 
% Spacetime meshing in computational science serves for finite element method (FEM), a class of numerical methods based on discretization of the domain that are frequently used for building partial differential equation (PDE) solvers in mathematically modeling~\cite{erickson2005building, ishii2019solving}. A typical way to create and update an $(n+1)$-D spacetime mesh is adding $(n+1)$-simplices on an $n$-D spatial mesh along time dimension~\cite{thite2009adaptive}, just like building a wall by adding bricks on foundation along height dimension.
% % while considering adaptivity of the algorithm, i.e., the ability to update mesh for varying phenomena and guarantee that the updated mesh is still suitable for PDE solvers in the meanwhile~\cite{abedi2004spacetime, thite2009adaptive}.
% This method forms a simplicial complex with a terrain-like surface as spacetime. There have been interests in different phenomena, such as wave propagation~\cite{thite2009adaptive} and rotation~\cite{von2021four}, which form spacetimes with different shapes (e.g., cone, prism) or time-variant topology of spatial mesh.
% % Traditional FEMs are based on time-stepping approach that builds spatial meshes and solves differential equations by marching sequentially in time~\cite{erickson2005building}. These methods achieve temporal parallelism but limit parallel scalability over space~\cite{ishii2019solving}. Spacetime discretization treats spatial and temporal dimensions as a whole and is therefore applied more recently to overcome this limitation~\cite{ishii2019solving, owrie1998space, thite2009adaptive}.

% \textit{Scientific visualization}. Spacetime mesh enables tasks including not only feature tracking as mentioned in~\Cref{sec:feature_tracking}, but also upsampling, volume rendering, and isosurfacing. With spacetime mesh, upsampling can be achieved by interpolation in cells. Volume data can be rendered by many algorithms such as ray-tracing and tools such as Paraview if it is sampled on mesh, no matter if the mesh is structured, unstructured, or from Adaptive Mesh Refinement (AMR)~\cite{sarton2023state}. The performance of our method will also be evaluated via these applications later.

% \HG{It may be useful to review staircase triangulation for explaining the ``straight'' extrusion.}

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.98\columnwidth]{figs/staircase.pdf}
  \caption{
  % \HG{It may be helpful to illustrate how each monotonous path correspond to a tetrahedron; see Fig. 4 in the FTK paper.  Although staircase is not new, it is necessary to introduce sufficient background for readers to understand.} 
  Subdivision of a 3D prism. (a) A Sch\"onhardt polyhedron, which is the simplest nonconvex nontriangulable polyhedron. (b) Monotone paths, corresponding tetrahedra, and simplicial subdivisions of a 3D prism by staircase triangulation, with the assumption that the 3 node indices satisfy $a_0<a_1<a_2$. Nodes in the lower mesh are only connected to nodes in upper mesh whose previous nodes have greater indices. (c) A simplicial subdivision of the 3D prism without considering staircase. (d) If all  3 quads are separated by connecting top-left corner and bottom-right corner, then there is no way to triangulate the prism without additional nodes. (e) Two adjacent prisms have conflictive triangulation ways on the common quadrilateral face $a_1-a_2-b_2-b_1$.
  % \HG{If you plan to explain Sch\"onhardt, it may be helpful to have an illustration of a Sch\"onhardt prism.}
  }
  \label{fig:staircase}
\end{figure}



% Our research is closely related to computational geometry, specifically triangulation. Starting with reviewing terms, we point out the similarity and difference between this work and previous triangulation methods in geometry domain.

% \HG{I miss a lead sentence here; e.g., our research is closely related to computational geometry, specifically triangulation.}



% \HG{Also, how is our work different from existing triangulation research?  What are unique/special in our research?}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Formulation and Preliminaries}
\label{sec:preliminaries}

% This section reviews necessary background for formalizing the deformed spacetime meshing. % problem and the driver application.  % rudimentary knowledge and We start with a brief introduction to basic conceptions in computational geometry and then focus on \textit{staircase triangulation}, a triangulation rule of prism that highly related to and is used in our work. We also introduce the fusion energy application that drives our work.

% \subsection{Deforming Spacetime}

% \textbf{Deformation induced by a vector field}.
This paper considers the deforming spacetime induced by a vector field (e.g., magnetic field in fusion plasma) $\mathbf{v}: \mathbb{R}^{n+1}\to\mathbb{R}^{n}$ (or $\mathbb{R}^{n}\to\mathbb{R}^{n}$ if the deformation is time-invariant), where $n$ is the dimensionality of the spatial domain.  Assuming Lipschitz continuity of $\mathbf{v}$, the \emph{deformation} $\Phi: \mathbb{R}^{n+2}\to\mathbb{R}^{n}$ (also known as \emph{flow map}) is a continuous function representing the solution of the following initial value problem: 
% 
\begin{equation}
  \frac{\partial\Phi(\mathbf{x}_0, t_0, t)}{\partial t} = \mathbf{v}(\Phi(\mathbf{x}_0, t_0, t))~\mbox{and}~\Phi(\mathbf{x}_0, t_0, t_0)=\mathbf{x}_0,
\end{equation}
% 
where $\mathbf{x}_0$ is a spatial location at time $t_0$ and $\Phi(\mathbf{x}_0, t_0, t)$ denotes the spatial location of $\mathbf{x}_0$ after deformation at $t$. With the deformation, one can define the field-following interpolation as 
\begin{equation}\label{eq:ff}
  f(\mathbf{x}, t) = \beta f(\Phi(\mathbf{x}, t, t_0)) + (1-\beta) f(\Phi(\mathbf{x}, t, t_1)),
\end{equation}
where $f: \mathbb{R}^{n+1}\to\mathbb{R}$ is a scalar function (e.g., density and temperature in fusion plasma) and the interpolation weight $\beta$ is $(t_1 - t) / (t_1 - t_0)$, $t_0\leq t\leq t_1$. 

In a discrete sense, we rephrase the objective of spacetime meshing in Section~\ref{sec:intro} as follows. Assuming each timestep $t_i$ in the sequence $t_0\le t_1 \le \ldots t_i \le \ldots$ is associated with an $n$-dimensional simplicial mesh $M_i=\left<V_i, E_i\right>$, where $V_i=\{v_0^{(i)},v_1^{(i)},...\}$ is the set of nodes in $M_i$ and $E_i=\{e_0^{(i)},e_1^{(i)},...\}$ is the set of edges in $M_i$, $v_j^{(i)}$ denotes the $j$th node and $e_j^{(i)}$ denotes the $j$th edge in $M_i$. A deformed spacetime mesh is an $(n+1)$-dimensional simplicial mesh $\mathcal{M}=\left<\mathcal{V}, \mathcal{E}\right>$ that contains all nodes $\mathcal{V}=\cup_i V_i$ and edges $\mathcal{E}=\left(\cup_i E_i\right) \cup \left(\cup_i\cup_{j,k} v^{(i)}_j v^{(i+1)}_k \right)$.  The temporal connectivity $v^{(i)}_j v^{(i+1)}_k$ between $M_i$ and $M_{i+1}$ is defined by the deformation $\Phi$ such that $v^{(i+1)}_k \approx \Phi(v^{(i)}_j, t_i, t_{i+1})$.  In practice, we only need to consider spatial meshes in two consecutive timesteps (namely, the \emph{lower} and \emph{upper meshes}); spacetime meshing can be trivially applied to every two adjacent timesteps for multiple timesteps.

In the rest of this section we first review the trivial case, triangulation of nondeforming spacetime ($\mathbf{v}=\mathbf{0}$), and then consider challenges in deforming spacetime.



% \subsection{Staircase Triangulation: A Method to Mesh the Non-Deforming Spacetime}
\subsection{Meshing Nondeforming Spacetime with Staircase Triangulation}
% \subsection{Non-Deforming Spacetime}
\label{sec:staircase}

With nondeforming spacetime, assuming the underlying spatial mesh $M_i$ is identical for all timesteps, the problem is reduced to the following: each node in the lower mesh connects to the same node in the upper mesh.  As a result, each lower-mesh triangle extrudes into a prism, which can be further subdivided into three tetrahedra with \emph{staircase triangulation}~\cite{de2010triangulations}, as illustrated in~\Cref{fig:staircase} (b).
% 
As such, staircase triangulation can help mesh $(n+1)$-dimensional nondeforming spacetime~\cite{guo2021ftk} and serve as a basis of our methodology. % and is further generalized for meshing the deforming spacetime.  
% We review key definitions and concepts below.

% the spacetime meshing can be achieved with \textit{staircase triangulation}~\cite{guo2021ftk, de2010triangulations}.  

% rudimentary knowledge and We start with a brief introduction to basic conceptions in 
% We start with \textit{staircase triangulation}, a general scheme to triangulate prismal cells~\cite{de2010triangulations} and can help mesh $(n+1)$-dimensional non-deforming spacetime~\cite{guo2021ftk}.  Staircase triangulation serves as a basis of our methodology and is further generalized for meshing the deforming spacetime.  We review key definitions and concepts below.

% triangulation rule of prism that highly related to and is used in our work. We also introduce the fusion energy application that drives our work.

\textbf{Staircase triangulation} provides an arbitrary rule to triangulate a prism based on node indices. Consider a prism $a_0a_1a_2-b_0b_1b_2$ extruded by a triangle $a_0a_1a_2$. With staircase triangulation, each tetrahedron corresponds to a \textit{monotone path} where both alphabets and subscripts are ascending. For example, all monotone paths for prism $a_0a_1a_2-b_0b_1b_2$ are $a_0a_1a_2b_2$, $a_0a_1b_1b_2$, and $a_0b_0b_1b_2$, as shown in~\Cref{fig:staircase}~{(b)} left.  Each path represents a tetrahedron separated from original prism by staircase triangulation (\Cref{fig:staircase}~{(b)} middle), and all tetrahedra corresponding to monotone paths constitute a tetrahedralization of the prism (\Cref{fig:staircase}~{(b)} right).  


% Note that a prism may be triangulated in multiple ways.  Consider the triangular prism with node indices of $0$, $1$, and $2$ as $a_0a_1a_2-b_0b_1b_2$; the prism can be divided into 3 tetrahedra by separating the 3 lateral faces of the prism into 6 triangles, and thus there are $2^3=8$ possible ways to subdivide 3 lateral quadrilaterals in total. 
% % Consider 3D triangular prisms that are related to 3D deforming spacetime. Every prism can be divided into 3 tetrahedra by separating the 3 lateral faces of the prism into 6 triangles, and thus there are $2^3=8$ separations of the 3 lateral quadrilaterals in total. 
% However, not all schemes can divide a prism into tetrahedra~(\Cref{fig:staircase}~{(d)}). 

% We can see that each node is immediately above or to the right of the previous node on the path in staircase triangulation~\cite{de2010triangulations}.

% For example, if we denote nodes on lower face by $a_0$, $a_1$, $a_2$ and their adjacent nodes on upper face by $b_0$, $b_1$, $b_2$,

% Tetrahedra derived by staircase triangulation can be described by \textit{monotone paths}  if we denote nodes on lower face by $a_0$, $a_1$, $a_2$ and their adjacent nodes on upper face by $b_0$, $b_1$, $b_2$, then all monotone paths for the 3D prism are $a_0a_1a_2b_2$, $a_0a_1b_1b_2$, and $a_0b_0b_1b_2$, where both alphabets and subscripts are ascending, as shown in~\Cref{fig:staircase}~{(e)} left. Each path represents a tetrahedron separated from original prism by staircase triangulation (\Cref{fig:staircase}~{(e)} middle).  We can see that each node is immediately above or to the right of the previous node on the path in staircase triangulation~\cite{de2010triangulations}.


\textbf{Meshing the nondeforming spacetime}.  Staircase triangulation directly generalizes to multiple prisms and makes it possible to triangulate prismatic meshes (\Cref{fig:spacetime} (a)) extruded from a 2D spatial mesh into 3D. % and involves multiple adjacent prisms.  
Note that although staircase triangulation is not the only way to triangulate each prism (see \Cref{fig:staircase} (c) as an example), indices that give a unique ordering of nodes define an elegant and consistent schema to represent a nondeforming spacetime as a simplicial complex.
% With the node indices available from 2D spatial meshes, one can triangulate each prism independently.  
For example, for the two adjacent prisms in~\Cref{fig:staircase}~{(e)}, both prisms triangulate the shared quadrilateral face in a consistent manner by staircase, creating a shared edge $a_1b_2$ that connects the lower index node ($a_1$) to a higher index node ($b_2$); the resulting triangulation includes no conflicting edge, such as $b_2a_1$.  

% there might be a conflict on their common quadrilateral face in prismatic mesh if we separate the quadrilateral in different ways when triangulating the two prisms containing it.

% Even if every prism is tetrahedralized successfully, for two adjacent prisms, there might be a conflict on their common quadrilateral face in prismatic mesh if we separate the quadrilateral in different ways when triangulating the two prisms containing it~(\Cref{fig:staircase}~{(e)}). To avoid these issues in simplicial subdivision for prismatic mesh, we apply staircase triangulation to coordinate subdivision on common faces, where connections between only nodes in lower mesh to nodes in upper mesh whose previous nodes has greater indices are allowed in both two prisms.



% \subsection{Problem Formulation}
% We formulate the problem of meshing the deforming spacetime introduced in~\Cref{sec:xgc} into a general simplicial meshing problem for 3D deformed spacetime generated by extruding 2D mesh along time dimension.
% % The triangulation in the 2D spatial mesh is invariant over time,
% % \HG{Again, do we have to assume the 2D mesh is invariant? We need to keep the statement consistent in the paper}
% Time-dependent physical behavior introduces twisting connectivity
% % \HG{change wording; avoid using topology}
% between consecutive timesteps, and thus deforms the shape of spacetime.

\subsection{2.5D Representation of a Deforming Spacetime: Lower/Upper Meshes and Temporal Connectivity}
% \subsection{Defining Deforming Spacetime with Lower/Upper Spatial Meshes and Temporal Connectivity}
% \subsection{Deforming Spacetime}
\label{sec:spacetime_formulation}



% \HG{I am not sure I understand what the subsection title means.}
% How a spacetime is induced from a lower-D spatial mesh depends on how next node is defined. 

% Without loss of generality, one can define an $(n+1)$-D deforming spacetime with (1) two $n$-D time-varying spatial meshes, namely the \textit{lower} and \textit{upper meshes} representing two consecutive timesteps, $t$ and $t+1$, respectively, and (2) connectivity of nodes across the two meshes.  We assume that both upper and lower meshes are simplicial, for example, 2D meshes containing only triangular cells.  

% This paper considers the temporal connectivity defined by \textit{next node}, that is, 
% Formally, how to build the correspondences between nodes in two successive timesteps depends on transporting behaviors of underlying physics. 

\begin{figure}[!t]
    \centering
    \includegraphics[width=\columnwidth]{figs/XGC_interpolation.pdf}
    \caption{Illustration of problem formulation. (a) Magnetic-following interpolation. It traces magnetic lines (shown as blue line) passing a given location $p$ at time $t$ in both forward and backward directions and locates intersections on two neighboring timesteps at $t_i$ and $t_{i+1}$. The function value at $(p,t)$ is interpolated by function values at two intersecting points. Gray lines show the magnetic lines passing nodes in $t_i$. (b) Next nodes (connected by black lines) of nodes in $t_i$ derived by magnetic lines and a triangulation based on next nodes (gray lines). Only connections whose two ending points are both in the 1D region are shown. 
    % \HG{To help explaining next-node in 4.1, I would split this figure into two (or optionally three) subfigures: (a) magnetic lines (draw a few more gray streamline as a context), (b) the next node of each node in $\phi_i$, (c) a triangulation based on next node; keep $p$ in each subfigure.}
    }
    \label{fig:xgc_interp}
\end{figure}

This paper considers the temporal connectivity defined by \textit{next node}; that is, each node in the lower mesh connects to another node (or null) in the upper mesh.  While alternative definitions may be valid (e.g., one-to-many) for our algorithms, we exclusively use the next-node definition because of the simplicity of deriving deformation induced by the vector field $\mathbf{v}$. 
% flow transport in the deformation induced by a vector field.
For example, in tokamaks where the magnetic field defines the deformation, 
% particles move under the influence of magnetic field, 
by placing a particle at a mesh node $v$ in a lower mesh, scientists can build a connection between $v$ and the node $v^\prime$ (called \textit{next node} of $v$) closest to the particle's position in the upper mesh (\Cref{fig:xgc_interp} (b)). 
% $n$ is also said to be \textit{previous node} of $n'$. 
% The way in which we build node connectivity tries to ensure that the same spatial features or objects at different timesteps are sewed together by the connectivity. 
% This way to build connectivity can also be easily generalized to spacetime construction for other applications with time-varying spatial quantities, although we use tokamaks as an example.

% This paper considers partial connectivity defined by \emph{next node}. Formally, for an $n$-D time-varying spatial mesh, its $(n+1)$-D spacetime is generated by connecting every node in lower mesh to another node in upper mesh. The correspondences between nodes in two successive timesteps are given by transporting behaviors of underlying physics, such as the movement of particles caused by magnetic field in tokamaks, where we build a connection between a node $n$ in lower mesh and the node $n'$ (called \textit{next node} of $n$) closest to its position in next timestep $t+1$ in upper mesh (\Cref{fig:xgc_interp} (b)). $n$ is also said to be \textit{previous node} of $n'$.

% The mapping from node to next node is neither injective nor surjective. % By our definition of next node for XGC data, the mapping is injective but not surjective. 
% In other words, each node in lower mesh has exactly one next node, but different nodes can have the same node in upper mesh as their next node. Nodes in upper mesh might have zero, one, or multiple previous nodes. As a result of the mapping, a $k$-simplex ($k=0,1,2$) in lower mesh may correspond to a simplex in a different dimension (e.g., triangle-to-node, triangle-to-edge, edge-to-edge, and edge-to-node) or a simplicial complex in upper mesh.  % , which will be further discussed below. % in the next section.
% % Intuitively, the correspondence reflects local deformation in the spacetime.
% % \HG{I would add a few sentences to explain the mapping of next-node is injective (\emph{not} bijective); also explain the reason.  A 1D illustration of the injection may be helpful as well.}

% \HG{We probably need some explanation on why fwe choose to study connectivity defined by next nodes.  Availability in XGC is a reason, but is this general enough?  There may exist other ways to define the connectivity.}

% \HG{I guess the connectivity may be defined in many different form, not necessarily specific to XGC's way of giving a injective mapping from lower plane nodes.}


% \subsubsection{Local Spacetime Deformed by Next Nodes} % A Speculation into Seven Scenarios in 3D Spacetime}
% \label{sec:seven_scenarios}

% 
% \textbf{Local spacetime deformed by next nodes}. 


\textbf{2.5D independent partition}. Denote the lower and upper spatial meshes by $S$ and $S^\prime$, respectively. We define a \emph{2.5D independent partition} $PQ$, where $P$ and $Q$, respectively, is a simply connected component in the lower and upper meshes such that (1) any node $v\in P$ maps to another node $v^\prime\in Q$ and (2) any node $v\in \partial P\setminus\partial S$ maps to another node $v^\prime\in \partial Q\setminus\partial S^\prime$; the symbol $\partial$ denotes the boundary of a partition.  Note that $P$ or $Q$ may be degenerate, such as a single node or one or multiple edges.  For example, the entire lower and upper meshes are a 2.5D independent partition. For another example, assuming $P$ contains one single triangle, the next-node mapping implicitly determines diverse situations on how $P$ maps to upper-mesh elements, which could be a single node, one or multiple edges, or one or multiple triangles as finite scenarios illustrated in~\Cref{fig:seven_scenarios}. 
The seven scenarios indicate different local deformation behaviors in spacetime. If we consider the deformation in the scale of edge, then two adjacent nodes may still be mapped to two adjacent nodes (small change of relative distance between the two nodes, e.g., $n_4$ and $n_5$ in~\Cref{fig:xgc_interp} (b)), merge to one node (shrinking behaviors, e.g., $n_5$ and $n_6$ in~\Cref{fig:xgc_interp} (b)), or break into a path (expanding behaviors, e.g., $n_7$ and $n_8$ in~\Cref{fig:xgc_interp} (b)) in the upper mesh.  

\textbf{2.5D domain decomposition}.  A decomposition of $PQ=\cup_i P_iQ_i$ is defined by nonoverlapping 2.5D independent partitions $P_iQ_i$.  Formally, for any two different partitions $P_iQ_i$ and $P_jQ_j$, the intersection is either null or lower-dimensional simplicies such as nodes and edges.  The union and intersection of two 2.5D independent partitions $P_iQ_i$ and $P_jQ_j$ are defined by $(P_i\cup P_j)(Q_i\cup Q_j)$ and $(P_i\cap P_j)(Q_i\cap Q_j)$, respectively.  We will further discuss how the 2.5D domain is decomposed into smaller partitions for divide-and-conquer processing in the next section.


% Two \textit{static nodes}, i.e., two adjacent nodes that are still mapped to two adjacent nodes in next timestep, represent small change of relative distance between the two nodes, such as $n_4$ and $n_5$ in~\Cref{fig:xgc_interp} (b).
% \HG{Have we defined the velocity of a node?}. 
% \HG{To help readers memorize and find out what are static nodes, what about using illustration?  A good place may be Fig.3(b).}
% Two \textit{merged nodes} (two adjacent nodes that merge to one node, such as $n_5$ and $n_6$ in~\Cref{fig:xgc_interp} (b)) mean congregation of particles around that location, while two \textit{split nodes} (e.g., $n_7$ and $n_8$ in~\Cref{fig:xgc_interp} (b)) that break into a path tell diffusing behavior of particles.

% and it seems easy to triangulate it by staircase triangulation. How we apply staircase triangulation to meshing deforming spacetime will be discussed in~\Cref{sec:face_triangulation}.
% Every \textit{triangle-to-node} cell directly gives us a tetrahedron. Tetrahedralization of \textit{triangle-to-triangle} is also solved and will be discussed in next section. Other scenarios need further discussion.

% Conceptually, the mapping implies diverse behaviors, including shrinking/expansion, rotation, and translation, as further explained below.

% In deforming spacetime, injective node-to-next-node mapping gives triangle cells in lower mesh different fates in next timestep. If every node in lower mesh is connected to itself in upper mesh, then every 2D triangle cell forms a triangle-to-triangle correspondence, which gives us a prismatic mesh for spacetime~(\Cref{fig:spacetime} (a)). But due to the variant property, a node $n$'s next node may not still be $n$ in upper mesh, which not only causes deforming shape of spacetime but also makes it impossible to find corresponding cells in upper mesh for cells in lower mesh~(\Cref{fig:spacetime} (c)). Two adjacent nodes in lower mesh might still form an edge, degenerate into a node, or become non-adjacent (i.e., form a \textit{path}) in upper mesh, which leads to 7 \textit{triangle-to-what} scenarios, as shown in~\Cref{fig:seven_scenarios}.


\begin{figure}[!t]
    \centering
    \includegraphics[width=0.49\textwidth]{figs/seven_scenarios.pdf}
    \caption{Seven \textit{triangle-to-what} scenarios. Connected components in lower and upper meshes are linked by next-node correspondence. Adjacent nodes might still be adjacent (in \textit{triange-to-triangle} scenario), merge (in \textit{triange-to-edge}, \textit{triange-to-node}, and \textit{triange-to-path} scenarios) or become non-adjacent (in \textit{triange-to-2edge1path}, \textit{triange-to-1edge2path}, and \textit{triange-to-3path} scenarios) in the next  timestep.}
    \label{fig:seven_scenarios}
\end{figure}

% \subsection{Assumptions}
% \label{sec:assumption}


\textbf{Triangulation of a 2.5D independent partition}.  A 2.5D independent partition $PQ$ consists of only triangles, edges, and nodes from 2D complexes $P$ and $Q$ and does not bound a 3D volume; one must first transform $PQ$ into a proper polyhedron before triangulating into tetrahedra.  We refer to deriving a polyhedron from a 2.5D independent partition as \emph{lateral triangulation}.  To form a polyhedron, because $P$ and $Q$ are already triangles or degenerate triangles, one has to define a surface/triangular mesh that connects 2D boundaries $\partial P$ and $\partial Q$; the lateral surface must contain edges between a node $v\in \partial P$ and its next node $v^\prime\in \partial Q$. However, multiple choices exist to find a lateral triangulation.  For example, two ways exist to triangulate a quadrilateral face on a prism.  The non-uniqueness leads to challenges: (1) the resulting polyhedron may be nontriangulable, as discussed below, and (2) two 2.5D independent partitions should have the same triangulation on their common lateral face, or otherwise one polyhedron will overlap or form a void with its neighboring polyhedron, as addressed in the next section.

\textbf{Triangulability of a 2.5D independent partition}.  Some partitions cannot find a triangulation, which we call \emph{ill-posed partitions}.  For example, the prism in~\Cref{tab:ill_cases} is twisted so much that any lateral triangulation will lead to self-intersections.  Although the prism is a triangle-to-triangle case, the ill-posed prism does not legitimately enclose a volume.  Ill-posed partitions are separate from  \emph{quasi-ill-posed partitions}, which can triangulate with a different lateral triangulation.  For example, a prism has eight possible ways to tessellate three quadrilateral faces (\Cref{tab:ill_cases}); two of them lead to a Sch\"onhardt prism and cannot triangulate, but choosing a different scheme will avoid this situation. We further discuss the adjustment of lateral triangulation later.

Should an ill-posed partition appear while no alternative way exists to decompose the domain, the spacetime cannot be triangulated. Unfortunately, ill-posed partitions suggest that temporal discretization is insufficient. In this case one can refine the temporal resolution with a smaller temporal gap $\Delta t$ so that no ill-posed partitions appear.\footnotemark~In other words, one can upsample the temporal resolution with field-following interpolation (Eq.~\eqref{eq:ff}) before attempting spacetime meshing.
\footnotetext{In general, the distortion decreases with a finer temporal resolution $\Delta t$; the distortion characterized by displacement $||\Phi(\mathbf{x}, t_0, t) - \mathbf{x}||$ is bounded by $L\cdot \Delta t$, $L$ being the Lipschitz constant of $\mathbf{v}$.} 

% ; one can reduce the distortion with a smaller $\Delta t$.}

% So far, our only assumption is the Lipschitz continuity of $\mathbf{v}$, and we must further assume any local spacetime is triangulatable.  As a counterexample, the prism in \HG{Figure X}, which we call \emph{ill-posed prism}\footnotemark, is non-triangulable.  Although the prism seems to be a triangle-to-triangle case, the prism is deformed so much that the prism does not legitimately enclose a volume.  
% \footnotetext{Ill-posed prism is not to confuse with \emph{quasi-ill-posed prisms}, as defined later, which are artifacts introduced by staircase triangulation and can resolve.}
% 
% Ill-posed cases suggest that temporal discretization is insufficient.  Luckily, one can find a sufficiently small temporal gap $\Delta t$ such that no ill-posed cases appear.\footnotemark~In other words, one can upsample the temporal resolution with field-following interpolation (Eq.~\eqref{eq:ff}) before attempting spacetime meshing.
% 
% \footnotetext{Intuitively, the distortion characterized by displacement $||\Phi(\mathbf{x}, t_0, t) - \mathbf{x}||$ is bounded by $L\cdot \Delta t$, $L$ being the Lipschitz constant of $\mathbf{v}$; one can reduce the distortion with a smaller $\Delta t$.}



% seems straightforward with staircase triangulation, the deformed prism is twisted
% a non-triangulatable prism is illustrated in \HG{Figure X}.  




% We make the following assumptions for the deforming spacetime to be tetrahedralized by our algorithm.

% The difficulty triangulating the spacetime formed by a single cell also varies by scenario. For example, every \textit{triangle-to-triangle} scenario creates a \textit{deformed prism}, which seems straightforward with staircase triangulation but is not always applicable, as discussed in later sections.

% \begin{itemize}
%     \item Time discretization of the data is fine enough such that it yields only triangulable things. If not, upsampling the data over temporal scale is needed. See~\Cref{sec:face_triangulation} for more details.
%     \item The spatial meshes at different timesteps may or may not be identical.
%     \item The mapping from node to next node could be neither injective nor surjective.
% \end{itemize}

\begin{figure*}[!t]
    \centering
    \includegraphics[width=\textwidth]{figs/pipeline.pdf}
    \caption{Pipeline of our tetrahedralization algorithm. The algorithm goes through two stages, \textit{divide} and \textit{conquer}. In the divide stage, the input is a 3D spacetime formed by lower mesh, upper mesh, and connectivities between them, and the output is a set of 2.5D independent partitions. Then we triangulate their lateral faces. The conquer stage take every face-triangulated polyhedron as an input. Note that the original 3D spacetime is also a 2.5D independent partition and one can skip the divide stage without complexity consideration. Indivisible prisms are filtered before dividing them, while more complex polyhedra need tentatively dividing first and retriangulating lateral faces if they fail to be divided. The output of the conquer stage as well as the whole algorithm is a set of tetrahedra that form a simplicial mesh of the input 3D spacetime.}
    \label{fig:pipeline}
\end{figure*}

\section{Methodology}
\label{sec:method}

%, as discussed in Section~\ref{sec:spacetime_formulation}.  
% (so called \emph{cutting faces} in this paper).  
%, which are also illustrated in \HG{Figure X}:

\Cref{fig:pipeline} illustrates our workflow to triangulate a 3D deforming spacetime with a \textit{divide-and-conquer} strategy, with the rationale that the domain may be decomposed into 2.5D independent partitions (or simply partitions when there is no ambiguity), each of which could lead to a polyhedron for further triangulation. Multitude challenges exist, including (1) defining a proper domain decomposition that could lead to triangulatable partitions, (2) maintaining compatibility, that is, ensuring no overlaps or voids exist between neighboring partitions, (3) finding a proper lateral triangulation scheme that leads to a possible 3D triangulation, and (4) triangulating a nonconvex polyhedron.  Each challenge is nontrivial and requires careful design.  To these ends, we design our methodology as three major phases with multiple refinements and trial-and-error steps: 
\begin{itemize}\setlength{\parskip}{0pt}
  \item \textbf{Spacetime decomposition} (divide stage): Decomposing the domain into as many  partitions as possible based on \emph{closed cutting paths} and inner-node refinements (\Cref{sec:divide}).
  % two spatial partitions that represent a local spacetime; %are matched by temporal connectivity of nodes; 
  % each spatial partition is a connected component of lower/upper mesh simplices separated by \emph{cutting faces} (\Cref{sec:divide}); 
  % may be degenerate, e.g., a path or a node (\Cref{sec:divide});
  \item \textbf{Lateral triangulation}: For each partition, forming a polyhedron by meshing the surface that connects the partition's lower and upper boundaries.  By default, staircase triangulation rules are applied when possible but subject to change if later volume triangulation is impossible (\Cref{sec:ftp}).
  
  % \item \textbf{Face-triangulated partition}: a spacetime polytope with all lateral faces triangulated (\Cref{sec:ftp});
  \item \textbf{Volume triangulation} (conquer stage): For each polyhedron, searching for a possible triangulation scheme with a two-tiered decision-tree algorithm (\Cref{sec:conquer}).
\end{itemize}


%  multiple refinements and trial-and-error steps to guarantee compatibility so that one can further handle partitions independently.
% To ease description, we define terminologies below:
% \begin{itemize}\setlength{\parskip}{0pt}
%   % \item \textbf{Space partition (ST)}: a connected component of spatial simplices separated by \emph{cutting faces} (defined in \Cref{sec:cutting_faces}) in lower/upper mesh;
%   \item \textbf{Spacetime partition}, or simply \textbf{partition}: two spatial partitions that represent a local spacetime; %are matched by temporal connectivity of nodes; 
%   each spatial partition is a connected component of lower/upper mesh simplices separated by \emph{cutting faces} (\Cref{sec:divide}); 
%   % may be degenerate, e.g., a path or a node (\Cref{sec:divide});
%   \item \textbf{Face-triangulated partition}: a spacetime polytope with all lateral faces triangulated (\Cref{sec:ftp});
%   \item \textbf{Volume-triangulated partition}: resulting tetrahedra of the face-triangulated spacetime partition (\Cref{sec:conquer}).
% \end{itemize}




% A key difficulty in our workflow is maintaining compatibility between partitions; that is, the shared spacetime boundary, namely \emph{lateral faces}, between two neighboring partitions must share the same triangulation.  




% We subdivide the 3D spacetime into 3-simplices by a \textit{divide-and-conquer} paradigm (\Cref{fig:pipeline}). 
% In our workflow, the ``divide'' stage splits the spacetime into smaller and independent partitions; each representing a local spacetime and sharing boundary with neighboring partitions (\Cref{sec:divide}).

% % ---such that triangulation in one partition does not influence any other partitions 
% % 
% Once the spacetime decomposes into partitions, the ``conquer'' stage triangulate each partition with a search-tree algorithm (Section~\ref{sec:conquer}).
% 





% \HG {The following sentence is too long to understand.}
% In our workflow, the \textit{divide} stage splits the spacetime into smaller \textit{partitions} that are polyhedra whose edges are all elements in original 3D spacetime and volumes are subsets of original 3D spacetime, such that tetrahedralization in one partition does not influence that in all other partitions, while the \textit{Conquer} part tetrahedralizes every partition by continuously eliminating mesh nodes and edges incident on it. It is obvious that the original spacetime is also a partition and can directly be input of Conquer stage, but Divide stage largely reduces the complexity of searching by decision tree in Conquer stage.
% \HG{Can you formally define what is a piece?}

% \HG{The above is a good lead paragraph of what you do.  Can you explain why do we introduce a divide-and-conquer method?  Because of complexity or what?}
% \HG{From reviewers' eye, deeper questions may be: (1) is divide-and-conquer guaranteed to find a feasible solution? (2) is the solution unique? (3) if not, is the algorithm robust, and what solution does our method favor? (4) if our algorithm does not guarantee to find a solution, what should one do?}

\subsection{Spacetime Decomposition} % (the ``Divide'' Statge)}
\label{sec:divide}

% \HG{Similar to the comments I posted in 5.2, can you formally state what are the inputs and outputs for this step?}

% The divide stage 
% Spacetime partitioning aims to decompose the problem into smaller ones so that we do not need to triangulate the entire spacetime with our search tree.


With the given lower and upper meshes and temporal connectivity, the first step is partitioning the spacetime into nonoverlapping 2.5D independent partitions.  We aim to decompose the domain into much smaller partitions so that we do not need to triangulate the entire spacetime with our decision-tree algorithm. % in~\Cref{sec:conquer}.

% The spacetime partitioning is conducted in two phases to ensure these partitions are as small as possible so that we do not need to handle the entire spacetime with our search tree.  During the division, we need to (1) well choose cutting faces to guarantee that further subdivision of every piece is independent of those of other pieces; (2) coordinate the subdivision on cross sections to avoid conflicts on shared faces as shown in \Cref{fig:staircase}~{(e)}.

% satisfying both union and intersection properties as stated in~\Cref{sec:staircase}.
% Input of this stage is 3D spacetime derived in~\Cref{sec:spacetime_formulation}, and output is a set of 3D partitions satisfying both union and intersection properties as stated in~\Cref{sec:staircase}. 


\subsubsection{Domain Decomposition by Cutting Graphs}
% % \subsubsection{Dividing by cutting faces}
% \label{sec:cutting_faces}

We introduce the definition of \emph{cutting edges}, \emph{cutting graphs}, and \emph{cutting paths} and prove that the domain can be decomposed into 2.5D independent partitions (as defined in~\Cref{sec:spacetime_formulation}) as connected components isolated by cutting graphs.

\textbf{Cutting edges}.  Formally, with a partition $PQ$, edges $uv\in P$ and $u^\prime v^\prime$ form a pair of cutting edges if there exists edge $u^\prime v^\prime\in Q$, where $u^\prime$ and $v^\prime$ are the next node of $u$ and $v$, respectively.  % $u^\prime v^\prime$ is also referred to as a separating edge in the upper mesh.  
As shown in an example in \Cref{fig:divide}~{(a)}, a pair of cutting edges is a 4-node cycle formed by four edges: a lower-mesh edge $n_1n_9$, an upper-mesh edge $n_1^\prime n_9^\prime$, $n_1n_1^\prime$, and $n_9n_9^\prime$.  
\begin{restatable}{assumption}{myassumption}\label{myassumption}
The lower/upper mesh and time are sufficiently refined such that for any point $x$ on a cutting edge $v_0v_1$ in the lower mesh, the nearest node of $\Phi(x)$ on the upper mesh is either $v_0^\prime$ or $v_1^\prime$.
\end{restatable}
Here, we omit time in $\Phi$ for clarity.  Conceptually, we assume that a cutting edge reflects distinct transport behaviors of the underlying flow $\mathbf{v}$ that drives the deformation; for example, in a laminar flow that is parallel to the cutting edge, no particle shall cross the edge.  


% Noting that a cutting edge may connect to multiple cutting edges, we generalize cutting paths to \textbf{cutting graph}, a 1-connected component(s) of cutting edges.  
% For example, in the case that deforming spacetime is driven by a laminar flow, no particles cross particles ahead of them in previous timestep, and thus no connections between nodes and their next nodes cross any separating edges. No crossing connections guarantees tetrahedralizations in any two partitions do not interfere each other.

\textbf{Cutting graphs and cutting paths}. A cutting graph is a graph that connects nodes that are connected by cutting edges.  A cutting path is an arbitrary path in the cutting graph.  For example, path $u_0u_1u_2u_3$ is a cutting path if $u_0u_1$, $u_1u_2$, and $u_2u_3$ are cutting edges; the counterpart consisting of next nodes $u_0^\prime u_1^\prime u_2^\prime u_3^\prime$ is also referred to a cutting path in the upper mesh.  
% We further define \textbf{closed cutting path} if the path and its counterpart is closed (connecting to boundaries $\partial P$ and $\partial Q$ or forming a loop), and we have
%
\begin{restatable}{lem}{mylemma}\label{mylemma}
Partition $PQ$ is independent if the boundary $\partial P$ is a closed cutting path.
% Partition $PQ$ can be subdivided into two nonoverlapping 2.5D independent partitions by a cutting path if the path is closed (connecting boundaries of $PQ$ or forming a loop). 
\end{restatable}
\begin{proof}
In a trivial case, if all nodes on $P$ are on boundary $\partial P$, $PQ$ is independent because every node on $P$ maps to $Q$.  Now, assume $v_k$ is an inner node of $P$ such that $v_k\notin\partial P$ and the next node is $v_k^\prime\notin Q$.  Because the next node is the closest point to $\Phi(v_k)$, one can find a sufficiently small disk $B(v_k, \epsilon), \epsilon>0$ so that the flow map $\Phi$ of any point on the disk has the closest point $v_k^\prime$.  Now, the domain $P$ maps to at least two connected components by $\Phi$; one is the image of the disk, and another is the region that intersects $Q$.  Having two connected components contradicts Rudin's Theorom 4.22~\cite{rudin}; the continuous mapping $\Phi$ of $P$ shall be connected.  This terminates the proof.
\end{proof}
%
% \begin{proof}
% Denote the cutting path as $p$ and its upper-mesh counterpart as $q$.  Because both $p$ and $q$ are closed paths, each path will lead to two spatial complexes $P_0\cup P_1=P$ and $Q_0\cup Q_1=Q$, and we have $P_0\cap P_1=p$ and $Q_0\cap Q_1=q$. Without loss of generality, for node $u\in P_0$ and $n\notin p$, let's assume its next node $u^\prime\in Q_0$. If there exists another node $v\in P_0$, $v\notin p$ and $v^\prime\in Q_1$, then $vv\prime$ indicates that the particle at $v$ crosses a cutting edge, which contradicts the fact that no edge could cross cutting edges in laminar flow. Thus, any node in $P_0$ is mapped to a node in $Q_0$; and similarly for $P_1$ and $Q_1$. By definition, both $P_0Q_0$ and $P_1Q_1$ are 2.5D independent partitions.  % The proof still holds if we switch the indices of $Q_0$ and $Q_1$.
% \end{proof}
% 
% With the intuition previously explained for a cutting edge, a closed cutting path further help segment the spacetime into smaller consistent regions.  Basically, an edge formed by a node to its next node guarantees not to cross a closed cutting path.  
% Note that a closed cutting path may be degenerate, for example, a loop with

% Note that a non-closed cutting path will not lead to a subdivision, e.g., if the path terminates at an internal node.  

% \textbf{Domain decomposition with connected component labeling}.  


The assumption and lemma make it possible to decompose the domain via connected component labeling.  Specifically, cutting graphs break the lower and upper meshes into nonoverlapping 2.5D independent partitions for further processing. We use a union-find implementation to join edge-sharing triangles in each component and then match the components between the lower and upper meshes.  Special treatment may be needed for degenerate cases on domain boundaries; for example, if the cutting path in the upper mesh is on the boundary, a connected component in the lower mesh may not find a matching component in the upper mesh. In this case the resulting partition is degenerate, and we assign a path as the partition's upper component, as illustrated in~\Cref{fig:divide}.  


% With the cutting graph, one can subdivide the domain into 2.5D independent partitions

% with connected component labeling.  
%  While it is possible to iteratively bisect the entire domain by following the lemma, we reformulate the problem as a CCL problem so that one can find all partitions in one pass with a union-find or breadth-first search (BFS) implementation.  In this formulation, each triangle with be assigned with a unique partition.  The connectivity is defined as follows: two edge-sharing triangles belong to the same connected component (partition) if the shared edge is not a cutting path.  

% \HG{To be continued.}



% The CCL formulation will lead to the equivalent 








% \begin{restatable}{lem}{mylemma}\label{mylemma}
% Partition $PQ$ can subdivide into multiple nonoverlapping partitions, each as a connected components of triangles


% by a closed separating path.
% \end{restatable}






% \HG{I like explanations of rationales like below.}

% We leave the proof in Appendix, but the rationales are below.  
Decomposing the domain by cutting graphs provides several benefits.  % we can earn from cutting faces defined above.
First, a closed cutting path helps isolate trivial scenarios such as deformed prisms  (i.e., cells with triangle-to-triangle scenario in~\Cref{sec:spacetime_formulation}) from  more complex scenarios, because all quadrilaterals in deformed prisms are 4-node cycles formed by two cutting edges. 
Second, the triangulation of a 4-node cycle of a cutting edge can be directly given by staircase triangulation as introduced in~\Cref{sec:staircase}, which makes it easy to coordinate cross sections. 
Third, each partition resulting from our decomposition could be processed independently, according to Lemma 1.  % without considering other 2.5D independent partitions.
% The partition scheme conforms the underlying physics, i.e., transport behaviors, of the deforming spacetime. For example, i
% For example, in the case that deforming spacetime is driven by a laminar flow, no particles cross particles ahead of them in previous timestep, and thus no connections between nodes and their next nodes cross any cutting edges. No crossing connections guarantees tetrahedralizations in any two partitions do not interfere each other.
% That is why these systems are described by ``independent''.
% Also, separating these cutting faces are unavoidable in tetrahedralization, because only cycles with 3 mesh nodes are allowed in a valid tetrahedralization but these 4-node cycles exists without any operation in spacetime.





% There are several benefits we can earn from cutting faces defined above.
% \HG{What do you mean by hyperplanes?} 

% That is also why they are called ``real''.



% Formally, with a partition $PQ$, we define a separating edge as a 4-node cycle formed by four edges: a lower-mesh edge $uv\in P$, an upper-mesh edge $u^\prime v^\prime\in Q$, $uu^\prime$, and $vv^\prime$, where $u^\prime$ and $v^\prime$ is the next node of $u$ and $v$, respectively (see an example in \Cref{fig:divide}~{(a)}).  A separating edge also refers to edges $uv$ or $u^\prime v^\prime$ if both edges exist.  


% we define a separating edges as a 4-node cycle formed by four edges: an edge in lower mesh, an edge in upper mesh, and two edges between two nodes and their next nodes (see an example in \Cref{fig:divide}~{(a)}).  A separating path is a connected component of separating edge that share nodes.  In  

% \textbf{Segmenting spacetime by separating edges}.  We divide the whole spacetime by a set of \textit{cutting faces} first. % (or cutting faces for short) first.
% % We first divide it into pieces with smaller number of cells by several ``cutting faces.''
% A \textit{cutting face} is a 4-node cycle formed by four edges: an edge in lower mesh, an edge in upper mesh, and two edges between two nodes and their next nodes (see an example in \Cref{fig:divide}~{(a)}). Cutting faces can only be formed by two adjacent nodes in lower mesh mapped to another two adjacent nodes in upper mesh.
% % \HG{(as defined in 4.1)}.
% Every partition isolated by cutting faces might have triangles from only lower mesh, triangles from only upper mesh, or triangles from both meshes. See \Cref{fig:divide} as a mini example.
% % We divide the whole spacetime using all cutting faces (see \Cref{fig:divide}~{(b)} as a mini example).
% % A 3D piece isolated by cutting faces is called an \textit{independent system}, which notes that simplicial subdivision in this piece is unrelated to other pieces if the way to triangulate all cutting faces surrounding it is fixed. The lower and upper cells included in an independent system form \textit{lower patch} and \textit{upper patch} of it, respectively.



% \HG{As a rule of thumb, avoid using pronouns in all occasions unless there is an absolute need.  Pronouns take a reviewer extra seconds (or minutes) to find out what a pronoun refers to.  Even worse, there may be multiple interpretations of what a pronoun refers to, leading to a confusion.  You definitely want to do your reviewers a favor to save time and confusions.}

% The partitions derived from spacetime of XGC data in this step all have triangle strip and fan in lower mesh, \HG{No explaination?} i.e., all the nodes in lower mesh of a partition are on the boundary of it (\Cref{fig:indSys}~{(a)}).

% Some 2.5D independent partitions have \textit{inner node(s)} in upper mesh (\Cref{fig:indSys}~{(b)}), which corresponds to scenarios with some edge(s) mapped to path(s) in \Cref{fig:seven_scenarios}. These systems are further divided in next section.

% \HG{I remember there exists a situation when we have conflicting subdivision scheme of a real face, leading to a tet.  Do you plan to talk about the situation?}

\begin{figure}[!t]
    \centering
    \includegraphics[width=\columnwidth]{figs/divide.pdf}
    \caption{Mini-example showing process of segmenting spacetime by cutting paths in the divide stage. (a) 4-node cycle formed by two cutting edges highlighted in red in the spacetime. The cycle may or may not be a plane. (b) All the cutting paths in the spacetime. (c) Segmenting lower and upper meshes into spatial partitions by cutting paths. (d) Matching lower and upper spatial partitions by node-to-next-node correspondences. Only correspondences between the blue spatial partitions are shown. Some of them are matched with degenerate spatial partitions, that is, edges, such as the gray partition in the lower mesh and the pink partition in the upper mesh. (e) After matching the spatial partitions, we add node-to-next-node correspondences back and get 2.5D independent partitions.}
    \label{fig:divide}
\end{figure}

\begin{figure}[!t]
    \centering
    \includegraphics[width=\columnwidth]{figs/indSys.pdf}
    \caption{(a) and (b) show two 2.5D independent partitions derived from XGC data. The indices of mesh nodes in all figures agree with those in XGC data. (a) has no inner nodes in upper mesh, while (b) has one inner node 50153 on upper mesh. (c) shows two partitions of (b) derived by optional refinement in \textit{divide} step.
    % \HG{Now that the examples are from XGC, can we put the actual node IDs instead of $n_i$?}
    }
    \label{fig:indSys}
\end{figure}


% \subsubsection{Splitting by inner nodes}

\subsubsection{(Optionally) Refining a Partition by Inner Nodes}

Some 2.5D independent partitions isolated by cutting graphs could still include a large number of triangles or even inner nodes in the lower or upper mesh. To derive smaller partitions and further reduce complexity in the latter processes, we decompose every 2.5D independent partition with inner node(s).
% \textbf{(Optionally) refining a partition with inner nodes}.  
% There is an opportunity to further decompose a 2.5D independent partition with an inner node(s) into smaller partitions to further reduce complexity in the latter processes.
% We refine some 3D partitions that have inner nodes in upper mesh in this stage to make partitions smaller and further reduce complexity in latter processes.
% 
% , then there must be some inner node(s) in its upper mesh having no previous nodes, because \textit{boundary nodes} are in cutting faces and have previous nodes.
These partitions are split by \textit{penta-faces}
% \HG{Why/does it have to be a penta-face?  I remember you mentioned even larger faces.}
that include an edge in the lower mesh, two edges between two nodes and their next nodes, and a path with two edges in the upper mesh passing the inner node (e.g., $50155-50156-50154-50153-50152$ in \Cref{fig:indSys}~{(b)}). \Cref{fig:indSys}~{(c)} shows the split result of the partition in \Cref{fig:indSys}~{(b)}. There might be more than one penta-face crossing one inner node; we use the one that can balance the number of cells in both the lower and upper meshes in two split partitions.

Note that one can further refine a 2.5D independent partition with multiple adjacent inner nodes by multinode cycles, not only by penta-faces with 5-node cycles. Splitting by multinode cycles yields smaller partitions but makes triangulation and compatibility on cross sections more complex.

% \HG{What about the backfilling process?  With the two rounds mentioned above, we could get holes from the upper plane.}

% \HG{Another comment is about the (sub)systems/partitions that are heavily referred to in the rest of the paper.  A figure explaining the hierarchy may be helpful: spacetime --> partition --> (optional) subpartition}


\subsection{Lateral Triangulation} % (Bridging ``Divide'' and ``Conquer'' Stages)}
\label{sec:ftp}

% \subsection{Conquer: Piece to Tetrahedra\HG{Tetrahedralization within Spacetime Partitions}}

\begin{figure}[!t]
    \centering
    \includegraphics[width=0.35\textwidth]{figs/penta-face.pdf}
    \caption{All possible lateral triangulations for a penta-face. Each of them can be used as long as the subdivisions of two partitions containing the same penta-face agree on this common penta-face. Cutting by faces with more nodes in cycle is feasible but incurs more complexity of subdividing and coordinating cross sections. For a convex polygon with $n_p$ nodes, there are $\frac{1}{n_p-1}{2n_p-4 \choose n_p-2}$  triangulations for it~\cite{de2010triangulations}.}
    \label{fig:penta-face}
\end{figure}

\begin{figure}[!t]
    \centering
    \includegraphics[width=0.3\textwidth]{figs/fake_ill_case.pdf}
    \caption{Example of a quasi-ill-posed prism. (a) Staircase triangulation makes the deformed prism nontriangulable without Steiner point $s$, because  newly added line segment $a_0b_1$ intersects with newly created triangle $a_1a_2b_2$ at point $s$. Visible faces of the polyhedron include $a_1b_1s$, $b_1b_2s$, $a_1a_2b_2s$, and $a_0a_2b_2$. (b) The deformed prism is successfully triangulated after replacing $a_0b_1$ by $a_1b_0$ to separate 4-node cycle $a_0a_1b_1b_0$.}
    \label{fig:fake_ill}
\end{figure}


% \HG{To be continued.}

% When implementing the division of spacetime by cutting faces... It is noticeable that cutting faces may not be planes spanned by two linearly independent vectors, so separating the original 3D spacetime by cutting faces is hard to implement directly in 3D sense. Instead, w
% We first segment upper and lower mesh by upper and lower edges of all cutting faces, respectively (\Cref{fig:divide}~{(b)}). Segmenting a 2D mesh is achieved by finding connected components in the dual graph of the 2D mesh where links across cutting faces are broken. Resulting segmented pieces of meshes are called \textit{spatial partitions}. Then we match every spatial partition in lower mesh with a spatial partition in upper mesh by correspondences between nodes and their next nodes (\Cref{fig:divide}~{(c)}). Matching spatial partitions in reverse yields the results. Two matched spatial partitions together with node connectivity between them form a \textit{2.5D independent partition}. Some spatial partitions in one mesh fail to match with any spatial partitions and are actually mapped to a path in the other mesh (\Cref{fig:divide}~{(d)}), which form \textit{degenerated 2.5D independent partitions} with only triangles from one mesh.


% The goal of the conquer stage is to triangulate each 2.5D independent partition in an independent manner as much as possible.  We follow two steps for each 2.5D independent partition: (1) determine the triangulation scheme of lateral faces, which lead to a face-triangulated spacetime partition; (2) search for a possible 3D triangulation of the spacetime partition.  Note that the first step may require coordination with neighboring partitions if the default staircase rule failed to lead to a valid triangulation.  The second step, which is independent to neighboring partitions, is based on decision-tree search.

% , and the output is a tetrahedralization of this polyhedron. 
% We first triangulate lateral faces of the polyhedron by refining staircase triangulation in~\Cref{sec:staircase} to avoid nontriangulable deformed prisms. Then we search for possible triangulation of the polytope with triangulated faces by a decision-tree paradigm.

% \subsubsection{Face triangulation}
% \subsubsection{Triangulating lateral faces}
\label{sec:face_triangulation}
% \subsection{Simplicial Subdivision for Prismatic Meshes}

Before triangulating a 2.5D independent partition partition, one must first triangulate its lateral faces to form a polyhedron, namely, a face-triangulated partition.  Note that components from the lower/upper mesh are already simplicial; thus, we  need to consider only lateral faces.  The lateral faces created in the divide stage include 4-node and 5-node cycles.  There are two subdivisions for 4-node cycles and three subdivisions for 5-node cycles (shown in~\Cref{fig:penta-face}).  The subdivision of a lateral face can be arbitrary if the subdivisions of two adjacent partitions agree on it.  As stated in~\Cref{sec:staircase}, we follow the staircase triangulation rule---choosing the splitting edge with ascending node indices---to coordinate the separation of cross sections caused by cutting edges in the divide stage.  However, triangulation in deforming spacetime has more problems and possibly needs more than one iteration.




% It is nontrivial to 

% partition 3D spacetime by cutting faces because (1) a cutting face is usually nonplanar, (2) a cutting face has multiple possible triangulations, and (3) different cutting face triangulation will lead to different partitions, which in extreme cases, may result in a nontriangulable partition.


% A 2.5D independent partition with triangulated lower, upper, and lateral faces is called \textit{face-triangulated spacetime partition}.  

% \textbf{Ill-posed prism and quasi-ill-posed prism}. 
\textbf{Lateral triangulation of deformed prisms}. % Ill-posed prism and quasi-ill-posed prism}. 
As we discussed in~\Cref{sec:preliminaries}, not all deformed prisms can be triangulated, which we refer to as \textit{ill-posed prisms} in this paper.  Take a step back to consider a regular prism with eight potential ways to triangulate its lateral quadrilaterals, as illustrated in~\Cref{sec:staircase}; six triangulation schemes can successfully split the prism into three tetrahedra, except for the two schemes that add only opposite diagonals.  Note that staircase triangulation is a particular case of the six successful schemes. Now consider a deformed prism; not all of the six separations can triangulate a deformed prism.  Some deformed prisms (called \textit{quasi-ill-posed prisms}) 
fail to triangulate by the staircase triangulation but are successfully triangulated by another triangulation scheme (see \Cref{fig:fake_ill} as an example). Unfortunately, ill-posed prisms cannot be triangulated by any scheme, as shown by the example in~\Cref{tab:ill_cases}.

% \textbf{Lateral face triangulation of deformed prisms}. Unlike a regular prism, not all deformed prisms can be triangulated, which are referred to as \emph{ill-prisms} in the rest of this paper.  
% % could be nontriangulable. 
% Every prism, as well as deformed prism, has eight separations of the three lateral quadrilaterals, as stated in~\Cref{sec:staircase}. Except the two separations that add only opposite diagonals, other six separations (as listed in~\Cref{tab:ill_cases}) can successfully divide a prism into three tetrahedra, and staircase triangulation provides one possible separation. However, not all of the six separations can tetrahedralize a deformed prism. Some deformed prisms (called ``fake'' nontriangulable prism) fail to be triangulated by the separation given by staircase triangulation but are successfully triangulated by another separation (see \Cref{fig:fake_ill} as an example). Some deformed prisms (called ``real'' nontriangulable prism) even cannot be triangulated by any separations, as shown by the example in~\Cref{tab:ill_cases}.

\begin{table}[!t]
  \caption{%
  	Check list for valid subdivision of deformed prisms. If a triangulation is valid, then the two corresponding conditions regarding whether two points are on different sides of a plane shown in the last three columns should  both be satisfied. If one of the points is in the plane, then there is a tetrahedron degenerating into a quadrilateral. The example shown in the third column is an ill-posed prism, which is created by two equilateral triangles where one can be aligned with the other by rotating it by $\pi/2$ around its centroid. The ill-posed prism cannot be triangulated by any of the six triangulations.
  }
  \label{tab:ill_cases}
  % \scriptsize%
  \centering%
  \begin{tabular}{m{0.1in}|m{0.45in}|m{0.8in}|m{0.2in}|m{0.2in}|m{0.6in}}
  	\toprule
  	\multirow{2}*{ID} & \multirow{2}*{\tabincell{Face\\separations}} & \multirow{2}*{Subdivision} & \multicolumn{3}{c}{Conditions to be checked} \\ \cline{4-6}
   & & & pt 1 & pt 2 & plane \\
   \midrule
   \end{tabular}
  \renewcommand{\arraystretch}{2.3}
  \begin{tabular}{m{0.1in}|m{0.45in}|m{0.8in}|m{0.2in}|m{0.2in}|m{0.6in}}
  	\multirow{3}*{1} & $a_0-b_1$ & \multirow{3}*{\includegraphics[width=0.8in]{figs/ill_cases/ill_case1.pdf}} & $a_2$ & $b_1$ & $a_0-a_1-b_2$ \\
   & $a_1-b_2$ & & & & \\
   & $a_0-b_2$ & & $a_1$ & $b_0$ & $a_0-b_1-b_2$ \\ \hline
  	\multirow{3}*{2} & $a_0-b_1$ & \multirow{3}*{\includegraphics[width=0.8in]{figs/ill_cases/ill_case2.pdf}} & $a_1$ & $b_2$ & $a_0-a_2-b_1$ \\
   & $a_2-b_1$ & & & & \\
   & $a_0-b_2$ & & $a_2$ & $b_0$ & $a_0-b_1-b_2$ \\ \hline
  	\multirow{3}*{3} & $a_1-b_0$ & \multirow{3}*{\includegraphics[width=0.8in]{figs/ill_cases/ill_case3.pdf}} & $a_2$ & $b_0$ & $a_0-a_1-b_2$ \\
   & $a_1-b_2$ & & & & \\
   & $a_0-b_2$ & & $a_0$ & $b_1$ & $a_1-b_0-b_2$ \\ \hline
  	\multirow{3}*{4} & $a_0-b_1$ & \multirow{3}*{\includegraphics[width=0.8in]{figs/ill_cases/ill_case4.pdf}} & $a_1$ & $b_0$ & $a_0-a_2-b_1$ \\
   & $a_2-b_1$ & & & & \\
   & $a_2-b_0$ & & $a_0$ & $b_2$ & $a_2-b_0-b_1$ \\ \hline
  	\multirow{3}*{5} & $a_1-b_0$ & \multirow{3}*{\includegraphics[width=0.8in]{figs/ill_cases/ill_case5.pdf}} & $a_0$ & $b_1$ & $a_1-a_2-b_0$ \\
   & $a_2-b_1$ & & & & \\
   & $a_2-b_0$ & & $a_1$ & $b_2$ & $a_2-b_0-b_1$ \\ \hline
  	\multirow{3}*{6} & $a_1-b_0$ & \multirow{3}*{\includegraphics[width=0.8in]{figs/ill_cases/ill_case6.pdf}} & $a_0$ & $b_2$ & $a_1-a_2-b_0$ \\
   & $a_1-b_2$ & & & & \\
   & $a_2-b_0$ & & $a_2$ & $b_1$ & $a_1-b_0-b_2$ \\
  	\bottomrule
  \end{tabular}%
\end{table}

% To triangulate faces of a deformed prism, we go following steps. First, we check if it can be triangulated by staircase triangulation. Second, if not, then the deformed prism is at least a ``fake'' nontriangulable prism. We check all separations in~\Cref{tab:ill_cases} to see if one of the separations can triangulate the deform prism. If two line-plane intersections in the last two columns in~\Cref{tab:ill_cases} for a separation all turn out to be false, then the triangulation is valid. Third, if all separations cannot triangulate the deform prism, then it is a ``real'' nontriangulable prism. It also means that the temporal resolution of data should be increased to reduce the rotation between two successive timesteps. In XGC, the toroidal coordinates is originally evenly separated into 16 pieces. We upsample it to 48 pieces to avoid ``real'' nontriangulable prism.

We follow three steps to triangulate the lateral faces of a deformed prism.   
First, we check whether staircase triangulation could lead to a valid triangulation.  
If not, second, the deformed prism is at least a quasi-ill-posed prism, which requires checking all separation schemes in~\Cref{tab:ill_cases} to see whether one of the separations can triangulate the deformed prism.  There are two conditions to check for whether two points are on different sides of a plane defined by a triangle, as shown in the last three columns in~\Cref{tab:ill_cases}, to see whether a triangulation is valid. Third, if all separations cannot triangulate the deformed prism, it is an ill-posed prism. It also means that the temporal resolution of data should be increased to reduce the deformation between two successive timesteps, as we previously discussed in~\Cref{sec:spacetime_formulation}.
% In XGC, the toroidal coordinates is originally evenly separated into 16 pieces. We upsample it to 48 pieces to avoid ill-posed prism to meet the assumption in~\Cref{sec:assumption}.

\textbf{Lateral triangulation of complex partitions beyond prisms}.
% \textbf{Lateral face triangulation of more complex 2.5D independent partitions}. 
We choose an arbitrary sequence of separations of the 2.5D independent partition's lateral faces that are not constrained by triangulation of quasi-ill-posed prisms, and we try to triangulate it. If we fails to triangulate the 2.5D independent partition with given triangulated faces, then we change the separation of one lateral face and try to triangulate the newly created polyhedron, as shown in~\Cref{fig:pipeline}. Different from deformed prisms, it is hard to determine whether a polyhedron is triangulable without trying to triangulate it first.

% \HG{Looks like you directly jumped into technical details without explaining the rationale behind.  Can you outline a few reasons why we favor staircase triangulation for trivial cases (assuming trivial cases are triangle-to-\{triangle|edge|node\}?  You may also want to state that trivial cases are handled in preprocessing; only nontrivial cases are handled by the search tree algorithm.}
% Staircase triangulation can directly solve three scenarios in previous section, \textit{triangle-to-triangle}, \textit{triangle-to-edge}, and \textit{triangle-to-node}. Every \textit{triangle-to-triangle} pair forms a \textit{prism}, whose simplicial subdivision can be solved by following \textit{staircase triangulation}~\cite{de2010triangulations}. We review terms and staircase triangulation in this section.
% \HG{What about Sch\"onhardt prisms?  If we assumed no Sch\"onhardt prisms, we need to justify why/how they can be avoided.
% \HG{Nevertheless, I would expect to see a list of objectives/constraints in a formal manner; similar to the list in the intro but more strict/mathematical.}

\begin{table}[!t]
  \caption{%
  	Types of neighbors around pivot nodes, their triangulation ways, and corresponding separated tetrahedra. Nodes with superscripts are on the upper mesh while those without superscripts are on the lower mesh. Red lines represent newly added edges in triangulation. ``Pviot'' in the last column is replaced by ``p'' for short.
   % \HG{The symbol $pivot$ is recurring in the last column; use $p$ instead?}
  }
  \label{tab:pivot_cases}
  \scriptsize%
  \centering%
  \begin{tabular}{m{0.1in}|m{0.2in}|m{0.8in}|m{0.85in}|m{0.5in}}
  	\toprule
  	Min deg & Pivot node & \tabincell{pivot node and\\its neighbors} & Volume triangulations & Separated tetrahedra \\
   \midrule
      % \includegraphics[width=0.3in]{figs/prism.pdf}} & \rotatebox{90}{Vis/SciVis} &   \rotatebox{90}{SciVis conf} &   \rotatebox{90}{InfoVis} &   \rotatebox{90}{VAST} &   \rotatebox{90}{VAST conf} &   \rotatebox{90}{TVCG @ VIS} &   \rotatebox{90}{CG\&A @ VIS} &   \rotatebox{90}{VIS/VisWeek} \rotatebox{90}{incl.\ TVCG/CG\&A}   &   \rotatebox{90}{VIS/VisWeek} \rotatebox{90}{w/o TVCG/CG\&A}   \\
  	\multirow{2}*{3} & \tabincell{in\\lower\\mesh} & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case1.pdf} & - & \{$p$, $n_0$, $n_1$, $n_x'$\} \\ \cline{2-5}
  	& \tabincell{in\\upper\\mesh} & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case2.pdf} & - & \{$p$, $n_0'$, $n_1'$, $n_x$\} \\ \hline
  	\multirow{10}*{4} & \multirow{5}*{\tabincell{in\\lower\\mesh}} & \multirow{2}*{\includegraphics[width=0.7in]{figs/pivot_cases/pivot_case3.pdf}} & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case3-1.pdf} & \tabincell{\{$p$, $n_0$,\\$n_x'$, $n_y'$\},\\\{$p$, $n_0$,\\$n_1$, $n_y'$\}}\\ \cline{4-5}
  	& & & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case3-2.pdf} & \tabincell{\{$p$, $n_0$,\\$n_1$, $n_x'$\},\\\{$p$, $n_1$,\\$n_x'$, $n_y'$\}} \\ \cline{3-5}
  	& & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case4.pdf} & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case4-1.pdf} & \tabincell{\{$p$, $n_0$,\\$n_x'$, $n_y'$\},\\\{$p$, $n_0$,\\$n_y'$, $n_z'$\}} \\ \cline{3-5}
  	& & \includegraphics[width=0.6in]{figs/pivot_cases/pivot_case5.pdf} & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case5-1.pdf} & \tabincell{\{$p$, $n_0$,\\$n_2$, $n_x'$\},\\\{$p$, $n_1$,\\$n_2$, $n_x'$\}} \\ \cline{3-5}
  	& & \includegraphics[width=0.6in]{figs/pivot_cases/pivot_case6.pdf} & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case6-1.pdf} & \tabincell{\{$p$, $n_0$,\\$n_2$, $n_x'$\},\\\{$p$, $n_1$,\\$n_2$, $n_x'$\}} \\ \cline{2-5}
  	& \multirow{5}*{\tabincell{in\\upper\\mesh}} & \multirow{2}*{\includegraphics[width=0.7in]{figs/pivot_cases/pivot_case7.pdf}} & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case7-1.pdf} & \tabincell{\{$p$, $n_x$,\\$n_0'$, $n_1'$\},\\\{$p$, $n_x$,\\$n_y$, $n_1'$\}}  \\ \cline{4-5}
  	& & & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case7-2.pdf} & \tabincell{\{$p$, $n_x$,\\$n_y$, $n_0'$\},\\\{$p$, $n_y$,\\$n_0'$, $n_1'$\}}  \\ \cline{3-5}
  	& & \includegraphics[width=0.6in]{figs/pivot_cases/pivot_case10.pdf} & \includegraphics[width=0.6in]{figs/pivot_cases/pivot_case10-1.pdf} & \tabincell{\{$p$, $n_x$,\\$n_y$, $n_0'$\},\\\{$p$, $n_y$,\\$n_z$, $n_0'$\}} \\ \cline{3-5}
  	& & \includegraphics[width=0.6in]{figs/pivot_cases/pivot_case8.pdf} & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case8-1.pdf} & \tabincell{\{$p$, $n_x$,\\$n_0'$, $n_2'$\},\\\{$p$, $n_x$,\\$n_1'$, $n_2'$\}} \\ \cline{3-5}
  	& & \includegraphics[width=0.6in]{figs/pivot_cases/pivot_case9.pdf} & \includegraphics[width=0.7in]{figs/pivot_cases/pivot_case9-1.pdf} & \tabincell{\{$p$, $n_x$,\\$n_0'$, $n_2'$\},\\\{$p$, $n_x$,\\$n_1'$, $n_2'$\}} \\
  	% \midrule
  	\bottomrule
  \end{tabular}%
\end{table}

% \HG{You may need to clearly state the inputs/outputs/constraints/assumptions for the conquer stage.  My understandings are as follows.  
% %
% Inputs: a (nonconvex) polytope (spacetime partition) with all faces triangulated (upper/lower triangulations mandated by the input 2D mesh, lateral faces by staircase triangulation).  
% %
% Outputs: a full tetrahedralization of the polytope.  
% %
% Constraints: (1) the partition boundary bounds all simplices, (2) no degenerate tets/simplices and thus no new edges in the upper/lower mesh, (3) for any lower-mesh node $n$, newly added edge must be connecting to a 1-hop neighbor of the next node $n^\prime$.  
% % 
% We did not mention constraint (3) anywhere, but it may be a good idea to mention it earlier.}


\subsection{Volume Triangulation} %  (the ``Conquer'' Stage)} % : Triangulating 2.5D independent partitions}
\label{sec:conquer}

% \subsubsection{Two-tier decision-tree for triangulation search} % triangulating 3D partitions}
% Node eliminate in a decision tree}
% \subsubsection{Node eliminate in a decision tree}

Once each 2.5D independent partition is face-triangulated into a polyhedron, we further triangulate each polyhedron into tetrahedra.  We propose a \textit{node elimination} algorithm that removes nodes one by one as well as edges incident on it as tetrahedra until all nodes are removed.  In searching for the sequence of nodes to be removed, there are multiple decisions to make, and we need to trace back to previous steps if we meet with nontriangulable remaining polyhedron. Thus, we formulate the node elimination algorithm into a decision-tree paradigm.

\textbf{Two-tiered decision tree for volume triangulation}. The decision tree makes two types of decisions alternately. In levels with odd numbers (odd levels), the tree makes a decision on which node to be eliminated (so-called \textit{pivot node}); in levels with even numbers (even levels), the tree chooses one way to triangulate the polyhedron formed by pivot node and edges incident on the node (so-called \textit{isolated polyhedron}).
% \HG{You didn't mention odd/even levels before}. 
% \HG{the term ``volume'' may be ambiguous; residue polytope?},
Every tree node contains information about separated tetrahedra, remaining polyhedron, pivot node, and unvisited children of the tree node.
% This algorithm can tetrahedralize any 3D volume whose surface (including base faces and lateral faces) has already been triangulated.

% \HG{I feel directly jumping into the search tree may be difficult for readers.  I recalled the discussion with Zeyu about the SAT solver.  Although we are not using SAT, for readers' sake, it may be better to define the problem as finding a node-elimination sequence and then define the decision/search tree.  Note that the search tree is a tool, not the goal.}

\textbf{Odd levels: decision on choosing a pivot node}. The pivot node can be any node in the polyhedron, and one can introduce any heuristic to choose pivot node. Here, we prioritize all mesh nodes with minimum degree
% \HG{Why?}
in the current remaining polyhedron as pivot nodes and create a tree node for every pivot node as child for current tree node. One advantage of this heuristic is that nodes with minimum degree are incident to fewer edges and thus incur simpler isolated polyhedra. See~\Cref{alg:pivot_choose} for pseudo-code of this step.
% \HG{In theory, any node could be chosen as a pivot; min-deg is just a heuristics but requires explanation.}


% \HG{I think this paragraph may need some work.  I had trouble understanding the details, which means a reviewer will surely be confused.  Clarify what are the goals and constraints, and explain the purpose/use of the look-up table.}
\textbf{Even levels: decision on triangulating polyhedron around the pivot node}.
% \HG{I think the decision should be the triangulation about the pivot; we already decided to remove the pivot node at this step.}
We continue with a pivot node chosen in the last level and eliminate the pivot. We first isolate the pivot, its neighbors, and edges between the pivot and its neighbors.
% \HG{``these'' is bad here; spell out?}.
All the isolated elements form an isolated polyhedron that is a subset of the inputted polyhedron.
% This kind \HG{I had trouble figuring out what is ``this kind.''} of piece
The isolated polyhedron usually contains 4--5 nodes; all possible isolated polyhedra with no more than 5 nodes are shown in the third column of \Cref{tab:pivot_cases}. Then we triangulate the isolated polyhedron whose tetrahedralization is shown in~\Cref{tab:pivot_cases}. It is easy to exhaust all possibilities for an isolated polyhedron,
%GAIL - d you mean easy to exhaust all possibilities for the polyhedra
remove derived tetrahedra from the current polyhedron, and update the remaining polyhedron. See~\Cref{alg:pivot_remove} for pseudo-code of this step.
%GAIL - do you mean update the remaining polyhedra
Note that not all the tetrahedralization ways in~\Cref{tab:pivot_cases} are feasible for a polyhedron; two conditions need checking: (1) if every newly created link is inside of the polyhedron (constraints introduced by nonconvexity) and (2) if any of separated tetrahedra has four coplanar mesh nodes. We create a tree node for every feasible tetrahedralization as a child for the tree node representing the pivot node in the upper level, and we update separated tetrahedra and remaining polyhedron 
%GAIL polyhedron or polyhedra in the preceding and following
for these children. In the next iteration, we repeat operations in the two levels for remaining polyhedron while the remaining polyhedron is non-empty. \Cref{fig:node-elimination} shows the node elimination algorithm in the structure of the decision tree, taking an arbitrary polyhedron with triangulated lateral faces as an example. See~\Cref{alg:pivot_remove} for pseudo-code of the whole-node elimination algorithm.

\begin{figure}[!ht]
    \centering
    \includegraphics[width=0.458\textwidth]{figs/node_elimination.pdf}
    \caption{Node elimination algorithm shown in a decision-tree paradigm. In odd levels, it makes a decision on ``which mesh node to eliminate'' among mesh nodes with minimum degree; in even levels, it makes a decision on ``how to eliminate'' based on the isolated polyhedron formed by the pivot node and its neighbors.
    % There are ten nodes in original volume, and thus this tree has $2(n-3)+1=15$ levels. The algorithm starts with a volume and an empty separated tetrahedra set, and terminates with an empty volume with all separated tetrahedra.
    }
    \label{fig:node-elimination}
\end{figure}

\textbf{Complexity of the algorithm}. The best and most common case is to terminate with a single branch ($O(n)$); the worst and nearly impossible case is to fully expand the tree ($O(n2^n)$). This tree has $2(n_v-3)+1$ levels ($n_v$ stands for the number of nodes in the face-triangulated polyhedron to be triangulated)
% \HG{you already used the symbol $n$ as a node; change to another symbol?})
because nearly every node (except the last three nodes that can be directly removed with the last tetrahedron) will be a pivot exactly once in any branch that connecting root and a leaf, and every pivot generates two levels. To find a tetrahedralization way, there is no need to traverse all possible tree node. The algorithm terminates when it finds a tree node whose remaining polyhedron is empty. When it meets with an indivisible isolated polyhedron around a pivot node (e.g., \Cref{fig:staircase}~{(b)}), it traces back and continues with the next unvisited child of the current tree node's parent.
% \HG{what is $e$?}. \HG{Move the conclusion to the top of the paragraph.}

\textbf{Completeness of the algorithm}. The algorithm is guaranteed to find a volume triangulation that allows only connections between nodes and their 2-hop neighbors if there is one, because this decision tree covers all possible connections within 2-hop neighbors. If there are no unvisited nodes in the tree, the algorithm terminates with no solution, which means that the inputted polyhedron is indivisible if only links between 2-hop neighbors are allowed to create. It is reasonable to consider links between 2-hop neighbors because links between greater than n-hop neighbors degrade accuracy of interpolation in tetrahedra.
% One solution is to find out whether other face triangulations can make the polyhedron triangulable.

\textbf{Optimality of the algorithm.} If cost of a solution in the algorithm is defined to be the number of visited nodes when finding a feasible volume triangulation, then our heuristics for choosing pivot may not direct the algorithm to an optimal solution. Choosing mesh nodes with minimum degree as a pivot may not be the most complexity-saving heuristics: removing a mesh node with minimum degree might introduce nonconvexity to the remaining polyhedron, while removing mesh nodes with higher degrees could maintain convexity. Nonconvexity brings more constraints on feasible connections that one can create and raises the risk of making the remaining polyhedron indivisible, which requires tracing-back in the decision tree and thus increases cost.

% \HG{Refer to the algorithms earlier.  Even if few readers read the pseudocode, no one will ever notice the pseudocode after you described everything.}
% The full algorithm is shown in~\Cref{alg:node_elimination}. Two functions, one for making decision on choosing pivot in odd levels (\Cref{alg:pivot_choose}) and the other for making decision on removing pivot in even levels (\Cref{alg:pivot_remove}), are recursively called in it.

\begin{algorithm}[!ht]
\caption{Node Elimination Algorithm}\label{alg:node_elimination}
\KwData{a face-triangulated polyhedron}
\KwResult{tetrahedralization result for the polyhedron}
root $\gets$ Node(\textit{separated\_tets} = $\emptyset$, \textit{remaining\_volume} = volume,
 \textit{pivot} = \texttt{None}, \textit{unvisited\_children} = $\emptyset$)\;
tree $\gets$ Tree(root=root)\;
current\_node $\gets$\ root\;
divisible $\gets$ \texttt{True}\;
\While{current\_node.remaining\_volume$\neq\emptyset$}{
\If{divisible == \texttt{False}}{
    current\_node $\gets$ the next unvisited child of current\_node's parent\;
    \If{no unvisited nodes in the tree}{ \textit{Throw}(``This polyhedron is indivisible.'')
    % \Comment{this volume is indivisible}
}
\eIf{current\_node.current\_node.pivot==\texttt{None}}{
    current\_node, tree, divisible $\gets$ \textit{PivotChoose}(current\_node, tree)\;
    }{
    current\_node, tree, divisible $\gets$ \textit{PivotRemove}(current\_node, tree)\;
    }
}}
\textbf{return} current\_node.\textit{separated\_tets}
\end{algorithm}

\begin{algorithm}[!ht]
\caption{PivotChoose}\label{alg:pivot_choose}
\KwData{tree\_node, tree}
find mesh nodes with minimum degree\;
generate children representing pivot nodes for tree\_node\;
divisible $\gets$ \texttt{False} \;
\While{divisible==\texttt{False}}{
\eIf{tree\_node.unvisited\_children$\neq\emptyset$}{
current\_node $\gets$ tree\_node.\textit{unvisited\_children}[0]\;
remove current\_node from tree\_node.\textit{unvisited\_children}\;
current\_node, tree, divisible $\gets$ \textit{PivotRemove}(current\_node, tree)\;
}{
\textbf{return} current\_node, tree, \texttt{False}\;
}
\textbf{return} current\_node, tree, \texttt{True}\;
}
\end{algorithm}

\begin{algorithm}[!ht]
\caption{PivotRemove}\label{alg:pivot_remove}
\KwData{tree\_node, tree}
check the cycles formed by pivot and its neighbors; determine its type\;
divisible $\gets$ if at least one divisible way exists\;
\eIf{divisible}{
generate children representing dividing ways for tree\_node\;
current\_node $\gets$ the first unvisited child of tree\_node\;
remove current\_node from tree\_node.\textit{unvisited\_children}\;
update current\_node.\textit{separated\_tets} and current\_node.\textit{remaining\_volume}\;
return current\_node, tree, \texttt{True}
}{
return current\_node, tree, \texttt{False}
}

\end{algorithm}

\begin{figure}[!t]
    \centering
    \includegraphics[width=\columnwidth]{figs/XGC_geometry.pdf}
    \caption{
    % \HG{I like this figure.  Is it possible to make the gray curves look more twisted like magnetic lines?}
    Illustration
    % of geometric modeling \HG{This is not about geometric modeling.  Saying the cylindrical coordinate system in tokamak is fine.}
    of a tokamak with a cylindrical coordinate system $(R,Z,\phi)^\intercal$, where $R$, $Z$, and $\phi$ are radial, axial, and toroidal coordinates, respectively. Two 2D poloidal planes with $\Delta\phi=\pi$ are shown, but actually we use XGC data discretized into 16 poloidal planes. The gray lines show trajectories of arbitrarily chosen particles. Each poloidal plane is meshed into an identical unstructured triangular grid. Part of the triangular mesh in the circled region is shown on the right. The whole mesh has 56,980 nodes and 112,655 cells.}
    \label{fig:xgc_geometry}
\end{figure}


% \section{Case Studies with XGC Fusion Plasma Simulations}
\section{Case Studies: Deforming Spacetime Induced by a Magnetic Field in Fusion Plasma} % s and Problems}
\label{sec:case-studies}

A primary application driver for our research is tokamak fusion; beyond fusion, the spacetime meshing capability can also benefit the visualization and analysis of various science applications that involve a deformed spacetime.
% 
% \subsection{Deforming spacetime induced by a magnetic field in fusion plasma}
% 
In fusion energy science, tokamaks---machines in a donut shape that confine burning plasma with strong magnetic fields---are the leading concept for future fusion power generation.  However, tokamaks are vulnerable to disruptions caused by plasma turbulence, and a critical mission in fusion research is to understand the plasma turbulence structures and how to reduce turbulence.  To this end, scientists developed large-scale simulations such as XGC (X-Point Included Gyrokinetic Code) to study tokamaks.  

Scientists favor 3D visualizations such as volume rendering, isosurfacing, and feature curve rendering of \emph{blobs}---filament structures of high turbulence in tokamaks that may cause disruptions and damage billion-dollar devices.  The specific scalar function for visualizing blobs is \emph{normalized electrostatic potential perturbation} ($\delta n_e / \delta n_{e0}$), denoted as $f$ in the rest of this paper. To date, visualization of blobs has mainly focused on 2D cross sections of blobs because of the prohibitive cost of field-following interpolation, as explained below. 

% A key challenge in visualizing and analyzing XGC data is how field variables (e.g., temperature, density, and magnetic potentials) are 

Scalar functions (e.g., temperature, density, and magnetic potentials) in XGC are interpolated via a \emph{particle shape function}, that is, field-following interpolation in Eq.~\eqref{eq:ff}. \Cref{fig:xgc_geometry} illustrates a cylindrical coordinate system (radial, axial, and toroidal coordinates $R$, $Z$, and $\phi$) representing the computational domain in XGC.  XGC uniformly subdivides the toroidal direction into a finite number of poloidal planes (i.e., $RZ$-planes).  Each poloidal plane uses the same triangular mesh.  For example, a typical XGC simulation discretizes the domain with 16 poloidal planes, each with a triangular mesh of $O(10^5)$ nodes and $O(10^5)$ triangles.


% have drawn much attention to the plasma physics community.  Although blobs are filament structures across multiple poloidal planes (cutting plane of the torus), up to date, blob extraction and tracking have been mainly focused on 2D structures.  For example, ...

In this case study we interpret XGC's spatial domain as a deformed spacetime induced by magnetic fields. \Cref{fig:xgc_interp} (a) illustrates the magnetic-following interpolation scheme.  The field variable $f$ is given at every  node of triangular meshes in all poloidal planes.  To interpolate $f$ at an arbitrary location $(R, Z, \phi)^\intercal$, one must first calculate a streamline of the magnetic field $\mathbb{B}$ (i.e., a magnetic line), in both directions, seeded from the given location.  The magnetic line normally intersects two poloidal planes at $\phi_i$ and $\phi_{i+1}$ with $\phi_i\leq\phi\leq\phi_{i+1}$, where $i$ and $i+1$ are the indices of poloidal planes in XGC's toroidal discretization.  Assuming the function value $f$ is $f_i$ and $f_{i+1}$ at the streamline's intersection on $i$th and $(i+1)$th poloidal plane, respectively, one can approximate $f$ as the linear combination of $f_i$ and $f_{i+1}$. 
% \begin{equation}
%   f(r, z, \phi) \approx
%   \frac{\phi_{i+1} - \phi}{\phi_{i+1}-\phi_i}f_i + 
%       \frac{\phi - \phi_i}{\phi_{i+1}-\phi_i}f_{i+1}.\label{eq:interpolation}
% \end{equation}
% 

The magnetic-following interpolation poses challenges in both forward and inverse evaluations. Here, we refer to \emph{forward evaluation} as the interpolation scheme and \emph{inverse evaluation} as  finding  locations that satisfy given criteria, for example, finding the zero-level-set where $f=0$. Both forward and inverse evaluations are expensive.  Forward evaluations---frequently used by volume rendering and particle tracing---require magnetic line tracing and are thus expensive.  Inverse evaluations such as isosurfacing and critical point extraction are complicated and may involve numerical optimization for root-finding.  Note that We regard the toroidal coordinates as  temporal coordinates in this work, which is different from real temporal dimensions in XGC data.

% Without loss of generality, we exemplify our algorithms with XGC meshes and data in this paper.  In later sections, we will demonstrate how we simplify forward and inverse evaluation with simplicial complexes representing the deformed spacetime induced by the magnetic-following interpolation scheme.

% \subsection{Deforming spacetime induced by a velocity field in rotary mechanics}
\subsection{Evaluation Schemes}

We demonstrate and evaluate the use of our methodology with XGC fusion plasma simulations, whose 2D triangular mesh has 56,980 nodes and 112,655 cells. Our algorithm yields 339,626 tetrahedra for its 3D spacetime with $ 56,980\times2=113,960$ nodes.  We define and use two types of approximations:
\begin{itemize}\setlength{\parskip}{0pt}
    \item Forward approximation: interpolation of function values at given 3D locations.  Examples include volume rendering and particle tracing, which require sampling at arbitrary locations for ray/path integration.
    \item Inverse approximation: root-finding of 3D locations with given function values.  Examples include isosurfacing in scalar-valued functions and critical-point-finding in vector-valued functions.
\end{itemize}
We compare three types of interpolation schemes:
\begin{itemize}\setlength{\parskip}{0pt}
    \item Magnetic-following (MF) interpolation, that is, the scheme defined in Eq.~\eqref{eq:ff} that requires magnetic-line tracing
    \item Our method: Piecewise linear (PL) interpolation induced by our simplicial meshing scheme that models deforming space
    \item A naive baseline: Straight linear (SL) interpolation induced by simplicial meshing on uniform spacetime, where every node's next node is itself.
    % \HG{I think it may be worth comparing with the straight linear (SL) interpolation, which is a wrong interpolation but worth to demonstrate how things look like if the deforming spacetime is not considered.}
\end{itemize}
The three interpolation schemes are measured by mean squared error (MSE), peak signal-to-noise ratio (PSNR), and running time. Compared with magnetic-following interpolation, the PL interpolation significantly reduces the computational cost of all visualization tasks, because no magnetic-line tracing is involved for the interpolation.  We recognize that the PL approximation introduces error, and thus we treat the magnetic-following interpolation as the ground truth and quantitatively evaluate the error resulting from our PL interpolation.  We also work together with domain scientists to evaluate our results qualitatively.


\subsection{Toroidal Upsampling (Forward Approximation)}
\label{sec:upsampling}

% is to upsample a scalar function (\emph{normalized perturbation of electric potential}, namely $dneOverne0$) in XGC output data, 

The objective of the toroidal upsamplling case study is to upsample XGC's 3D scalar field data, which are sampled with relatively low resolution (e.g., $n_\phi=16$) uniformly along the toroidal ($\phi$) direction at $\phi=2\pi i/16$ ($i=0,1,...,15$).  With each given $\phi$, science variables (e.g., density and temperature) are given at nodes on 2D triangular meshes.  We interpolate simulation outputs to obtain values of nodes at two arbitrary $\phi$'s, $(5/32)2\pi$ and $(27/32)2\pi$, which are both in the middle of two consecutive sampled toroidal angles, to retrieve a much higher toroidal resolution of the data by MF, PL, and SL interpolations. These three interpolations are implemented on a machine with an Apple M1 chip and 64 GB system memory with no parallel technology. Results of three kinds of interpolations and quantitative evaluation are shown in \Cref{fig:unsampling} and \Cref{tab:performance}, respectively. PL interpolation results show a similar pattern of values, low MSE and high PSNR, and significantly less running time ($O(10^4)$ faster) when compared with MF interpolation, as shown in~\Cref{fig:phi_vs_time}.
% (e.g., $n_\phi^\prime=1,024$ or $360^\circ/1,024=0.352^\circ$) of the data.
% \HG{Need accuracy and performance data.}

\begin{figure}
    \centering
    \includegraphics[width=0.8\columnwidth]{figs/unsampling.pdf}
    \caption{Upsampling results of MF, PL, and SL interpolation methods. The overall patterns shown by different interpolation methods are similar. However, PL interpolation has closer values to MF while SL has a smaller range of values compared with MF and PL. In several local regions (such as those marked by red boxes), SL even shows  shapes different from those of the other two.}
    \label{fig:unsampling}
\end{figure}

\begin{table}
    \centering
    \caption{Performance of MF, PL, and SL
    % \HG{What about ``straight'' linear interpolation?}
    interpolation methods on XGC data with 56,980 nodes and 112,655 cells. MF interpolation results are taken as ground truth. The two $\phi$'s are selected because they show the best and worst approximation quality in the middle of two sampled $\phi$'s.
    % \HG{Also, there is a gap of PSNR (46 vs 51 dB).  Can you explain?  I remember we discussed about sweeping different phi to demonstrate how the error would change.}
    }
    \begin{tabular}{c|c|c|c|c}
        \toprule
        $\phi$ & Interp & MSE & PSNR (dB) & Time (sec) \\ \midrule
        \multirow{3}*{$2\pi\times\frac{5}{32}$} & MF & - & - & 183,704 \\ \cline{2-5}
        & PL & 2.714e-5 & 48.569 & 13.280 \\ \cline{2-5}
        & SL & 2.632e-4 & 38.704 & 10.398 \\ \hline
        \multirow{3}*{$2\pi\times\frac{27}{32}$} & MF & - & - & 183,703 \\ \cline{2-5}
        & PL & 2.544e-5 & 38.304 & 13.283 \\\cline{2-5}
        & SL & 2.601e-4 & 28.208 & 10.398 \\ \bottomrule
    \end{tabular}
    \label{tab:performance}
\end{table}

We also uniformly took 1,024 poloidal planes over $\phi\in[0,2\pi)$ to show the tendency of PSNR with varying $\phi$ (\Cref{fig:phi_vs_psnr}), taking MF interpolation results as ground truth. Cells in SL interpolation do not conform with the variation in time-dependent data and thus interpolate the value of a position with four irrelevant mesh nodes. Therefore, PL is expected to have better performance than SL. As shown in \Cref{fig:phi_vs_psnr}, PL interpolation results are close to ground truth when $\phi$ approaches sampled toroidal angles and show worst accuracy when $\phi$ is in the middle of two sampled angles. Also, the local minima of PL are greater than those of SL.

\begin{figure}[!th]
    \centering
    \includegraphics[width=\columnwidth]{figs/time_comp.pdf}
    \caption{
    $\phi$ vs. time for MF, PL, and SL interpolations on XGC data. SL is slightly faster than PL, while MF usually needs $O(10^4)$ times compared with PL and SL. On the sampled poloidal planes, MF takes less than one second because only a 2D cell locating and a 2D barycentric interpolation are required for MF on sampled poloidal planes.
    }
    \label{fig:phi_vs_time}
\end{figure}

\begin{figure}[!th]
    \centering
    \includegraphics[width=\columnwidth]{figs/psnr_comp.pdf}
    \caption{$\phi$ vs. PSNR for PL and SL interpolations on XGC data taking MF interpolation results as ground truth. Both of the curves have 48 infinity values because we upsampled the toroidal coordinate from 16 into 48 poloidal planes to avoid ill-posed prisms, and all three methods do the same interpolation on poloidal planes. The local minima are always attained at the middle of two poloidal planes for both PL and SL, while the PSNR of PL is usually greater than 40 dB while that of SL drops below 30 dB sometimes.}
    \label{fig:phi_vs_psnr}
\end{figure}

\subsection{Volume Rendering (Forward Approximation)}
\label{sec:vol_ren}

We implement and compare three ways to volume render the scalar function $f$ in XGC data with the toroidal coordinate straightened, which transfers $(R,Z,\phi)^\intercal$ into $(x,z,y)^\intercal$, respectively, and renders in Cartesian coordinate system (\Cref{fig:vol_ren}). We can see a rotating pattern of values in both rendering images. However, SL fails to show a continuous pattern. The reason is that PL and SL use different nodes in interpolation at the same location. PL is based on a simplicial mesh of deforming spacetime that relates each mesh node to those close to its location at the previous or next time slot, while SL always uses nodes at the same location for interpolation, which fails to incorporate any deforming property in the underlying physics.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/volRen.pdf}
    \caption{Volume rendering result of $f$. Transfer function is shown on the bottom. PL interpolation shows continuous trajectories of value while SL interpolation breaks the features into dashed lines, because SL builds unmeaningful connectivities between two consecutive timesteps.}
    \label{fig:vol_ren}
\end{figure}

\subsection{Isosurfacing (Inverse Approximation)}

% \HG{Can you talk more here?  What are the challenges without spacetime mesh?  What are the science story?}

As shown in~\Cref{fig:isosurfacing}, isosurfacing results of $f$ show the same rotating pattern as volume rendering results in~\Cref{fig:vol_ren}. Again, PL gives continuous isosurfaces while SL yields ``dashed'' lines. Also, when we increase the isovalue to 0.2, then isosurfaces are mainly concentrated in only a small area that is around the saddle point in the magnetic field in XGC data, which also verifies the correctness of isosurfacing results in~\Cref{fig:isosurfacing}. In order to show more interesting features such as blobs, well-chosen isovalues are required, which is beyond the scope of this work.


\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/isosurface.pdf}
    \caption{Isosurfaces with isovalues $\pm 0.1$ and $\pm 0.2$. Red lines show isosurfaces with positive isovalues while blue lines show isosurfaces with negative isovalues. Well-chosen values can show more interesting blobs, which is beyond the scope of this work.
    % \HG{Why are there artifacts in the first figure?  We do not need to see the zero-levelset.}
    }
    \label{fig:isosurfacing}
\end{figure}

% \subsection{\poincare~plots (Forward approximation) \CR{To be continued..}}

% We demonstrate using our deformed meshing scheme to significantly reduce the computational cost of \poincare~plots in XGC simulations by replacing the magnetic-following interpolation with a piecewise linear approximation induced by our 3D simplicial mesh.  

% In general, \poincare~plots are essential tools for analyzing recurrence patterns in periodical domains.  In magnetic fusion applications, \poincare~plots are essential to evaluate the magnetic confinement of plasma and are notoriously expensive in computation because of the magnetic-following interpolation scheme.  In a tokamak, \poincare~plots are computed by tracing many particles along the total magnetic field $\mathbf{B}=\mathbf{B}_0+\Delta\mathbf{B}$, where $\mathbf{B}_0$ is the static magnetic field that does not change over time, and $\Delta\mathbf{B}$ is the time-dependent magnetic perturbation.  Specific to XGC, $\Delta\mathbf{B}$ is indirectly given by the ``parallel'' component\footnotemark of the magnetic potential $A_{par}$ (or simply $A$), such that the magnetic potential can be written as $\mathbf{A}=A\mathbf{b}_0$, where $\mathbf{b}_0$ is the unit vector of $\mathbf{B}_0$, and we have
% \footnotetext{The magnetic perturbation always aligns with the direction of $\mathbf{B}_0$ (a.k.a the ``parallel'' direction) in XGC.}
% \begin{equation}
%   \Delta \mathbf{B} = \nabla \times \mathbf{A} = \nabla \times (A\mathbf{b}_0) = (\nabla A)\times \mathbf{b}_0 + A\nabla\times\mathbf{b}_0,\label{eq:magnetic}
% \end{equation}
% where $\nabla$ is the gradient operator in 3D (cylindrical coordinate system). %; note that $\mathbf{b}_0$ and $\nabla\times\mathbf{b}_0$ can be precomputed.  
% We denote \emph{magnetic map} $\Phi(R_0, Z_0, \phi_0, \phi)$---the function that maps $(R_0, Z_0, \phi_0)^\intercal$ to the end location at $\phi$---as the solution of the below initial value problem
% \begin{align}
%    \frac{\partial\Phi(R_0, Z_0, \phi_0, \phi)}{\partial\phi} &= \frac{R}{B_\phi}(B_R, B_Z, B_\phi)^\intercal\\
%    \Phi(R_0, Z_0, \phi_0, \phi_0) &= (R_0, Z_0, \phi_0)^\intercal,
% \end{align}
% where $(R_0, Z_0, \phi_0)^\intercal$ is the starting position, and $B_R$, $B_Z$, and $B_\phi$ is the radial, axial, and toroidal component of $\mathbf{B}$, respectively.  % In practice, Runge-Kutta integration methods are used to solve the initial value problem numerically.  
% Formally, the \poincare~plot of one single seed $(R_0, Z_0, 0)^\intercal$ is the set of points $\{\Phi(R_0, Z_0, 0, 2k\pi )\}$, where $k=0, 1, \ldots, n_c-1$ samples $n_c$ cycles.  In practice, one also needs to compute the plot with many (denoted as $n_p$) initial locations to visualize the periodical patterns of particles in a tokamak.

% In this study, we replace the expensive magnetic-following interpolation with a piecewise linear (PL) approximation.  Note that in Eq.~\eqref{eq:magnetic}, both $A$ and $\nabla A$ are evaluated with the magnetic-following interpolation scheme for arbitrary locations during the computation of \poincare~plots, resulting in a prohibitive computational cost that requires a supercomputer to accomplish.  For the query of $A$ and $\nabla A$ at an arbitrary location $(R, Z, \phi)^\intercal$ during the numerical integration, normally, one must interpolate the function values with the magnetic-following interpolation (Eq.~\eqref{eq:interpolation}), which requires tracing a magnetic line of $\mathbf{B}_0$ in both directions and calculating the function values at the intersections at two poloidal planes.  With our 3D simplicial mesh that models the deforming space, one can directly estimate $A$ and $\nabla A$ at arbitrary locations without the need to trace magnetic lines of $\mathbf{B}_0$.

% We compare the resulting \poincare~plots with magnetic-following interpolation and the PL approximation.  We see a significant performance improvement with the PL approximation \HG{TBA}.  We also evaluate the accuracy of our approach by comparing $\Phi(R_0, Z_0, 0, 2\pi)$ (location after the integration with $2\pi$) with two interpolation schemes.  We use each node on XGC's triangular mesh for the evaluation and record the mean-square-root error (MSE) of the output locations.  As a result, we see an MSE of XX and a peak-noise-signal ratio (PSNR) of XXdB.  


\subsection{Extremum Lines (Inverse Approximation)}

\begin{figure}
    \centering
    \includegraphics[width=0.65\columnwidth]{figs/filament.png}
    \caption{Extremum line extraction in XGC data with our deformed mesh.} % \HG{TODO: visualize the same thing with non-deforming mesh}}
    \label{fig:filament}
\end{figure}

To study 3D blob filaments, we co-designed with scientists the definition of \emph{blob core lines} as the extremum lines---loci of local minimum/maximum where radial and axial gradients of $f$ vanish and the radial-axial Hessian $\mathbf{H}_f$ is positive-/negative-definite:
\begin{equation}
  \frac{\partial f}{\partial R} = \frac{\partial f}{\partial Z} = 0 \mbox{~and~} \lambda_1\lambda_2>0,
  % \mathbf{H}(f(\mathbf{x})) \mbox{positive-/negative-definite},
\end{equation}
where $R$ and $Z$ are the radial and axial coordinates, respectively, and 
$\lambda_1$ and $\lambda_2$ are the eigenvalues of the radial-axial Hessian $\mathbf{H}_f$, which considers only $R$ and $Z$ axes:
\begin{equation}
  \mathbf{H}_f = 
  {\begin{pmatrix} \frac{\partial^2 f}{\partial R^2} & \frac{\partial^2 f}{\partial R\partial Z}\\
  \frac{\partial^2 f}{\partial Z\partial R} & \frac{\partial^2 f}{\partial Z^2}
  \end{pmatrix}}.
\end{equation}
% $\mathbf{x}\in\mathbb{R}^4$ are the radial, axial, azimuth, and time coordinates of a point of the blob filament; 
% 
We reformulate the extraction of extremum lines as a critical point tracking problem by treating $\phi$ as time and use the feature tracking kit FTK~\cite{guo2021ftk} to extract and visualize the curves.  FTK assumes that the gradient vector field is piecewise linear so that critical point trajectories are reconstructed directly in a spacetime mesh, which is by default nondeforming; we made minor changes to FTK software to support deformed spacetime meshes.  Figure~\ref{fig:filament} demonstrates extremum line extraction results with our deformed mesh.  Qualitatively speaking, the resulting curves reflect the same trends that scientists observe in isosurface and volume rendering; we plan to investigate further the accuracy evaluation of extremum lines in future work.

% While we observe a significant benefit of using deformed spacetime meshes for feature tracking, 

% and non-deformed mesh.  
% We observe artifacts resulting from non-deforming meshes, and extremum lines appear closer to the isosurfaces of $f$.  

% With assumption that both radial gradient $\partial f / \partial R$ and axial gradient $\partial f / \partial R$ are piecewise linear in $R$, $Z$, and $\phi$ directions. \HG{To be continued.} 


% \subsection{Inverse evaluation in centrifugal pumps: Tracking critical points}
\subsection{Additional Case Study with Synthetic Data}

The additional case study demonstrates the generality of our methodology beyond fusion. The data are synthesized by rotating two Gaussian blobs with phase difference $\pi$ around the center of a circular field. We set the radius of the circular field to be 1 and the full width at half maximum of two blobs to be 0.3, and we uniformly choose 16 timesteps as sampled times in one cycle ($2\pi$). 
The synthetic function is 
% Then the value at position $(x,y)$ and timestep $t$ is given by the following closed form:
\begin{align}
\begin{split}
    f(x,y,t)&=Ce^{-B((x-A\cos\omega t)^2+(y-A\sin\omega t)^2))}\\
    &+Ce^{-B((x+A\cos\omega t)^2+(y+A\sin\omega t)^2))},
\end{split}
\end{align}
% \begin{align}
% \begin{split}
%     f(x,y,t)&=9.806e^{-30.807((x-0.5cos(0.393t))^2+(y-0.5sin(0.393t))^2))} \\
%     &+9.806e^{-30.807((x+0.5cos(0.393t))^2+(y+0.5sin(0.393t))^2))}.
% \end{split}
% \end{align}
% 
where $C=9.806$, $B=30.807$, and $\omega=0.393$.  
Values at several timesteps are shown in \Cref{fig:synthetic}. Discretization of this 2D field gives 1,000 nodes and 1,947 triangle cells (\Cref{fig:synthetic_mesh}~{(a)}). Next nodes of nodes are determined by their positions at the next timestep, as illustrated in \Cref{sec:spacetime_formulation}, which forms a spacetime with 2,000 nodes (\Cref{fig:synthetic_mesh}~{(c)}). Our algorithm yields 5,961 tetrahedron cells for this spacetime. We conduct the same evaluation as in Sections~\ref{sec:upsampling} and \ref{sec:vol_ren} on this synthetic dataset.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/synthetic.pdf}
    \caption{Synthetic dataset showing rotation of two Gaussian blobs, denoted by A and B, around the center of a circular field with cycle $T=16$. This time-varying dataset also introduces a deforming spacetime.}
    \label{fig:synthetic}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/synthetic_mesh.pdf}
    \caption{
    % \HG{It appears that boundary triangles are very thin.  Let me try to find a good circular mesh for the synthetic case.}
    Mesh for synthetic data. (a) and (b) show the original 2D triangular mesh of a circular field with 1,000 nodes and 1,903 cells. Spatial-connected components in lower and upper meshes are differentiated by a categorical colormap in (a) and (b), respectively, while white edges highlight cutting graphs. Every red line connects two matched components. (c) shows the deformed spacetime mesh with 2,000 nodes. Lateral edges indicate the rotation of nodes. All edges are straight lines, although they create a curved visual effect on two sides. The deformed mesh generated by our algorithm yields 5,841 cells for these 2,000 nodes.}
    \label{fig:synthetic_mesh}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figs/synthetic_unsampling.pdf}
    \caption{Ground truth (GT) and upsampling results of our method, piecewise linear (PL) interpolation on the synthetic dataset. PL at $\phi=(21/32)2\pi$ has MSE 0.074 and PSNR 50.718 dB; that at $\phi=(31/32)2\pi$ has MSE 0.191 and PSNR 50.593 dB.}
    \label{fig:synthetic_unsampling}
\end{figure}

\textbf{Temporal upsampling}. We again choose two arbitrary timesteps that are in the middle of two successive sampled timesteps, corresponding to $\phi=(21/32)2\pi$ and $\phi=(31/32)2\pi$ in toroidal coordinates. The comparison of PL interpolation and ground truth is shown in \Cref{fig:synthetic_unsampling}, which shows a high PSNR. \Cref{fig:synthetic_phi_vs_psnr} shows the variation of PSNR over toroidal coordinates on both PL and SL. Since there are 16 sampled timesteps, there are 16 arcs in \Cref{fig:synthetic_phi_vs_psnr}. Each arc attains its maxima at two sampled toroidal angles and its minimum at the middle of these two angles. Also, PL has a higher minimum PSNR than does SL, meaning SL gives wrong interpolation results.

\begin{figure}[!th]
    \centering
    \includegraphics[width=\columnwidth]{figs/denoised_comp.pdf}
    \caption{$\phi$ vs. PSNR of PL and SL interpolations on synthetic dataset. We evenly sample 16 timesteps, so there are 16 infinity values on both curves. All other observations agree with those in~\Cref{fig:phi_vs_psnr}.}
    \label{fig:synthetic_phi_vs_psnr}
\end{figure}

\begin{figure}[!th]
    \centering
    \includegraphics[width=0.75\columnwidth]{figs/synthetic_volRen.pdf}
    \caption{Volume rendering results and isosurfaces for PL and SL interpolations on synthetic data. Similar to~\Cref{fig:vol_ren} and \Cref{fig:isosurfacing}, PL yields continuous features while SL breaks the features.}
    \label{fig:synthetic_vol}
\end{figure}

\textbf{Volume rendering and isosurfacing}. Results for both PL (our method) and SL (naive interpolation) are shown in \Cref{fig:synthetic_vol}. % It is obvious that 
PL reflects the correct rotation of two blobs, but SL gives a broken translation along the time axis.


\section{Discussion}
\label{sec:discussion}

% \textbf{Observations}. 
Simplicial mesh for deforming spacetime makes it possible to derive function values at any given 2D position and any given time $t$ simply by linear interpolation. It gives results with comparable quality (measured by MSE and PSNR) compared with those given by physical derivation while taking significantly less time. It also shows why the naive SL interpolation is not applicable: our method, PL interpolation, builds reasonable connectivities between two successive timesteps, whereas  SL interpolation relates a 3D point with irrelevant mesh nodes. Furthermore, it enables the inverse approximation, a solution for root-finding problem, which cannot be achieved by tracing magnetic lines.

%GAIL - enables  a problem seems odd; perhaps enable inveerse approximation or root-finding, which...? 
\textbf{Limitations}. Although this work does not require the identity of spatial meshes over time or any properties regarding the connectivity of next nodes, there are still several limitations.  First, the decision tree search can be improved, such as by extending the search of new connections to $n$-hop neighbors, finding a better heuristic for choosing a pivot node, and parallelizing the branch search to reduce its time complexity for the worst case. Second, the indivisible problem is not completely solved. For example, we do not triangulate ill-posed prisms but instead increase the temporal resolution to avoid them. Also, this work provides no methods to check indivisible polyhedra except for prisms without dividing them first. Third, we do not continue with figuring out the optimal balance between the number of nodes in cycles formed by cutting faces and the complexity of triangulating and coordinating the cutting faces. Fourth, the strict complexity of the whole algorithm is challenging to analyze because of the uncertainty of how many loops one would take to deal with the indivisible polyhedra.

\textbf{Future work}. Besides solving the previously mentioned limitations, this work has several possible extensions. One possible way is to extend the problem from 3D deforming spacetime to multi-D deforming spacetime. For example, if we have a 3D tetrahedron mesh as a spatial mesh, then it extrudes a 4D spacetime. 
%GAIL - what does they refer to?
Triangulating a 4-polytope into 4-simplices with all similar constraints involves more problems. Another promising work could be extending the triangular spatial mesh to mesh with any polygon as cells. Additional work is triangulating the polygonal cells while guaranteeing that the whole spacetime is still triangulable. A third future work could be using other interpolation methods, such as neural networks, to increase interpolation quality further. Piecewise-linear interpolation may not be the best model of value distribution in cells. Although we can increase temporal resolution to improve interpolation accuracy, upsampling via magnetic-following interpolation is still time-consuming.
% Also, considering spacetime generated by more kinds of spatial meshes can incur new problems, such as adaptive mesh refinement (AMR)~\cite{buschmann2022dark}, scale-space~\cite{agustsson2020scale}, and stellarator~\cite{moritaka2021improving}.

% \HG{I would expect to see a discussion section, which summarize the observations from the case studies.  Also add more discussion on limitations.}

\section{Conclusion}
\label{sec:conclusion}

We propose an algorithm based on the divide-and-conquer paradigm to triangulate 3D nonconvex deforming spacetime with geometric and connectivity constraints. Specifically, in the divided stage we split the spacetime into smaller 3D partitions to reduce time complexity. Cutting faces are well chosen to guarantee the independence of subdividing any partition; communication between partitions on cross sections is also taken care of. In the conquer stage, we eliminate nodes in a 3D partition one by one through a decision-tree-based search. This tree makes decisions on which node to eliminate and how to eliminate it. Our algorithm provides a scheme for and reduces the complexity of many visualization and analysis tasks, such as upsampling, volume rendering, and isosurfacing, on time-dependent datasets, especially those with rotating features. We evaluate the algorithm both quantitatively and qualitatively on XGC data and a synthetic dataset, which verifies its high accuracy and low time complexity compared with the magnetic-line-following method and naive spacetime interpolation.


% % \appendices
% \section{Proof of the First Zonklar Equation}
% Appendix one text goes here.

% % you can choose not to have a title for an appendix
% % if you want by leaving the argument blank
% \section{}
% Appendix two text goes here.


% use section* for acknowledgment
\ifCLASSOPTIONcompsoc
  % The Computer Society usually uses the plural form
  \section*{Acknowledgments}
\else
  % regular IEEE prefers the singular form
  \section*{Acknowledgment}
\fi

The authors thank Drs. Choong-Seock Chang, Michael Churchill, Robert Hager, Seung-Hoe Ku, Zeyu Guo, and Rephael Wenger for insightful discussions. 
This research is supported by DOE DE-SC0022753, NSF OAC-2311878, NSF OAC-2313123, NSF IIS-1955764.  This research is also supported by the Exascale Computing Project (ECP), project number 17-SC-20-SC, a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration.  It is also supported by the U.S. Department of Energy, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357.


% Can use something like this to put references on a page
% by themselves when using endfloat and the captionsoff option.
% \ifCLASSOPTIONcaptionsoff
%   \newpage
% \fi

% \newpage
% \nocite{*}
\bibliographystyle{abbrv-doi-hyperref}
%\bibliographystyle{abbrv-doi-hyperref-narrow}
%\bibliographystyle{abbrv-doi}
%\bibliographystyle{abbrv-doi-narrow}

\bibliography{references}

\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{photos/cren.eps}}]{Congrong Ren} is a PhD student at the Department of Computer Science and Engineering in The Ohio State University. She received her B.E. degree in Electrical Engineering from Xi'an Jiaotong University. Her research interests include feature tracking, topological data analysis, and data compression for scientific data.
\end{IEEEbiography}

\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{photos/hguo.eps}}]{Hanqi Guo} is an Associate Professor at the Department of Computer Science and Engineering in The Ohio State University.  He also holds a joint appointment at the Mathematics and Computer Science Division at Argonne National Laboratory.  His research interests include data analysis, visualization, and machine learning for scientific data.  He is an awardee of the U.S. DOE Early Career Research Program (ECRP) in 2022 and received multiple best paper awards in premiere visualization conferences. He received the BS degree in mathematics and applied mathematics from the Beijing University of Posts and Telecommunications in 2009 and the PhD degree in computer science from Peking University in 2014.
\end{IEEEbiography}

% You can push biographies down or up by placing
% a \vfill before or after them. The appropriate
% use of \vfill depends on what kind of text is
% on the last page and whether or not the columns
% are being equalized.

%\vfill

% Can be used to pull up biographies so that the bottom of the last one
% is flush with the other column.
%\enlargethispage{-5in}


% that's all folks
\end{document}


