% ---------------------------------------------------------------------------
% Author guideline and sample document for EG publication using LaTeX2e input
% D.Fellner, v1.20, Jan 18, 2023

\documentclass{egpubl}
\usepackage{eg2026}
 
\Tutorial               % uncomment for Tutorial contribution

% !! *please* don't change anything above
% !! unless you REALLY know what you are doing
% ------------------------------------------------------------------------
\usepackage[T1]{fontenc}
\usepackage{dfadobe}  
\usepackage{bm}
\usepackage{dirtytalk}

\usepackage{cite}  % comment out for biblatex with backend=biber

\usepackage{amssymb, amsmath, mathtools
% , amsthm
}
\usepackage{latexsym}
\usepackage{framed,multirow, lineno, url}
\usepackage{epsfig}

% ---------------------------
%\biberVersion
\BibtexOrBiblatex
%\usepackage[backend=biber,bibstyle=EG,citestyle=alphabetic,backref=true]{biblatex} 
%\addbibresource{egbibsample.bib}
% ---------------------------  
\electronicVersion
\PrintedOrElectronic
% for including postscript figures
% mind: package option 'draft' will replace PS figure by a filename within a frame
% \ifpdf \usepackage[pdftex]{graphicx} \pdfcompresslevel=9
% \else \usepackage[dvips]{graphicx} \fi

\usepackage{egweblnk} 
% end of prologue
% ---------------------------------------------------------------------
% EG author guidelines plus sample file for EG publication using LaTeX2e input
% D.Fellner, v2.04, Dec 14, 2023

\usepackage{cleveref}

\newcommand{\mathR}{\mathbb{R}}
\newcommand{\XI}{{\mathcal X}}
\newcommand{\bu}{{\mathbf u}}
\newcommand{\bp}{{\mathbf p}}
\newcommand{\bn}{{\mathbf n}}
\newcommand{\bX}{{\mathbf X}}
\newcommand{\bY}{{\mathbf Y}}
\newcommand{\Lag}{\mathrm Lag}
\newcommand{\entreguillemets}[1]{``{#1}''}

\newcommand\bx{{\mathbf x}}
\newcommand\by{{\mathbf y}}
\newcommand\bv{{\mathbf v}}
\newcommand\bq{{\mathbf q}}
\newcommand\bg{{\mathbf g}}
\newcommand\bbf{{\mathbf f}}
\newcommand\bK{{\mathbf K}}
\newcommand\bpsi{\psi}

\newcommand{\cph}[1]{\textcolor{red}{#1}}
% \newcommand{\cph}[1]{}
\newcommand{\newtext}[1]{\textcolor{blue}{#1}}

%\newcommand{\eqref}[1]{\ref{#1}}

\newtheorem{observation}{Observation}
\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}

\def \Vor {\mbox{Vor}}
\def \Pow {\mbox{Pow}}
\def \Del {\mbox{Del}}
\def \Reg {\mbox{Reg}}


\title[Optimal Transport for Fluid Simulation: New and Old]%
      {Optimal Transport for Fluid Simulation: New and Old}

% for anonymous conference submission please enter your SUBMISSION ID
% instead of the author's name (and leave the affiliation blank) !!
% for final version: please provide your *own* ORCID in the brackets following \orcid; see https://orcid.org/ for more details.
\author[C. Plateau Holleville \& B. Lévy]
{\parbox{\textwidth}{\centering C. Plateau Holleville$^{1}$\orcid{0000-0003-1510-557X}
        and B. Lévy$^{1, 2}$\orcid{0000-0002-7007-3219}
        }
        \\
% For Computer Graphics Forum: Please use the abbreviation of your first name.
{\parbox{\textwidth}{\centering 
         $^1$INRIA Saclay, France\\
         $^2$ Université Paris-Saclay, CNRS, Laboratoire de Maths d'Orsay\\
       }
}
}
% ------------------------------------------------------------------------

% if the Editors-in-Chief have given you the data, you may uncomment
% the following five lines and insert it here
%
% \volume{36}   % the volume in which the issue will be published;
% \issue{1}     % the issue number of the publication
% \pStartPage{1}      % set starting page


%-------------------------------------------------------------------------
\begin{document}

\maketitle
%-------------------------------------------------------------------------
\begin{abstract}
Recently introduced, numerical methods for fluid simulation relying on Optimal Transport offer new guarantees as compared to classical algorithms. The fluid is modeled with a Lagrangian point of view, directly tracking the motion of the matter. In contrast with previous methods, they are able to enforce volume preservation at each step of the simulation through the resolution of an Optimal Transport problem. Interestingly, this formulation enables the use of an unfamiliar geometric representation, characterized by a generalized Voronoi diagram, offering new possibilities. This tutorial will start with the theoretical foundations of these numerical methods and present the main computational strategies and applications.

\begin{CCSXML}
<ccs2012>
<concept>
<concept_id>10010147.10010371.10010352.10010379</concept_id>
<concept_desc>Computing methodologies~Physical simulation</concept_desc>
<concept_significance>500</concept_significance>
</concept>
<concept>
<concept_id>10003752.10010061.10010063</concept_id>
<concept_desc>Theory of computation~Computational geometry</concept_desc>
<concept_significance>300</concept_significance>
</concept>
</ccs2012>
\end{CCSXML}

\ccsdesc[500]{Computing methodologies~Physical simulation}
\ccsdesc[300]{Theory of computation~Computational geometry}

\printccsdesc   

\end{abstract}  

\section{Introduction}

This document is a short introduction to Optimal Transport for the simulation of fluids. It accompanies the hands-on tutorial given at Eurographics 2026, in addition to the programming material available from \url{https://github.com/PlathC/pygeogram}. 

It focuses on the theoretical and algorithmic machinery behind the scenes. It is mostly based on \cite{Levy2015} %(preliminary tech. report \cite{Levy2015a}), \cite{Levy2018} %(preliminary tech. report \cite{Levy2018a}) 
and \cite{Levy2022}. %(preliminary tech. report \cite{Levy2022a}). 
This is made with the goal of letting the reader grasp an intuition of general Optimal Transport theory, how it translates to effective algorithms in the semi-discrete setting and how the specific context of free-surface fluids can be formulated as a so-called Partial Optimal Transport problem.

Before entering the heart of the matter, we mention that fluid simulation is a mature field of research, with many different numerical methods that can simulate complex behaviors of fluids and deformable materials, such as \emph{Smooth Particle Hydrodynamics}~(SPH)~\cite{Gingold1977}, \emph{Lattice-Boltzmann Methods}~(LBM)~\cite{McNamara1988} or \emph{Material Point Method}~(MPM)~\cite{Sulsky1995} to name a few. An extensive survey is available in two Eurographics state-of-the-art reports~\cite{Ihmsen14, Koschier2022}.

Recently, Optimal Transport appeared as an alternative to simulate fluids. It is not as mature as the families of methods mentioned above, but it has the interesting property of defining a decomposition of the fluid volume into volumetric cells that follow the motion of the fluid (Lagrangian) while precisely controlling the volume of those cells. Interestingly, in the so-defined decomposition, the same representation is used for all the interfaces (fluid-boundary interactions, contact between different fluids and free surfaces). Besides these exotic methods for fluid simulations, the set of algorithms presented here has a potential impact in geometry processing, machine learning, and other fields of numerical simulations such as topology optimization.

The presentation stays at an elementary level, which corresponds to a computer scientist's
vision of the problem. In the article, we stick to using standard notions of analysis
(functions, integrals) and linear algebra (vectors, matrices), and give an intuition of
the notion of measure. The main objective of the presentation is to understand the overall
structure of the reasoning (\say{Teach principles, not equations.} [R. Feynmann]), and
to follow a continuous path from the theory to an efficient
algorithm that can be implemented in a computer.  


\section{Optimal Transport}

Optimal transport, initially studied by Monge, \cite{Monge1784}, is a
very general mathematical framework that can be used to model a wide
class of application domains. In particular, it is a natural
formulation for several fundamental questions in computer graphics
\cite{DBLP:journals/focm/Memoli11,DBLP:journals/cgf/Merigot11,DBLP:journals/tog/BonneelPPH11},
because it makes it possible to define new ways of \emph{comparing} functions,
of measuring \emph{distances} between functions and
\emph{interpolating} between two (or more) functions~:

\begin{figure}[b]
      \centerline{
         \includegraphics[width=\columnwidth]{Figures/distance_fonction}
      }
      \caption{
        Comparing functions: one would like to say that $f_1$ is nearer to $f_2$ than $f_3$, but
        the classical $L_2$ norm \entreguillemets{does not see} that the graph of $f_2$ corresponds
        to the graph of $f_1$ slightly shifted along the $x$ axis.
      }
      \label{fig:distance_fonction}
\end{figure}  

\paragraph*{Comparing functions:} Consider the functions $f_1$, $f_2$ et $f_3$ in Figure \ref{fig:distance_fonction}.
Here we have chosen a function $f_1$ with a wildly oscillating graph, and a function $f_2$ obtained by translating the graph
of $f_1$ along the $x$ axis. The function $f_3$ corresponds to the mean value
of $f_1$ (or $f_2$). If one measures the relative distances between these functions using the classical $L_2$ norm,
that is  $d_{L_2}(f,g) = \int (f(x) - g(x))^2 dx$, one will find that $f_1$ is nearer to $f_3$ than $f_2$. Optimal transport
makes it possible to define a distance that will take into account that the graph of $f_2$ can be obtained from $f_1$ through
a translation (like here), or through a deformation of the graph of $f_1$. From the point of view of this new distance,
the function $f_1$ will be nearer to $f_2$ than to $f_3$.

\begin{figure}
      \centerline{
         \includegraphics[width=\columnwidth]{Figures/transport_1d}
      }
      \caption{
        Interpolating between two functions: linear interpolation makes a hump disappear
        while the other hump appears; displacement interpolation, stemming from optimal transport,
        will displace the hump as expected.
      }
      \label{fig:interpolation}
\end{figure}  

\paragraph*{Interpolating between functions:} Now consider the functions $u$ and $v$
in Figure \ref{fig:interpolation}. Here we suppose that $u$ corresponds
to a certain physical quantity measured at an initial time $t=0$ and that $v$
corresponds to the same phenomenon measured at a final time $t=1$. The problem
that we consider now consists in reconstructing what happened between $t=0$ and
$t=1$. If we use linear interpolation (Figure \ref{fig:interpolation},top-right),
we will see the left hump progressively disappearing while the right hump progressively
appears, which is not very realistic if the functions represent for instance a propagating
wave. Optimal transport makes it possible to define another type of interpolation
(Mc. Cann's displacement interpolation, Figure
\ref{fig:interpolation},bottom-right), that will progressively displace and
morph the graph of $u$ into the graph of $v$.

Optimal transport makes it possible to define a \emph{geometry} of
a space of functions (or more general objects, called probability
measures, more on this later), and thus gives a definition of \emph{distance}
in this space, as well as means of interpolating between different functions,
and in general, defining the barycenter of a weighted familly of functions, in
a very general context. Thus, optimal transport appears as a fundamental tool in many
applied domains. In computer graphics, applications were proposed, to compare and
interpolate objects of diverse natures \cite{DBLP:journals/tog/BonneelPPH11}, to generate
lenses that can concentrate light to form caustics in a prescribed manner
\cite{DBLP:journals/corr/abs-1708-04820,
  DBLP:journals/tog/SchwartzburgTTP14}. Moreover, optimal transport
defines new tools that can be used to discretize Partial Differential
Equations, and define new numerical solution mechanisms \cite{JKO2}.
This type of numerical solution mechanism can be used to simulate
for instance fluids \cite{Gallouet2017}, with spectacular
applications and results in computer graphics
\cite{DBLP:journals/tog/GoesWHPD15}.

\label{chap5:sec:ot:OT}

The two sections that follow are partly inspired by \cite{OTON}, \cite{OTintro}, \cite{MAEintro} and \cite{OTuserguide},
but stays at an elementary level. Here the main goal is to give an intuition of the different concepts,
and more importantly an idea of the way the relate together. Finally we will see how they can be
directly used to design a computational algorithm with very good performance, that can be used in practice
in several application domains.

% TODO: annonce du plan.

\begin{figure}[b]
      \centerline{
         \includegraphics[width=\columnwidth]{Figures/monge}
      }
      \caption{
        Given two terrains defined by their height functions $u$ and $v$, symbolized here
        as gray levels, Monge's problem consists in transforming one terrain into the other one
        by moving matter through an application $T$. This application needs to satisfy a mass
        preservation constraint.
      }
      \label{fig:monge}
\end{figure}  

\section{Monge's problem}

The initial optimal transport problem was first introduced and
studied by Monge, right before the French revolution \cite{Monge1784}.
We first give an intuitive idea of the problem, then quickly introduce
the notion of \emph{measure}, that is necessary to formally state the
problem in its most general form and to analyze it.

\subsection{Intuition} Monge's initial motivation to study this problem
was very practical~: supposing you have an army of workmen, how can you
transform a terrain with an initial landscape into a given desired target
landcape, while minimizing the total amount of work?

Monge's initial problem statement was as follows:
\begin{equation}
\begin{array}{l}
   \inf\limits_{T~: X \rightarrow X} \int\limits_X c(x, T(x)) u(x) dx \\[6mm]
   \mbox{\ \ \ \  subject to~: } \\[6mm]
   \forall B \subset X, \int\limits_{T^{-1}(B)} u(x)dx = \int\limits_B v(x)dx
\end{array}
\label{eqn:Monge}
\end{equation}
where $X$ is a subset of ${\mathbb R}^2$, $u$ and $v$ are two positive functions
defined on $X$ and such that $\int_X u(x) dx$ = $\int_Y v(x)dx$, and $c(.,.)$ is
a convex distance (the Euclidean distance in Monge's initial problem statement). 

The functions $u$ and $v$ represent the height of the current lanscape and the height
of the target landscape respectively (symbolized as gray levels in Figure \ref{fig:monge}).
The problem consists in finding (if it exists) a function $T$ from $X$ to $X$ that
transforms the current landscape $u$ into the desired one $v$, while minimizing the product
of the amount of transported earh $u(x)$ and the distance $c(x,T(x))$ to which it was transported.
Clearly, the amount of earth is preserved during transport, thus the total quantity of earth should
be the same in the source and target landscapes (the integrals of $u$ and $v$ over $X$ should coincide).
This \emph{global} matter preservation constraint needs to be completed by a \emph{local} one. The local
matter preservation constraint enforces that in the target landscape, the quantity of earh received in
any subset $B$ of $X$ corresponds to what was transported here, that is the quantity of earth initially
present in the preimage $T^{-1}(B)$ of $B$ under $T$. Without this constraint, one could locally create
matter in some places and anhihilate matter other places in a counterbalancing way. A map $T$
that satisfies the local mass preservation constraint is called a \emph{transport map}.

\begin{figure}
      \centerline{
         \includegraphics[width=5cm]{Figures/discrete1}
      }
      \caption{
         Transport from a function (gray levels) to a discrete pointset
         (bleue disks).
      }
      \label{fig:discrete1}
\end{figure}  

\subsection{Monge's problem with measures}
\label{sect:monge_measure}
We now suppose that instead of a \entreguillemets{target landscape}, we wish to
transport earth (or a resource) towards a set of points (that will be denoted by $Y$ for now on), that
represent for instance a set of factories that exploit a resource, see 
Figure \ref{fig:discrete1}. Each factory wishes to receive a certain quantity of resource
(depending for instance of the number of potential customers around the factory). Thus, the
function $v$ that represents the \entreguillemets{target landscape} is replaced with a function
on a finite set of points. However, if a function $v$ is zero everywhere except on a finite set
of points, then its integral over $X$ is also zero. This is a problem, because for instance
one cannot properly express the mass preservation constraint. For this reason, the notion
of function is not rich enough for representing this configuration. One can use instead
\emph{measures} (more on this below), and associate with each factory a \emph{Dirac mass} weighted by the quantity
of resource to be transported to the factory. 

From now on, we will use measures $\mu$ and $\nu$ to represent the
\entreguillemets{current landscape} and the \entreguillemets{target
  landscape}.  These measures are supported by sets $X$ and $Y$, that
may be different sets (in the present example, $X$ is a subset of
${\mathbb R}^2$ and $Y$ is a discrete set of points). Using measures
instead of function not only makes it possible to study our
\entreguillemets{transport to discrete set of factories} problem, but
also it can be used to formalize computer objects (meshes) and
directly leads to a computational algorithm. This algorithm is very
elegant because it is a \emph{verbatim computer translation of the
  mathematical theory} (see \S\ref{sect:algo}).  In this particular
setting, translating from the mathematical language to the algorithmic
setting does not require to make any approximation. This is made
possible by the generality of the notion of measure.

The reader who wishes to learn more on measure theory will refer to
the reference textbooks \cite{TaoMeasureTheory}. To keep the length
of this article reasonable, we will not give here the formal definition
of a measure. In our context, one can think of a measure as a
\entreguillemets{function} that can be only queried using integrals
and that can be \entreguillemets{concentrated} on very small sets (points).
The following table can be used to intuitively translate from the
\entreguillemets{langage of functions} to the \entreguillemets{language of measures}~:
$$
\begin{array}{|c|c|}
    \hline
    \mbox{Function $u$}    & \mbox{Measure $\mu$ } \\[2mm]
    \hline    
    \int_B u(x) dx      & \mu(B) \mbox{ or } \int_B d\mu \\[2mm]
    \int_B f(x) u(x) dx & \int_B f(x) d\mu \\[2mm]
    u(x)                & N/A \\
    \hline    
\end{array}
$$
(Note: in contrast with functions, measures cannot be evaluated at a point,
 they can be only integrated over domains).
 
\begin{figure*}
\includegraphics[width=\textwidth]{Figures/OTplans}
\caption{
  Four example of transport plans in 1D. A: a segment is
  translated. B: a segment is split into two segments.
  C: a Dirac mass is split into two Dirac masses; D: a Dirac
  mass is spreaded along two segments. The first two examples
  (A and B) have the form $(Id \times T)\sharp \mu$ where $T$ is a transport
  map. The third and fourth ones (C et D) have no corresponding transport map,
  because each of them splits a Dirac mass.} 
\label{fig:ch5:ot:OTplans}
\end{figure*}

In its version with measures, Monge's problem can be stated as follows:

\begin{equation}
  (M) \quad \quad
  \begin{array}{l}
    \inf\limits_{T~:X\rightarrow Y} \int\limits_X c(x,T(x)) d\mu
     \mbox{ subject to: } \nu = T\sharp\mu
    \end{array}
    \label{eqn:monge}
\end{equation}


where $X$ and $Y$ are Borel sets (that is, sets that can be measured),
$\mu$ and $\nu$ are two measures on $X$ and $Y$ such that $\mu(X)
= \nu(Y)$, and $c(.,.)$ is a convex distance.  The constraint $\nu =
T\sharp\mu$, that reads \entreguillemets{$T$ pushes $\mu$ onto $\nu$}
corresponds to the local mass preservation constraint. Given a
measure $\mu$ on $X$ and a map $T$ from $X$ to
$Y$, the measure $T\sharp\mu$ on $Y$, called \entreguillemets{the
pushforward of $\mu$ by $T$}, is such that $T\sharp\mu(B)
= \mu(T^{-1}(B))$ for all Borel sets $B \subset Y$. Thus, the local mass
preservation constraint means that $\mu(T^{-1}(B))
= \nu(B)$ for all Borel set $B$ $\subset$ $Y$.

The local mass conservation constraint makes the problem very
difficult: imagine now that you want to implement a computer
program that enforces it: the constraint concerns
\emph{all the subsets} $B$ of $Y$. Could you imagine an algorithm
that just \emph{tests} whether a given configuration satisfies it?
What about \emph{enforcing} it? We will see below a series of
transformations of the initial problem into equivalent problems,
where the constraint become \emph{linear}. We will finally end up
with a simple convex optimization problem, that can be solved
numerically using classical methods.

\begin{figure}[b]
      \centerline{
         \includegraphics[width=\columnwidth]{Figures/monge_pb}
      }
      \caption{A classical example of the existence problem: there is no optimal tansport between
        a segment $L_1$ and two parallel segments $L_2$ and $L_3$
      (it is always possible to find a better transport by replacing $h$ with $h/2$).      
      }
      \label{fig:monge_pb}
\end{figure}  


Before then, let us get back to examine the original problem. The local mass preservation
constraint is not the only difficulty: the functional optimized by Monge's problem is
non-symmetric, and this causes additional difficulties when studying the existence
of solutions for problem (M). The problem is not symmetric because $T$ needs to
be a map, therefore the source and target landscape do not play the same role.
Thus, it is possible to merge earth (if $T(x_1) = T(x_2)$ for two different points
$x_1$ and $x_2$), but it is not possible to split earth (for that, we would need a
\entreguillemets{map} $T$ that could send the same point $x$ to two different points $y_1$ and $y_2$).
The problem is illustrated in Figure \ref{fig:monge_pb}:
suppose that you want to compute the optimal transport between a segment
$L_1$ (that symbolizes a \entreguillemets{wall of earth}) and two parallel
segments $L_2$ and $L_3$ (that symbolize two \entreguillemets{trenches} with a depth that
correspond to half the heigth of the wall of earth). Now we want to transport the wall of earth
to the trenches, to make the landscape flat. To do so, it is possible to decompose $L_1$ into
segments of length $h$, sent alternatively towards $L_2$ and $L_3$
(Figure \ref{fig:monge_pb} on the left). For any length
$h$, it is always possible to find a better map $T$, that is a lower value
of the functional in \eqref{eqn:monge}, by subdividing $L_1$ into
smaller segments (Figure \ref{fig:monge_pb} on the right). The
best way to proceed consists in sending from each point of
$L_1$ \emph{half} the earth to $L_2$ and \emph{half} the earth to $L_3$, which
cannot be represented by a map. Thus, the best solution of problem (M) is not a
map. In a more general setting, this problem appears each time the source measure
$\mu$ has mass concentrated on a manifold of dimension $d-1$ \cite{McCannEUMMP05}
(like the segment $L_1$ in the present example).


\section{Kantorovich's relaxed problem}

To overcome this difficulty, Kantorovich stated a problem with a
larger space of solutions, that is, a relaxation of problem (M),
where mass can be both split and merged. The idea consists
in solving for the \entreguillemets{graph of $T$} instead of $T$.
One may think of the graph of $T$ as a function $g$ defined on $X \times Y$
that indicates for each couple of points $x\in X, y\in Y$ how much matter goes from
$x$ to $y$. However, once again, we cannot use standard functions to represent the
graph of $T$: if you think about the graph of a univariate function $x \mapsto f(x)$,
it is defined on $\mathbb{R}^2$ but \emph{concentrated} on a curve. For this reason, as in our
previous example with factories, one needs to use measures. Thus, we are now 
looking for a measure $\gamma$ supported by the
product space $X \times Y$. The relaxed problem is stated as follows:

\begin{equation}
(K) \quad \quad
   \begin{array}{l}
     \inf\limits_{\gamma} \left\{ \int\limits_{X \times Y} c(x,y) d \gamma \ | \ \gamma \ge 0 \mbox{ et } \gamma \in \Pi(\mu,\nu) \right\} \\[5mm]
   \mbox{where: } \\[3mm]
     \Pi(\mu,\nu) = \{ \gamma \in {\mathcal P}(X \times Y) \ | \
     (P_X)\sharp\gamma = \mu \ ; \ (P_Y)\sharp\gamma = \nu \}
   \end{array}
   \label{eqn:K}
\end{equation}
where $(P_X)$ and $(P_Y)$ denote the two projections $(x,y) \in X \times Y \mapsto x$ and
$(x,y) \in X \times Y \mapsto y$ respectively.

The two measures $(P_X)\sharp\gamma$ and $(P_Y)\sharp\gamma$ obtained by pushing forward $\gamma$ by the two projections
are called the \emph{marginals} of $\gamma$. The measures $\gamma$ in the admissible set $\Pi(\mu,\nu)$,
that is, the measures that have $\mu$ and $\nu$ as marginals, are called \emph{optimal transport plans}. Let us now have a
closer look at the two constraints on the marginals $(P_X)\sharp\gamma = \mu$ and $(P_X)\sharp\gamma$
that define the set of optimal transport plans $\Pi(\mu,\nu)$. Recalling the definition of the pushforward (previous subsection),
these two constraints can also be written as:
\begin{equation}
   \begin{array}{lcl}
       (P_X)\sharp\gamma = \mu & \iff & \forall B \subset X, \int_B d\mu = \int_{B \times X} d\gamma \\[2mm]
       (P_Y)\sharp\gamma = \nu & \iff & \forall B^\prime \subset Y, \int_{B^\prime} d\nu = \int_{X \times B} d\gamma       
   \end{array}
\label{eqn:ch5:ot:marginales}   
\end{equation}
Intuitively, the first constraint $(P_X)\sharp\gamma = \mu$ means that everything that comes from a subset $B$
of $X$ should correspond to the amount of matter (initially) contained by $B$ in the source landscape, and the second one
$(P_Y)\sharp\gamma = \nu$ means that everything that goes into a subset $B^{\prime}$ of $Y$ should correspond to the
(prescribed) amount of matter contained by $B^{\prime}$ in the target landscape $\nu$. 


\begin{figure*}
      \centerline{
         \includegraphics[width=\textwidth]{Figures/DK}
      }
      \caption{A discrete version of Kantorovich's problem.}
      \label{fig:DK}
\end{figure*}  

We now examine the relation between the relaxed problem (K)
\eqref{eqn:K} and the initial problem (M)
\eqref{eqn:monge}.  One can easily check that among the optimal
transport plans, those with the form $(Id \times T)\sharp \mu$
correspond to a transport map $T$:

\begin{observation}
  If $(Id \times T)\sharp\mu$ is a transport plan, then $T$ pushes $\mu$ onto $\nu$.
\end{observation}
\begin{proof}
 $(Id \times T)\sharp\mu$ is in $\Pi(\mu,\nu)$, thus
  $(P_2)\sharp(Id \times T)\sharp \mu = \nu$, or
  $\left((P_2)\circ(Id \times T)\right)\sharp \mu = \nu$, and finally
  $T\sharp\mu = \nu$
\end{proof}

We can now observe that if a transport plan $\gamma$ has the form $\gamma = (Id \times T)\sharp \mu$, then
problem (K) becomes:
$$
     \mbox{min} \left\{ \int\limits_{X \times Y} c(x,y)d\left( (Id \times T)\sharp \mu \right) \right\} 
\quad = \quad
     \mbox{min} \left\{ \int\limits_X c(x,T(x)) d \mu \right)
$$
(one retreives the initial Monge problem). \\

To further help grasping an intuition of this notion of transport plan,
we show four 1D examples in Figure \ref{fig:ch5:ot:OTplans} (the transport
plan is then $1D \times 1D = 2D$). Intuitively, the transport plan $\gamma$
may be thought of as a \entreguillemets{table} indexed by  $x$ and $y$
that indicates the quantity of matter transported from $x$ to $y$.
More exactly (recall that one cannot evaluate a measure $\gamma$ at a point $(x,y)$, we can just compute integrals with $\gamma$),
the measure $\gamma$ is non-zero on subsets of $X \times Y$ that contain
points $(x,y)$ such that some matter is transported from $x$ to $y$. Whenever
$\gamma$ derives from a transport map $T$, that is if $\gamma$ has the form
$(Id \times T)\sharp \mu$, then we can consider $\gamma$ as the
\entreguillemets{graph of $T$} like in the first two examples (A) and (B). For example (B), the transport map $T$ is undefined at the center of the segment, but it is not a problem, because one can remove from $X$ any subset with zero measure, and in this example there is no measure concentrated in the center of the segment.

The transport plans of the two other examples (C) and (D) have no
associated transport map, because they split Dirac masses.
The transport plan associated with Figure \ref{fig:monge_pb}
has the same nature (but this time in $2D \times 2D = 4D$). It cannot be written
with the form $(Id \times T)\sharp \mu$ because it splits the mass concentrated in $L_1$
into $L_2$ and $L_3$.


Now the theoretical questions are:
\begin{itemize}
  \item when does an optimal transport plan exist?
  \item when does it admit an associated optimal transport map?
\end{itemize}
A standard approach to tackle this type of
existence problem is to find a certain \emph{regularity} both in the functional
and in the space of the admissible transport plans, that is, proving that the
functional is sufficiently \entreguillemets{smooth} and finding a compact set
of admissible transport plans. Since the set of admissible transport plans
contains at least the product measure $\mu \otimes \nu$, it is non-empty, and
existence can be proven thanks to a topological argument that exploits the
regularity of the functional and the compactness of the set. Once the existence
of a transport plan is established, the other question concerns the existence of
an associated transport map. Unfortunately, problem (K) does not directly show the
structure required by this reasoning path. However, one can observe that (K) is a
\emph{linear} optimization problem subject to \emph{linear} constraints. This suggests
using certain tools, such as the dual formulation, that was also developed by Kantorovich.
With this dual formulation, it is possible to exhibit an interesting structure of the
problem, that can be used to answer both questions (existence of a transport plan,
and existence of an associated transport map).

\section{Kantorovich's dual problem}
\label{ch5:sec:ot:KD}

Kantorovich's duality applies to problem (K), in its most
general form (with measures). To facilitate understanding, we will
consider instead a discrete version of problem (K), where the involved
entities are vectors and matrices (instead of measures and operators).
This makes it easy to better discernate the structure of (K), that is
the linear nature of both the functional and the constraints. This
also makes it easier for the reader to understand how to construct
the dual by manipulating simpler objects (matrices and vectors), however
the structure of the reasoning is the same in the general case.

\subsection{The discrete Kantorovich problem}
In Figure \ref{fig:DK}, we show a discrete version of the 1D transport between
two segments of Figure \ref{fig:ch5:ot:OTplans}.
The measures $\mu$ and $\nu$
are replaced with vectors $U = (u_i)_{i=1\ldots m}$ and $V = (v_j)_{i=1\ldots n}$.
The transport plan $\gamma$ becomes a set of coefficients $\gamma_{ij}$. Each 
coefficient $\gamma_{ij}$ indicates the quantity of matter  that will be transported
from $u_i$ to $v_j$.
The discrete Kantorovich problem can be written as follows:
\begin{equation}
\quad \quad
   \min\limits_\gamma <C,\gamma> \mbox{ subject to }
      \left\{
       \begin{array}{l}
         P_1 \gamma = U \\
         P_2 \gamma = V \\
         \gamma_{i,j} \ge 0 \quad \forall i,j
       \end{array}  
      \right.
\label{chap5:eqn:ot:primal}      
\end{equation}
where $\gamma$ is the vector of $\mathbb{R}^{m \times n}$ with all
coefficients $\gamma_{ij}$ (that is, the matrix $\gamma_{ij}$ \entreguillemets{unrolled}
into a vector), and $C$ the vector of $\mathbb{R}^{m \times n}$ with the coefficients $c_{ij}$
indicating the transport cost between point $i$ and point $j$ (for instance, the
Euclidean cost). The objective function is simply the dot product, denoted by $<C,\gamma>$,
of the cost vector $C$ and the vector $\gamma$. The objective function is \textbf{linear in $\gamma$}.
The constraints on the marginals \eqref{eqn:ch5:ot:marginales} impose in this discrete version
that the sums of the $\gamma_{ij}$ coefficients over the columns correspond to the $u_i$ coefficients (Figure \ref{fig:ch5:ot:OTplans}-B)
and the sums over the rows correspond to the $v_j$ coefficients (Figure \ref{fig:ch5:ot:OTplans}-C).
Intuitively, everything that leaves point $i$ should correspond to $u_i$, that is the quantity
of matter initially present in $i$ in the source landscape, and everything that arrives at a point 
$j$ should correspond to $v_j$, that is the quantity of matter desired at $j$ in
the target landscape. As one can easily notice, in this form, both constraints are
 \textbf{linear in $\gamma$}. They can be written with two matrices $P_1$ and $P_2$, of dimensions
$m \times mn$ and $n \times mn$ respectively. 

\subsection{Constructing the Kantorovich dual in the discrete setting}

We introduce, arbitrarily for now, the following function $L$ defined by:
$$
    L(\varphi, \psi) = <C,\gamma> - <\varphi, P_1 \gamma - U> - <\psi, P_2 \gamma - V>
$$
that takes as arguments two vectors, $\varphi$ in ${\mathbb R}^m$ and $\psi$ in
${\mathbb R}^n$. The function $L$ is constructed from the objective function $<C,\gamma>$ from which
we subtracted the dot products of $\varphi$ and $\psi$ with the vectors that correspond to the
degree of violation of the constraints. One can observe that:
$$
\begin{array}{lcl}
  \sup\limits_{\varphi,\psi}[ L(\varphi,\psi) ] & = & <C,\gamma> \mbox{ si } P_1 \gamma = U \mbox{ et } P_2 \gamma = V \\
                                            & = & +\infty \mbox{ otherwise.}
\end{array}
$$
Indeed, if for instance a component $i$ of $P_1\gamma$ is non-zero, one can make $L$ arbitrarily large
by suitably choosing the associated coefficient $\varphi_i$.

Now we consider:
$$
  \inf\limits_{\gamma \ge 0}\left[ \sup\limits_{\varphi, \psi} [L(\varphi, \psi)] \right] =
  \inf\limits_{\begin{tiny}\begin{array}{l} \gamma \ge 0 \\ P_1 \gamma = U \\ P_2 \gamma = V \end{array}\end{tiny}} \left[ <C, \gamma>  \right].
$$
There is equality, because to minimize $\sup[ L(\varphi, \psi) ]$, $\gamma$ has no other choice than
satisfying the constraints (see the previous observation). Thus, we obtain a new expression (left-hand side)
of the discrete Kantorovich problem (right-hand side). We now further examine it, and replace $L$ by its expression:
\begin{eqnarray}
\vspace{-5mm}\label{chap5:eqn:ot:dual1}
\inf\limits_{\gamma \ge 0}\left[ \sup\limits_{\varphi, \psi} 
    \left(\begin{array}{ll}
      <C,\gamma> & - <\varphi, P_1 \gamma - U> \\
                 & - <\psi, P_2 \gamma - V> 
    \end{array}\right)
  \right] & = & \\[3mm]
\label{chap5:eqn:ot:dual2} 
\sup\limits_{\varphi, \psi} \left[ \inf\limits_{\gamma \ge 0}
  \left(\begin{array}{ll}
  <C,\gamma> & - <\varphi, P_1 \gamma - U> \\
             & - <\psi, P_2 \gamma - V> 
  \end{array}\right)
 \right] & = & \\[3mm]
\label{chap5:eqn:ot:dual3} 
\sup\limits_{\varphi, \psi} \left[ \inf\limits_{\gamma \ge 0}
  \left(
    \begin{array}{l}
      <\gamma, C-P_1^t \varphi - P_2^t \psi> + \\[1mm]
      <\varphi, U> + <\psi, V>
    \end{array} 
  \right)
  \right] & = & \\[3mm]
\label{chap5:eqn:ot:dual4} 
 \sup\limits_{\begin{tiny} \begin{array}{c} \varphi, \psi \\ P_1^t \varphi + P_2^t \psi \le C\end{array}\end{tiny}} \left[ <\varphi,U> + <\psi,V> \right]. \quad \quad 
\end{eqnarray}

The first step \eqref{chap5:eqn:ot:dual2} consists in exchanging the $\inf$ and $\sup$. Then we rearrange the terms \eqref{chap5:eqn:ot:dual3}.
By reinterpreting this equation as a constrained optimization problem (similarly to what we did in the previous paragraph), we finally
obtain the constrained optimization problem in \eqref{chap5:eqn:ot:dual4}. In the constraint $P_1^t \varphi + P_2^t \psi \le C$, the inequality
is to be considered componentwise. Finally, the problem \eqref{chap5:eqn:ot:dual4} can be rewritten as:

\begin{equation}
\begin{array}{l}
   \sup\limits_{\varphi, \psi} \left[ <\varphi, U> + <\psi, V>  \right] \\[2mm]
      \mbox{ subject to~: } \varphi_i + \psi_j \le c_{ij} \quad \forall i,j
\end{array}   
\end{equation}

As compared to the primal problem \eqref{chap5:eqn:ot:primal} that depends on $m \times n$ variables (all the coefficients
$\gamma_{ij}$ of the optimal transport plan for all couples of points $(i,j)$), this dual problem depends on
$m + n$ variables (the components $\psi_i$ and $\psi_j$ attached to the source points and target points).
We will see later how to further reduce the number of variables, but before then,
we go back to the general continuous setting (that is, with functions, measures and operators).
    
\subsection{The Kantorovich dual in the continuous setting}
The same reasoning path can be applied to the continuous Kantorovich problem
\eqref{eqn:K}, leading to the following problem (DK):

\begin{equation}
   \begin{array}{l}
   (DK) \quad \quad \sup\limits_{\varphi,\psi} \left[ \int\limits_X \varphi d\mu + \int\limits_Y \psi d\nu \right ] \\[4mm]
   \mbox{ subject to: } \\
     \varphi(x) + \psi(y) \le c(x,y) \quad \forall (x,y) \in X \times Y,
  \end{array}
  \label{eqn:DK}  
\end{equation} where $\varphi$ and $\psi$ are now functions defined on $X$ and $Y$. The functions
$\varphi$ and $\psi$ need to be taken in $L^1(\mu)$ and $L^1(\nu)$. The proof of the equivalence with problem (K)
requires more precautions than in the discrete case, in particular step \eqref{chap5:eqn:ot:dual2}
(exchanging $\mbox{sup}$ and $\mbox{inf}$), that uses a result of convex analysis (due to Rockafellar),
see \cite{OTON} chapter 5.

The classical image that gives an intuitive meaning to this dual problem is to consider that instead of
transporting earh by ourselves, we are now hiring a company that will do the work on our behalf. The
company has a special way of determining the price: the function $\varphi(x)$ corresponds to what they charge for
loading earth at $x$, and $\psi(y)$ corresponds to what they charge for unloading earth at $y$.
The company aims at maximizing its profit (this is why the dual problem is a \entreguillemets{$\mbox{sup}$} rather
than an \entreguillemets{$\mbox{inf}$)},
but it cannot charge more than what it would cost us if we were doing the work by ourselves (hence the constraint).

The existence of solutions for (DK) remains difficult to study, because the set of functions $\varphi,\psi$ that
satisfy the constraint is not compact. However, it is possible to reveal more structure of the problem,
by introducing the notion of \emph{c-transform}, that makes it possible to exhibit a set
of admissible functions with sufficient regularity:

\begin{definition}
  Given a function $f~: X \rightarrow \mathbb{R} \cup \{ -\infty, +\infty \}$, the c-transform $f^c$ is
  defined on $Y$ by:
$$
   f^c(y) := \inf\limits_{x \in X} \left[ c(x,y) - f(x) \right]
$$
\begin{itemize}
   \item If a function $\varphi$ is such that there exists a function $f$ such that $\varphi = f^c$, 
     then $\varphi$ is said to be \emph{c-concave}~;
   \item ${\bm \Psi}_c(X)$ denotes the set of c-concave functions on $X$.
\end{itemize}
\end{definition}

We now show two properties of (DK) that will allow us to restrict the problem to
the class of c-concave functions to search for $\varphi$ and $\psi$:

\begin{observation}
  If the pair $(\varphi,\psi)$ is admissible for (DK), then the pair $(\varphi,\varphi^c)$ is admissible as well.
\end{observation}
\begin{proof}
$$
\begin{array}{l}
    \left\{
        \begin{array}{l}
           \forall(x,y) \in X \times Y, \varphi(x) + \psi(y) \le c(x,y) \\
           \varphi^c(y) = \inf\limits_{x \in X} \left[ c(x,y) - \varphi(x) \right]
        \end{array}
    \right. \\[5mm]
    \begin{array}{lcl}
        \varphi(x) + \varphi^c(y) & = & \varphi(x) + \inf\limits_{x^\prime \in X}\left[  c(x^\prime,y) - \varphi(x^\prime)  \right] \\
                            & \le & \varphi(x) + c(x,y) - \varphi(x) \\
                            & \le & c(x,y)
    \end{array}
\end{array}
$$
\end{proof}


\begin{observation}
   If the pair $(\varphi,\psi)$ is admissible for (DK), then one obtains a \emph{better pair} by replacing $\psi$ with $\varphi^c$:
\end{observation}
\begin{proof}
$$
\left.
   \begin{array}{lcl}
       \varphi^c(y) & = & \inf\limits_{y \in Y} \left[ c(x,y) - \varphi(x) \right] \\
       \forall x \in X, \quad \psi(y) & \le & c(x,y) - \varphi(x)
   \end{array}
\right\} 
   \Rightarrow
   \psi(y) \le \varphi^c(y).
$$
\end{proof}

In terms of the previous intuitive image, this means that by replacing $\psi$ with $\varphi^c$, the company can charge more while
the price remains acceptable for the client (that is, the constraint is satisfied). Thus, we have:
$$
\begin{array}{ll}
   \inf(K)  & = \sup\limits_{\varphi \in {\bm \Psi}_c(X)} \quad \int\limits_X \varphi d\mu + \int\limits_Y \varphi^cd\nu \\[4mm]
            & = \sup\limits_{\psi \in {\bm \Psi}_c(Y)} \quad \int\limits_Y \psi d\mu + \int\limits_X \psi^cd\nu
\end{array}   
$$   

We do not give here the detailed proof for existence. The reader is 
referred to \cite{OTON}, Chapter 4. The idea is that we are now in
a much better situation, since the set of admissible functions
${\bm \Psi}_c(Y)$ is compact, provided that the value of
$\psi$ is fixed at a point of $Y$ in order to suppress invariance
with respect to adding a constant to $\psi$.

The optimal value gives an interesting information, that is the minimum cost of transforming $\mu$ into $\nu$.
This can be used to define a distance between distributions, and also gives a way to compare different distributions,
which is of practical interest for some applications. 

\section{From the Kantorovich dual to the optimal transport map}

\subsection{Subdifferential}

Suppose now that in addition to the optimal cost you want to know the associated way to transform $\mu$ into $\nu$,
in other words, when it exists, the map $T$ from $X$ to $Y$ which associated transport plan $(Id \times T)\sharp\mu$ minimizes
the functional of the Monge problem. A result characterizes
the support of $\gamma$, that is the subset $\partial_c \psi \subset X \times Y$ of the pairs of points $(x,y)$ connected by the transport plan:

\begin{theorem}
$$
\forall (x,y) \in \partial_c \psi, \nabla \psi(x) - \nabla_x c(x,y) = 0
$$
where $\partial_c \psi = \{ (x,y) | \varphi(x) + \psi(y) = c(x,y) \}$ denotes the so-called \emph{subdifferential} of $\psi$.
\label{thm:MongeSol}
\end{theorem}
\begin{proof}
See \cite{OTON} chapters 9 et 10.
\end{proof}

To give an idea of the relation between the subdifferential and the associated transport map $T$, we summarize below
the heuristic argument given at the beginning of the same chapter of \cite{OTON}:
Consider a point $(x,y)$ in the subdifferential $\partial_c \psi$. By definition of
the subdifferential, we have:
\begin{equation}
\varphi(y) + \psi(x) = c(x,y).
\label{eqn:ssdiff1}
\end{equation}

By definition of the c-conjugate, $\varphi(y) = \psi^c(y) = \inf\limits_x c(x,y) - \psi(x)$, thus
$\forall \tilde{x}, \varphi(y) \le c(\tilde{x},y) - \psi(\tilde{x})$, or:
\begin{equation}
\varphi(y) + \psi(\tilde{x}) \le c(\tilde{x},y).
\label{eqn:ssdiff2}
\end{equation}

By substituting \eqref{eqn:ssdiff2} into \eqref{eqn:ssdiff1},
we obtain that $\psi(\tilde{x}) - \psi(x) \le c(\tilde{x},y) - c(x,y)$ for all $\tilde{x}$.

Imagine now that $\tilde{x}$ follows a trajectory parameterized by $\epsilon$ with its origin at $x$. We can
compute the gradient at $x$ with respect to an arbitrary direction $w$ by considering the limit when
 $\epsilon$ tends to zero within the relation $\frac{\psi(\tilde{x}) - \psi(x)}{\epsilon} \le \frac{c(\tilde{x},y) - c(x,y)}{\epsilon}$.
Thus we see that $\nabla \psi(x) \cdot w \le \nabla_x c(x,y) \cdot w$. We can do the same derivation along direction $-w$
instead of $w$, and then we get:
$\forall w, \nabla \psi(x) \cdot w = \nabla_x c(x,y) \cdot w$, thus
$\forall (x,y) \in \partial_c \psi, \nabla \psi(x) - \nabla_x c(x,y) = 0$. \\
\emph{Note~: the derivations above are purely formal and do not consistute
  a proof. The complete proof requires a finer analysis, using a notion
  of generalized differentiability, as well as tools from convex analysis.}

In the case of the $L_2$ cost, that is with $c(x,y) = 1/2 \| x - y \|^2$,
this relation becomes $\forall
(x,y) \in \partial_c \psi, \nabla \psi(x) + y - x = 0$, thus, when
the optimal transport map $T$ exists, it is given by
$T(x) = x - \nabla \psi(x) = \nabla (\|x\|^2/2 - \psi(x))$. Not
only this gives an expression of $T$ in function of $\psi$, which is
of high interest to us if we want to compute the transport explicitly.
In addition, this makes it possible to characterize $T$ as the gradient of a \emph{convex} function (see also
Brenier's polar factorization theorem \cite{BrenierPFMR91}). This
convexity property is interesting, because it means that
two \entreguillemets{transported particles} $x_1 \mapsto T(x_1)$ et
$x_2 \mapsto T(x_2)$ will never collide. We now see how to prove
these two assertions ($T$ gradient of a convex function
and absence of collision) in the case of the $L_2$ transport (with $(c(x,y) = \| x - y \|^2$).

\begin{figure}
      \centerline{
         \includegraphics[width=50mm]{Figures/convex}
      }
      \caption{
         The upper enveloppe of a family of affine functions is a convex function.
      }
      \label{fig:convex}
\end{figure}

\begin{figure*}
      \centerline{
         \includegraphics[width=\textwidth]{Figures/transport_types}
      }
      \caption{
        Different types de transport, with continuous and discrete measures $\mu$ and $\nu$.
      }
      \label{fig:transport_types}
\end{figure*}

\begin{observation}
If $c(x,y)$ = $1/2 \| x - y \|^2$ and $\psi \in {\bm \Psi}_c(X)$, then $\bar{\psi}~: x \mapsto \bar{\psi}(x) = \| x \|^2/2 - \psi(x)$
is a convex function (it is an equivalence if $X = Y = \mathbb{R}^d$).
\label{obs:PhiConvexity}
\end{observation}
\begin{proof}
$$
\begin{array}{lcl}
   \psi(x) & = & \varphi^c(x) \\
           & = & \inf\limits_y \frac{\|x - y\|^2}{2} - \varphi(y) \\
           & = & \inf\limits_y \frac{\|x\|^2}{2} - x \cdot y + \frac{\|y\|^2}{2} - \varphi(y) \\
-\bar{\psi}(x) & = & \psi(x) - \frac{\|x\|^2}{2} = \inf\limits_y \left[ -x \cdot y + \left( \frac{\|y\|^2}{2} - \varphi(y)  \right) \right] \\[4mm]
\bar{\psi}(x) & = & \sup\limits_y \left[ x \cdot y - \left( \frac{\|y\|^2}{2} - \varphi(y) \right) \right]
\end{array}
$$
The function $ x \mapsto x \cdot y - \left( \frac{\|y\|^2}{2} - \varphi(y) \right)$ is affine in $x$, therefore the graph of
$\bar{\psi}$ is the upper envelope of a family of hyperplanes, therefore $\bar{\psi}$ is a convex function (see Figure \ref{fig:convex}).
\end{proof}

\begin{observation}
We now consider the trajectories of two particles parameterized by $t \in [0,1]$, $t \mapsto (1-t)x_1 + t T(x_1)$ and
$t \mapsto (1-t)x_2 + t T(x_2)$. If $x_1 \neq x_2$ and $0 < t < 1$ then there is no collision between the two
particles.
\end{observation}
\begin{proof}
  By contradiction, suppose there is a collision, that is there exists 
  $t \in (0,1)$ and $x_1 \neq x_2$ such that: 
$$
  \begin{array}{l}
  \begin{array}{lcl}
   (1-t)x_1 + tT(x_1) & = & (1-t)x_2 + tT(x_2) \\
   (1-t)x_1 + t \nabla \bar{\psi}(x_1) & = & (1-t)x_2 + t \nabla \bar{\psi}(x_2)
   \end{array} \\
  \begin{array}{lcl}   
   (1-t)(x_1 - x_2)  + t (\nabla \bar{\psi}(x_1) - \nabla \bar{\psi}(x_2)) & = & 0 \\
   \forall v, (1-t) v \cdot(x_1 - x_2)  + t  v \cdot (\nabla \bar{\psi}(x_1) - \nabla \bar{\psi}(x_2)) & = & 0 \\
   \quad \quad \quad \mbox{ take } v = (x_1 - x_2) \\
   (1-t)\| x_1 - x_2 \|^2 + t (x_1 - x_2)\cdot(\nabla \bar{\psi}(x_1) - \nabla \bar{\psi}(x_2)) & = & 0.
  \end{array}
  \end{array}
$$

The last step leads to a contradiction, between the left-hand side is the sum of two strictly positive numbers (recalling the definition of the convexity of   
$\bar{\psi}$~: 
   $\forall x_1 \neq x_2, (x_1-x_2) \cdot (\nabla \bar{\psi}(x_1) - \nabla \bar{\psi}(x_2)) > 0$)
Note that even if there is no collision, the trajectories can cross (see example in \cite{OTON}. It the cost is the Euclidean distance (instead of squared Euclidean distance), the property is stronger and trajectories cannot cross. This comes at the expense of losing the uniqueness of the optimal transport plan\cite{OTON}.
\end{proof}

\section{Continuous, discrete and semi-discrete transport}

The properties that we have presented in the previous sections are true for any couple of
source and target measures $\mu$ et $\nu$, that can derive from continuous functions or
that can be discrete empirical measures (sum of Dirac masses). Figure \ref{fig:transport_types} presents three
configurations that are interesting to study. These configurations have specific properties, that lead to
different algorithms for computing the transport. We give here some indications and references concerning
the \emph{continuous $\rightarrow$ continuous} and \emph{discrete $\rightarrow$ discrete} cases. Then we will
develop the \emph{continuous $\rightarrow$ discrete} case with mode details in the next section. 

\subsection{The continuous $\rightarrow$ continuous case and Monge-Amp\`ere equation}
We recall that when the optimal transport map exists, in the case of the $L_2$ cost
(that is, $c(x,y) = \| x - y \|^2$), it can be deduced from the function $\psi$ by using
the relation $T(x) = \nabla \bar{\psi} = x - \nabla \psi$. The change of variable formula for integration
over a subset $B$ of $X$ can be written as:
\begin{equation}
   \forall B \subset X, \int_B \mu(x) d\mu = \mu(B) = \nu(T(B)) = \int_B \left| \det J_T(x) \right| T(x) d\nu
   \label{eqn:changevar}
\end{equation}
where $J_T$ denotes the Jacobian matrix of $T$ and $\det$ the determinant.

If $\mu$ and $\nu$ have densities $u$ and $v$ respectively, that is
$\forall B, \mu(B) = \int_B u(x)dx$ and $\nu(B) = \int_B v(x)dx$, then one can (formally) consider
\eqref{eqn:changevar} pointwise in $X$:
\begin{equation}
    \forall x \in X, u(x) = \left|\det J_T(x)\right| v(T(x)).
   \label{eqn:pointwise}
\end{equation}
By injecting $T=\nabla\bar{\psi}$ and $J_T = H \bar{\psi}$ into \eqref{eqn:pointwise}, one obtains:
\begin{equation}
\forall x \in X, u(x) = \left| \det H \bar{\psi}(x) \right| v (\nabla \bar{\psi}(x)),
\label{eqn:MAE}
\end{equation}
where $H\bar{\psi}$ denotes the Hessian matrix of $\bar{\psi}$. Equation
\eqref{eqn:MAE} is known ad the \emph{Monge-Amp\`ere equation}. It is a highly non-linear
equation, and its solutions when they exist often present singularities.
This is similar to the eikonal equation, which solution corresponds to the distance field, that has a singularity on the medial axis. Note that the derivation above is purely formel,
and that studying the solutions of the Monge-Amp\`ere equation require using more sophisticated
tools. In particular, it is possible to define several types of weak solutions (viscosity solutions,
solution in the sense of Brenier, solutions in the sense of Alexandrov \ldots). Several algorithms
to compute numerical solutions of the Monge-Amp\`ere equations were proposed. As such, see for
instance the Benamou-Brenier algorithm \cite{ACFM:BB:2000}, that uses a dynamic formulation
inspired by fluid dynamics (incompressible Euler equation with specific boundary conditions).
See also \cite{papadakis:hal-00816211}. % TODO: comprendre Lagrangien augmenté et opérateur proximal.

\subsection{The discrete $\rightarrow$ discrete case}
If $\mu$ is the sum of $m$ Dirac masses and $\nu$
the sum of $n$ Dirac masses, then the problem boils down to finding the $m \times n$ coefficients $\gamma_{ij}$
that give for each pair of points $i$ of the source space and $j$ of the target space
the quantity of matter transported from $i$ to $j$. This corresponds to the transport plan
in the discrete Kantorovich problem that we have seen previously \S\ref{ch5:sec:ot:KD}.
This type of problem (referred to as an \emph{assignment problem}) can be solved by
different methods of linear programming \cite{AP:BDM:2009}. These method can be dramatically
accelerated by adding a regularization term, that can be interpreted as the entropy of
the transport plan \cite{ChristianLeonardSchroedinger}. This regularized version of 
optimal transport can be solved by highly efficient numerical algorithms \cite{DBLP:conf/nips/Cuturi13}.

\subsection{The continuous $\rightarrow$ discrete case}
This configuration, also called \emph{semi-discrete}, corresponds to a continuous function transported to a sum of Dirac masses (see the examples of c-concave functions
in \cite{gangbo1996}). This correspond to our example with factories that consume a resource, in \S\ref{sect:monge_measure}.
Semi-discrete transport has interesting connections with some notions of
computational geometry and some specific sets of convex polyhedra that were studied by Alexandrov
\cite{Alexandrov48} and later by Aurenhammer, Hoffman et Aranov \cite{DBLP:conf/compgeom/AurenhammerHA92}.
The next section is devoted to this configuration.

\section{Semi-discrete transport}
\label{chap5:sec:ot:semidiscreteOT}

\begin{figure}
      \centerline{ \includegraphics[width=50mm]{Figures/semi_discrete}}


\caption{
  Semi-discrete transport: gray levels symbolize the quantity of a resource that will
  be transported to 4 factories. Each factory will be allocated a part of the terrain
  in function of the quantity of resource that it should collect.
}
\label{fig:semi_discrete}
\end{figure}

We now suppose that the source measure $\mu$ is continuous, and that the target measure
$\nu$ is a sum of Dirac masses. A practicle example of this type of configuration corresponds
to a resource which available quantity is represented by a function $u$. The resource
is collected by a set of $n$ factories, as shown in Figure \ref{fig:semi_discrete}.
Each factory is supposed to collect a certain prescribed quantity of resource $\nu_j$. Clearly, the sum
of all prescriptions corresponds to the total quantity of available resource ($\sum_{j=1}^n \nu_j= \int_X u(x) dx$).

\subsection{Semi-discrete Monge problem}
\label{ss:semi-discrete-monge-problem}

In this specific configuration, the Monge problem becomes:

$$
  (M) \quad \quad
  \begin{array}{l}
    \inf\limits_{T:X\rightarrow Y} \int\limits_X c(x,T(x)) u(x)dx
     \mbox{\ \ subject to: } \forall j, \int_{T^{-1}(y_j)} = \nu_j.
    \end{array}
    \label{eqn:sdmonge}
$$

 A transport map $T$ associates with each point $x$ of $X$ one of the points $y_j$. Thus, it is possible
to partition $X$, by associating to each $y_j$ the region $T^{-1}(y_j)$ that contains all the points transported
towards $y_j$ by $T$. The constraint imposes that the quantity of collected resource over each region $T^{-1}(y_j)$
corresponds to the prescribed quantity $\nu_j$.

Let us now examine the form of the dual Kantorovich problem.
In terms of measure, the source measure $\mu$ is associated
with the density $u$, and the target measure $\nu$ is a sum
of Dirac masses $\nu = \sum_{j=1}^n \nu_j \delta_{y_j}$, supported by the
set of points $Y = \left\{ y_j \right\}$. We recall that in its general
form, the dual Kantorovich problem is written as follows:
\begin{equation}
   \sup\limits_{\psi \in {\bm \Psi^c}(Y)} \left[ \int_X \psi^c(x)d\mu + \int_Y \psi(y)d\nu \right].
   \label{eqn:sdk1}   
\end{equation}

In our semi-discrete case, the functional becomes a function of $n$ variables, with the following form:
\begin{eqnarray}
   \label{eqn:sdk00}
     F(\psi) = F(\psi_1, \psi_2, \ldots \psi_n) = \\
   \label{eqn:sdk2}   
     \int\limits_X \psi^c(x) u(x) dx + \sum_{j=1}^n \psi_j \nu_j = \\
   \label{eqn:sdk3}
     \int\limits_X \inf\limits_{y_j \in Y} \left[c(x,y_j) - \psi_j \right] u(x) dx + \sum_{j=1}^n \psi_j \nu_j = \\
   \label{eqn:sdk4}
     \sum\limits_{j=1}^n \quad \int\limits_{\Lag_\psi^c(y_j)} \left( c(x, y_j) - \psi_j \right) u(x) dx + \sum_{j=1}^n \psi_j \nu_j. 
\end{eqnarray}

The first step \eqref{eqn:sdk2} takes into account the nature of the measures $\mu$ and $\nu$. 
In particular, one can notice that the measure $\nu$ is completely defined by the scalars $\nu_j$ associated
with the points $y_j$, and the function $\psi$ is defined by the scalars $\psi_j$ that correspond to its value
at each point $y_j$. The integral $\int_Y \psi(y)d\nu$ becomes the dot product 
$\sum_j \psi_j \nu_j$. Thus, the functional that corresponds to the dual Kantorovich problem becomes a function $F$
that depends on $n$ variables (the $\psi_j$). Let us now replace the c-conjugate $\psi^c$ with its expression, which gives
\eqref{eqn:sdk3}. The integral in the left term can be reorganized, by grouping 
the points of $X$ for which the same point $y_j$ minimizes $c(x,y_j) - \psi_j$, which gives \eqref{eqn:sdk4},
where the \emph{Laguerre cell} $\Lag_\psi^c(y_j)$ is defined by:
$$
    \Lag_\psi^c(y_j) = \left\{ x \in X \quad | \quad c(x,y_j) - \psi_j \le c(x,y_k) - \psi_k \quad \forall k \right\}.
$$

The Laguerre diagram, formed by the union of the Laguerre cells, is a classical structure in computational geometry.
In the case of the $L_2$ cost $c(x,y) = \| x - y \|^2$, it corresponds to the \emph{power diagram},
that was studied by Aurenhammer at the end of the 80's \cite{DBLP:journals/siamcomp/Aurenhammer87}. One of its
particularities is that the boundaries of the cells are rectilinear, making it reasonably easy to design computational
algorithms to construct them.

\begin{figure}[b]
      \centerline{
         \includegraphics[width=.49\columnwidth]{Figures/concavite1}
         \includegraphics[width=.49\columnwidth]{Figures/concavite2}         
      }
      \caption{
        The objective function of the dual Kantorovich problem is concave, because its graph
        is the lower envelope of a family of affine functions.
      }
      \label{fig:concave}
\end{figure}

\subsection{Concavity of $F$}
The objective function $F$ is a particular case of the Kantorovich dual,
and naturally inherits its properties, such as its concavity (that we did
not discuss yet). This property is interesting both from a theoretical
point of view, to study the existence and uniqueness of solutions, and from
a practical point of view, to design efficient numerical solution mechanisms.
In the semi-discrete case, the concavity of $F$ is easier to prove than in
the general case. We summarize here the proof by 
Aurenhammer \emph{et. al} \cite{DBLP:conf/compgeom/AurenhammerHA92},
that leads to an efficient algorithm
\cite{DBLP:journals/cgf/Merigot11}, \cite{Levy2015}, \cite{DBLP:journals/corr/KitagawaMT16}.

\begin{theorem}
The objective function $F$ of the semi-discrete Kantorovich dual problem \eqref{eqn:sdk1} is concave.
\end{theorem}
\begin{proof}
Consider the function $G$ defined by:
\begin{equation}
G(A, \left[ \psi_1, \ldots \psi_n \right]) = \int\limits_X \left( c(x, y_{A(x)}) - \psi_{A(x)} \right) u(x) dx,
\end{equation}
and parameterized by an \emph{assignment} $A~: X \rightarrow [1 \ldots n]$, that is, a function that
associates with each point $x$ of $X$ the index $j$ of one of the points $y_j$.
If we denote $A^{-1}(j) = \left\{ x | A(x) = j \right\}$, then $G$ can be also written as:
\begin{equation}
\begin{array}{lcl}
G(A,\psi) & = & \sum\limits_j \int\limits_{A^{-1}(j)} \left(c(x,y_j) - \psi_j\right) u(x) dx \\[2mm]
          & = & \sum\limits_j \int\limits_{A^{-1}(j)} c(x,y_j) u(x) dx - \sum\limits_j \psi_j \int\limits_{A^{-1}(j)} u(x) dx.
\end{array}
\label{ch4:eqn:ot:Gaffine}
\end{equation}
The first term does not depend on the $\psi_j$, and the second one is a linear combination of the $\psi_j$ coefficients, thus,
for a given fixed assignment $A$, $\psi \mapsto G(A, \psi)$ is an affine function of $\psi$. Figure \ref{fig:concave} depicts
the appearance of the graph of $G$ for different assignment. The horizontal axis symbolizes the components of the
vector $\psi$ (of dimension $n$) and the vertical axis the value of $G(A,\psi)$. For a given assignment $A$, the graph
of $G$ is an hyperplane (symbolized here by a straight line). 

Among all the possible assignments $A$, we distinguish $A^\psi$ that associates with a point $x$ the index $j$ of the Laguerre
cell $x$ belongs to ($A$ is undefined on the set of Laguerre cell boundaries, this does not cause any difficulty
because this set has zero measure), that is:
$$
    A^\psi(x) = \arg \min\limits_j \left[ c(x,y_j) - \psi_j \right]
$$

For a fixed vector $\psi = \psi^0$, among all the possible assignments $A$,
the assignment $A^{\psi^0}$ minimizes the value $G(A, \psi^0)$,
because it minimizes the integrand pointwise (see Figure \ref{fig:concave} on the left.
Thus, since $G$ is affine with respect to $\psi$, the graph of the function $\psi \rightarrow
G(A^\psi,\psi)$ is the lower envelope of a family of hyperplanes (symbolized as straight lines
in Figure \ref{fig:concave} on the right), hence $\psi \rightarrow
G(A^\psi,\psi)$ is a concave function. Finally, the objective function $F$ of the
dual Kantorovich problem can be written as
$F(\psi) = G(A^\psi,\psi) + \sum_j \nu_j \psi_j$, that is, the sum of
a concave function and a linear function, hence it is also a concave function.
\end{proof}

\subsection{The semi-discrete optimal transport map}
Let us now examine the subdifferential $\partial_c \psi$, that is,
the set of points $(x \in X,y \in Y)$ connected by the optimal
transport map. We recall that the subdifferential $\partial_c \psi$ is defined by
$ \partial_c \psi = \left\{ x,y_j \quad | \quad \psi^c(x) + \psi_j = c(x,y_j) \right\}$. Consider
a point $x$ of $X$ that belongs to the Laguerre cell $\Lag^\psi(y_j)$. The 
subdifferential $[\partial_c\psi](x)$ at $x$ is defined by:
\begin{eqnarray}
\label{eqn:Tdisc1}
   [\partial_c\psi](x) = \left\{ y_k \quad | \quad \psi^c(x) + \psi_k = c(x,y_k) \right\} = \\
\label{eqn:Tdisc2}   
     \left\{ y_k \quad | \quad \inf_{y_l} \left[ c(x,y_l) - \psi_l \right] + \psi_k = c(x,y_k) \right\} = \\
\label{eqn:Tdisc3}     
     \left\{ y_k \quad | \quad c(x,y_j) - \psi_j + \psi_k = c(x,y_k) \right\} = \\
\label{eqn:Tdisc4}     
     \left\{ y_k \quad | \quad c(x,y_j) - \psi_j  = c(x,y_k) - \psi_k \right\} = \\
\label{eqn:Tdisc5}          
     \left\{ y_j \right\}. 
\end{eqnarray}

In the first step \eqref{eqn:Tdisc2}, we replace $\psi^c$ by its definition, then we
use the fact that $x$ belongs to the Laguerre cell of $y_j$ \eqref{eqn:Tdisc3}, and finally,
the only point of $Y$ that satisfies \eqref{eqn:Tdisc4} is $y_j$ because we have supposed 
$x$ \emph{inside} the Laguerre cell of $y_j$.

To summarize, the optimal transport map $T$ moves each point $x$ to the point $y_j$
associated with the Laguerre cell $\Lag^\psi(y_j)$ that contains $x$. The vector $\psi_1, \ldots \psi_n$
is the unique vector that maximizes the discrete dual Kantorovich function $F$ such that $\psi$ is c-concave.
It is possible to show that $\psi$ is c-concave if and only if no Laguerre cell is empty of matter, that is the integral
of $u$ is non-zero on each Laguerre cell.

We now proceed to compute the first and second order derivatives of the objective function $F$. These derivatives
are useful in practice to design computational algorithms.

\begin{figure}[t]
      \centerline{
         \includegraphics[width=6cm]{Figures/enveloppe}
      }
      \caption{
        Illustration of the \entreguillemets{envelope theorem}: at a given point,
        the tangent to the lower envelope of a family of functions corresponds to the tangent to the
        lowest function of the family at that point.
      }
      \label{fig:enveloppe}
\end{figure}

\subsection{First-order derivatives of the objective function}

Since it is concave,  $F$ admits a unique maximum $\psi^*$, characterized by $\nabla
F(\psi^*) = 0$ where $\nabla F$ denotes the gradient. Let us now examine the form of the gradient $\nabla F
= \nabla \left(G(A\psi, \psi) + \sum_{j=1}^n \psi_j \nu_j\right)$.  At
first sight, the term in $G$ will be delicate to compute, because it has integrals
over domains $A^{-1}(y_j)$ that vary in function of the vector of parameters $\psi$. In general
differentiaing integrals over varying domains requires using Reynold's formula (see for instance \cite{EmmanuelMaitreCours}
P. 7). However, in this particular case, recalling that the left term corresponds to the lower envelope of
a family of function, one can use the \emph{envelope theorem} \cite{RePEc:ecm:emetrp:v:70:y:2002:i:2:p:583-601}, illustrated
in Figure \ref{fig:enveloppe}: at a given point, the gradient of the lower envelope of a family of functions
corresponds to the gradient of the function with the smallest value at that point. Thus, in this particular case,
one can ignore the variations of the Laguerre cells with respect to $\psi$. By replacing $G(A_\psi, \psi)$ with its
expression \eqref{ch4:eqn:ot:Gaffine}, one obtains:
$$
\begin{array}{lcl}
\frac{\partial}{\partial \psi_j} G(A_\psi, \psi) & = & \\[2mm] \frac{\partial}{\partial \psi_j} \left(
    \sum\limits_k \int\limits_{A^{-1}_\psi(k)} c(x,y_k) u(x) dx - \sum\limits_k \psi_k \int\limits_{A^{-1}_\psi(k)} u(x) dx \right) & = & \\[6mm]
                                               -\int\limits_{A^{-1}_\psi(j)} u(x) dx & = &\\[6mm]
                                               -\int\limits_{\Lag^\psi(y_j)} u(x) dx.
\end{array}
$$
Recalling that the objective function $F$ is given by $F(\psi) = G(A_\psi, \psi) + \sum_j \nu_j \psi_j$, we finally obtain:
\begin{equation}
\frac{\partial F}{\partial \psi_j} = \nu_j - \int\limits_{\Lag^\psi(y_j)} u(x) dx.
\label{eqn:grad}
\end{equation}

Since the objective function $F$ is concave, it admits a unique maximum.
At this maximum, all the components of the gradient vanish, this implies
that the quantity of matter obtained at $y_j$, that is the integral of
$u$ on the Laguerre cell of $y_j$, corresponds to the prescribed quantity
of matter, that is $\nu_j$.

\subsection{Second order derivatives of the objective function}
The coefficients of the Hessian matrix $\partial^2 F / \partial \psi_i \partial \psi_j$
are slightly more difficult to compute, since here, we cannot invoke the envelope theorem.
We cannot avoid invoking Reynold's formula. We do not detail the computations for length
considerations, but give the final result.

In the particular case of the $L_2$ cost, that is $c(x,y) = \| x - y \|^2$, the second-order derivatives are given by:
\begin{equation}
\begin{array}{lcl}
  \displaystyle\frac{\partial^2 F}{\partial \psi_i \partial \psi_j}  & = &
       \displaystyle\frac{\int\limits_{\Lag^\psi(y_i) \cap \Lag^\psi(y_j)} u(x) dx}{\| y_j - y_i \|^2}
       \quad \forall j \neq i \\[6mm]
  \displaystyle\frac{\partial^2 F}{\partial \psi_i \partial \psi_i} & = &
    \displaystyle-\sum\limits_{j \neq i} \frac{\partial^2 F}{\partial \psi_i \partial \psi_j}
\end{array}       
\label{eqn:hess}    
\end{equation}

\subsection{A computational algorithm for $L_2$ semi-discrete optimal transport}
\label{sect:algo}
With the definition of $F(\psi)$, the expression of its first order derivatives
(gradient $\nabla F$) and second order derivatives (Hessian matrix $\nabla^2 F = \left(\partial^2 F / \partial \psi_i \partial \psi_j\right)_{ij}$),
we are now equipped to design a numerical solution mechanism that computes semi-discrete optimal transport by
maximizing $F$, based on a particular version \cite{DBLP:journals/corr/KitagawaMT16} of Newton's optimization method \cite{Nocedal2006NO}:
$$
\begin{array}{ll}
   \mbox{\textbf Input:}  & \mbox{a mesh that supports the source density } u \\
                        & \mbox{the points } (y_j)_{j=1}^n \\
                        & \mbox{the prescribed quantities } (\nu_j)_{j=1}^n \\[2mm]
   \mbox{\textbf Output:} & \mbox{the (unique) Laguerre diagramm } \Lag^\psi \mbox{ such that:}\\
                        & \quad \int\limits_{\Lag^\psi}(y_j) u(x)dx = \nu_j \quad \forall j\\[5mm]
                        \hline \\
   (1) & \psi \leftarrow [ 0 \ldots 0 ] \\
   (2) & \mbox{While convergence is not reached} \\
   (3) & \quad \mbox{Compute } \nabla F \mbox{ and } \nabla^2 F \\
   (4) & \quad \mbox{ Find } p \in {\mathbb R}^n \mbox{ such that } \nabla^2 F(\psi) p = -\nabla F(\psi) \\
   (5) & \quad \mbox{ Find the descent parameter } \alpha \\
   (6) & \quad \psi \leftarrow \psi + \alpha p \\
   (7) & \mbox{End while}
\end{array}
$$

The source measure is given by its density, that is a positive piecewise linear function $u$,
supported by a triangulated mesh (2D) or tetrahedral mesh (3D) of a domain $X$. The target measure
is discrete, and supported by the pointset $Y = (y_j)_{j=1}^n$. Each target point will receive the prescribed
quantity of matter $\nu_j$. Clearly, the prescriptions should be balanced with the available resource,
that is $\int_X u(x)dx = \sum_j \nu_j$. The algorithm computes for each point of the target measure
the subset of $X$ that is affected to it through the optimal transport, $T^{-1}(y_j) = \Lag^\psi(y_j)$,
that corresponds to the Laguerre cell of $y_j$. The Laguerre diagram is completely determined
by the vector $\psi$ that maximizes $F$.

Line (2) needs a criterion for convergence. The classical convergence
critertion for a Newton algorithm uses the norm of the gradient of
$F$.  In our case, the components of the gradient of $F$ have a
geometric meaning, since $\partial F / \partial \psi_j$ corresponds to
the difference between the prescribed quantity $\nu_j$ associated with
$j$ and the quantity of matter present in the Laguerre cell of $y_j$
given by $\int_{\Lag^\psi(y_j) u(x)dx}$. Thus, we can decide to stop
the algorithm as sonn as the largest absolute value of a component
becomes smaller than a certain percentage of the smallest
prescription. Thus we consider that convergence is reached if $\max |
\nabla F_j | < \epsilon \min_j \nu_j$, for a user-defined $\epsilon$
(typically 1\% in the examples below).

Line (3) computes the coefficients of the gradient and the Hessian matrix of $F$, using \eqref{eqn:grad} and
\eqref{eqn:hess}. These computations involve integrals over the Laguerre cells and over their boundaries.
For the $L_2$ cost $c(x,y) = \| x - y \|^2$, the boundaries of the Laguerre cells are rectilinear, which
dramatically simplifies the computations of the Hessian coefficients \eqref{eqn:hess}. In addition,
it makes it possible to use efficient algorithms to compute the Laguerre diagram
\cite{DBLP:journals/cj/Bowyer81,journals/cj/Watson81}. Their implementation is available in several
programming libraries, such as GEOGRAM(\url{https://brunolevy.github.io/geogram/}) et
CGAL (\url{http://www.cgal.org}). Then one needs to compute the intersection between each Laguerre cell
and the mesh that supports the density $u$. This can be done with specialized algorithms~\cite{Levy2015},
also available in GEOGRAM.

Line (4) finds the Newton step $p$ by solving a linear system. We use the Conjugate Gradient algorithm \cite{Stiefel1952}
with the Jacobi preconditioner. In our empirical experiments below, we stopped the conjugate iterations as soon
as $\| \nabla^2F p + \nabla F \| / \| \nabla F \| < 10^{-3}$.


Line (5) determines the descent parameter $\alpha$. A result due to Mérigot and Kitagawa \cite{DBLP:journals/corr/KitagawaMT16}
ensures the convergence of the Newton algorithm if the measure of the smallest Laguerre cell remains larger than a certain
threshold (that is, half the smallest prescription $\nu_j$). There is also a condition on the norm of the gradient $\|\nabla F\|$
that we do not repeat here (the reader is referred to Mérigot and Kitagawa's original article for more details). In our implementation,
starting with $\alpha = 1$, we iteratively divide $\alpha$ by two until
both conditions are satisfied.

Let us now make one step backwards and think about the original
definition of Monge's problem \eqref{eqn:monge}. We wish to stress
that the initial constraint (local mass preservation) that
characterizes transport maps was terribly difficult. It is remarkable
that after several rewrites (Kantorovich relaxation, duality,
c-convexity), the final problem becomes as simple as optimizing a
regular ($C^2$) concave function, for which computing the gradient and
Hessian is easy in the semi-discrete case and boils down to evaluating
volumes and areas in a Laguerre diagram. We wish also to stress that
the computational algorithm did not require to make any approximation
or discretization. The discrete, computer version is a particular
setting of the general theory, that fully describes not only transport
between smooth objects (functions), but also transport between less
regular objects, such as pointsets and triangulated meshes.  This is
made possible by the rich mathematical vocabulary (measures) on which
optimal transport theory acts. Thus, the computational algorithm is
an elegant, direct verbatim translation of the theory into a computer program.

\section{Fluid simulation based on optimal transport}
\label{sect:fluidot}

Based on the very general mathematical setting presented in previous sections, we now provide the link between Optimal Transport and fluid simulation, starting from classical discretization and then deriving this new representation. The reader interested in general fluid simulation is referred to Eurographics state-of-the-art report~\cite{Ihmsen14, Koschier2022}. We give here a short introduction to the subset, then focus on Optimal Transport-based methods.

The discretization of the motion of a fluid is classically performed from two different points of view:\begin{itemize}
    \item \emph{Eulerian}: the state of the fluid is represented by a set of fields, that map fixed locations in space to physical quantities (density, velocity). In terms of implementation, such fields are sampled in regular grids (that can be optionally adaptively subdivided). Impressive results were notably shown~\cite{Li2023, Li2024, Wang2025} using the \emph{Lattice-Boltzmann Method} (LBM)~\cite{McNamara1988};
    \item \emph{Lagrangian}: the physical quantities are attached to a representation that follows the motion of the fluid, such as discrete particles. To ease evaluating the differential quantities involved in the physics, it is possible to attach to each particle a function with compact support. Notably, \emph{Smoothed Particles Hydrodynamics} (SPH)~\cite{Gingold1977} is an efficient numerical scheme which is very popular in Computer Graphics~\cite{Ihmsen14, Koschier2022}.
\end{itemize} Despite their more accurate tracking of the matter, \emph{Lagrangian} discretizations often face two challenges: to robustly discretize differential operators and to accurately conserve volume across the numerical simulation. In particular, it can be difficult to precisely evaluate differential quantities relying on a meshless representation like SPH. On the other hand, other Lagrangian discretizations are based on an explicit mesh which directly provides an accurate definition of the neighborhood but also need complex topological edition across the simulation. Furthermore, ensuring precisely that the volume of matter linked to each particle remains constant (\textit{e.g.} to simulate an incompressible fluid) remains a challenge since it requires to accurately pair a subset of the matter to each particle. Several schemes called incompressible SPH~\cite{Ihmsen2014, Bender2017} were proposed to address this challenge indirectly, by updating the velocity field in such a way that the volume conservation constraint is satisfied. In this part, we describe an alternative scheme, ensuring volume conservation \emph{by design}.

Methods based on Optimal Transport take a different approach, still with a Lagrangian representation of the fluid. The difference is that instead of using smoothed particles, the Lagrangian representation is based on a set of $n$ cells $\left(\Lag^\psi(y_i)\right)_{i=1}^n$ parameterized by the coordinates $(y_i)_{i=1}^n$ with $y_i \in \mathbb R^d$, and the volumes $(\nu_i)_{i=1}^n \in \mathbb R^n$. Based on the results presented in the previous sections, it is possible to construct such a representation from a Laguerre diagram  \[
    \Lag^\psi(y_i) \coloneqq \left\{ x \mid \|y_i - x\|^2 - \psi_i \leq \|y_j - x\|^2 - \psi_j, \; i \neq j \right\},
\] where $(y_i)_{i=1}^n$ are the seeds of the Laguerre diagrams, and 
$\Psi \in \mathbb R^n$ the vector of weights defining the Laguerre diagram 
(see \Cref{ss:semi-discrete-monge-problem}). Now one needs to deduce the weight vector $\Psi$ from our parameters, that is, the $n$ seeds locations $(y_i)$ and the associated prescribed volumes $(\nu_i)$. We know, thanks to the previous developments, that the weight vector that satisfies the prescribed volumes $\nu_i = |\Lag^\psi(y_i)|\ \forall i$ exists and is unique. We also know that from a practical point of view, the weight vector can be determined by a Newton algorithm (\Cref{sect:algo}). To simulate a fluid with such a representation, one now needs to deduce from the model of the fluid (Navier Stokes) the laws of evolution of the parameters (the seed positions $y_i$). 

This family of methods were inspired by ideas from Brenier and Benamou, leveraging connections between optimal transport and fluid dynamics, presented in their seminal works~\cite{Brenier1989, Benamou2000}. In Computer Graphics, these have notably paved the way for the first fluid simulations based on a semi-discrete Optimal Transport by de~Goes~et~al.~\cite{Goes2015}. In their work, the evolution of the mesh is indirectly computed by advecting the seeds that define the Laguerre diagram. To ensure a uniform distribution of the seeds, they are then explicitly relocated to the centroids of their cells at each time step. An alternative discretization was proposed by Gallouet \& Merigot~\cite{Gallouet2017}, that provably converges to the solution of the continuous problem for incompressible Euler fluids. In a certain sense, the Gallouet and Merigot scheme implements Arnold's abstract point of view that
the fluid's configuration follows a geodesic on the manifold of incompressible motion. The fluid's trajectory is projected onto the manifold by a "spring" that connects each seed to its centroid. Inspired by these two works, multiple numerical schemes were proposed, enabling massively parallel implementations on GPU with adaptive splitting and merging of cells~\cite{Zhai2020}, as well as hybrid algorithms combining Lagrangian and Eulerian representations~\cite{Qu2022, Qu2023}.

Using a Laguerre diagram to represent the fluid is a natural way of constructing a Lagrangian mesh representation while strictly enforcing volume conservation, 
but it also comes with some constraints. Laguerre diagram, as an unbounded representation, covers the entire domain which is incompatible with free surface simulations. To address this issue, previous works either explicitly represent the air medium with dedicated cells~\cite{Goes2015, Gallouet2017, Zhai2020}, or use a discrete representation which can be easily restricted~\cite{Qu2022, Qu2023}. The next sections explain how the notion of Partial Optimal Transport can be used to formalize this type of configuration with free surfaces, and how it can be approximated for computer implementation. 

\section{Free boundaries through partial optimal transport}
%----------------------------------------------------------
\label{sect:SDPOT}
 
Let us now consider a different setting, where the simulated object does not fill the volume $\Omega$ entirely. That is, the sum of the prescribed volumes $\nu_i$ is smaller than the volume of $\Omega$. This problem is referred to as a \emph{partial optimal transport} problem. We show how
partial optimal transport can be considered as a particular instance of optimal transport. To do so, We still consider the optimal transport problem (Eq. \ref{eqn:Monge}) and its Kantorovich dual (Eq. \ref{eqn:sdk1}), with the difference that this time transport is done towards a set of "objects" (pointsets) $\mathcal{O}_i$:

\begin{equation}
    \begin{array}{l}
    \inf_I \left[ \int_{\Omega} d^2(\bx,\mathcal{O}_{I(\bx)})d\bx\right] \\[2mm]
    \mbox{subject to } | I^{-1}(\mathcal{O}_i) | = \nu_i \quad \forall i
    \end{array}
\end{equation}
where:
\begin{itemize}
\item the index map $I: \Omega \rightarrow [1..N]$ assigns an object $\mathcal{O}_i$ to each point $\bx$ of $\Omega$;
\item each $\mathcal{O}_i$ is a set of points of $\Omega$ (that can contain either a single point, or an integer number of points, or be a continuous subset of $\Omega$);
\item each $\mathcal{O}_i$ is supposed to receive a prescribed quantity of matter $\nu_i$ (constraint);
\item the distance between a point $\bx$ and an object $\mathcal{O}$ is defined by $d(\bx, \mathcal{O}) = \inf_{\by \in \mathcal{O}} \left[ \| \bx - \by \| \right]$;
\item $I^{-1}(\mathcal{O}_i) = \{ \bx \ | \ I(\bx) = i \}$;
\item Given an index map $I$, the corresponding map $T_I: \Omega \rightarrow \Omega$ is defined by $T_I(\bx) = \mbox{arg inf}_{\by \in \mathcal{O}_{I(\bx)}} \| \bx - \by \|^2$.
\end{itemize}

Consider now the dual problem:
$$
\begin{array}{l}  
   \sup_\psi \left[
   K(\psi) = \int\limits_{\Omega} \inf\limits_i \left[ d^2(\bx, \mathcal{O}_i) - \psi_i \right] d\bx + \sum\limits_i \psi_i \nu_i \right] \\[2mm]
   \mbox{subject to } \psi^{cc} = \psi
\end{array}.   
$$

\begin{figure}[t]
  \centerline{
    \includegraphics[width=\columnwidth]{Figures/POT3}
  }
  \caption{Left: approximation of partial optimal transport using a     
     discretization of the background ("ghost cells"). 
     Right: exact solution to partial optimal transport, where the background is considered as a continuum. Each cell $V_i^\psi$ is the intersection between the Laguerre cell $Lag_i^\Psi$ and the ball centered on $\bx_i$ of radius $\sqrt{\psi_i}$.
  }
  \label{fig:POT}
\end{figure}

It keeps the structure of our initial Kantorovich dual (Eq. \ref{eqn:sdk1}, \ref{eqn:sdk4}), and the argument for the convexity, smoothness, existence and uniqueness of the solution still hold, but now we have a more general setting that we can use to account for free boundaries: we consider now a set of $N+1$ objects $\mathcal{O}_{i=0}^N$:
\begin{itemize}
\item The first object $\mathcal{O}_0$ is a set of $M$ points $\{\by_1, \ldots \by_M\}$ uniformly distributed in $\Omega$;
\item each $\mathcal{O}_i$ for $1 \le i \le N$ is the singleton $\{\bx_i\}$; 
\item the object $\mathcal{O}_0$ is associated with $\nu_0 = |\Omega| - \sum_{i=1}^N \nu_i$, that is all the volume of $\Omega$ not affected to the points $\{\bx_i\}$.
\end{itemize}

This configuration is depicted in Figure \ref{fig:POT}-Left: the cells
of the fluid are displayed in grey. Each of them is a convex polygon. The background points $\by_{i=1}^M$ fill the entire domain $\Omega$, including the part occupied by the fluid. All the $\by_i$ points  share the same Lagrange multiplier $\psi_0$. This may be interpreted as the fact that it takes no cost for background cells to exchange matter (air). For a part occupied by the fluid, each cell $Lag^\psi_i$ that correspond to the fluid has a Lagrange multiplier $\psi_i$ larger than $\psi_0$, then the fluid cell "shadows" the background cells. On the boundary of the fluid, the presence of the background cells limits the size of the fluid cells. The external boundary of the fluid is formed by straight segment shared by a fluid cell and a background cell (closeup in Figure \ref{fig:POT}-Left).  

A similar technique was used in \cite{DBLP:journals/tog/GoesWHPD15}, where it is proposed to insert points $\by_i$ in the vicinity of the boundary of the fluid (in the cited article, the corresponding cells are called "ghost cells"). The main advantage it that it allows directly reusing a discrete optimal transport implementation to solve transport problems with free boundaries. However, the "ghost cells" technique has two drawbacks: 
\begin{itemize}
    \item first, it is difficult to know \emph{in advance} where to insert the ghost cells to make sure the boundary of the fluid is accurately represented, that is, there is no easy way of predetermining which ghost cells influence the fluid;
    \item second, the ghost cells have a significant computational cost: computing a Laguerre diagram costs $O(N\sqrt[d]{N})$ (where $d \in \{2,3\}$ is the dimension), and the Newton algorithm converges in $O(N \log(N))$ iterations (empirical results in \cite{levy:hal-03081581}). The "ghost cells" techniques computes a diagram with $N+M$ vertices (instead of $N$), where $M$ needs to be sufficiently large to accurately capture the free boundary.  
\end{itemize}

\begin{figure}[t]
  \centerline{
    \includegraphics[width=.45\columnwidth]{Figures/lag4}\hfill
    \includegraphics[width=.45\columnwidth]{Figures/lag3}
  }
  \caption{Left: The minimization diagram of the familly of functions
    $f_i(\bx) = \| \bx - \bx_i \|^2 - \psi_i$ augmented with the graph of the function
    $f(\bx) = 0$ (horizontal plane). Right: seen from below, the obtained cells
    $V_i^\psi$ are  the intersections between the
    Laguerre cells $Lag^\bpsi_i$ and the disks of radii $\sqrt{\psi_i}$ centered
    on the $\bx_i$'s (or the empty set if $\psi_i$ is negative).
  }
  \label{fig:laguerre2}
\end{figure}


To overcome these limitations, we propose an alternative technique: let us now imagine that the number of background points tends to infinity. More simply put, the background object $\mathcal{O}_0$ becomes the entire $\Omega$ domain (and it is still associated with $\nu_0 = |\Omega| - \sum_{i=1}^N \nu_i$). Let us also remember that the vector $\psi_{i=0}^N$ is independent on a translation, thus, w.l.o.g. it is possible to choose $\psi_0 = 0$. Let us now consider a cell $V_i^\psi$ associated with a fluid particle:
$$
\begin{array}{rl}
   V_i^\psi =& \{ \bx | d^2(\bx, \mathcal{O}_i) - \psi_i < d^2(\bx, \mathcal{O}_j) - \psi_j) \; \forall j \neq i \} \\[3mm]
            =& \{ \bx | \| \bx - \bx_i \|^2 - \psi_i < \| \bx - \bx_j \|^2 - \psi_j) \; \forall 1 \le j \neq i \le N \\
             & \mbox{ and} \; \| \bx - \bx_i \|^2 - \psi_i < d^2(\bx, \mathcal{O}_0) - \psi_0 \} \\ [3mm]
            =& \{\bx | \| \bx - \bx_i \|^2 - \psi_i < \| \bx - \bx_j \|^2 - \psi_j) \; \forall 1 \le j \neq i \le N \\
             & \mbox{and} \; \| \bx - \bx_i \|^2  <  \psi_i \}. 
\end{array}
$$

The last line is obtained by remembering that $\psi_0 = 0$, and by noticing 
that $d^2(\bx, \mathcal{O}_0) = 0$ for all $\bx$ in $\Omega$ since $\mathcal{O}_0 = \Omega$. In other words, $V_i^\psi$ is the intersection between the Laguerre cell 
$Lag^\psi_i$ and a disk of radius $\sqrt{\psi_i}$ centered on $\bx_i$, as shown in Figure \ref{fig:POT}-Right, or the empty set if $\psi_i$ is negative. Note that a negative $\psi_i$ cannot happen in our context, since it would contradict the convexity of $\psi$ ($\psi^{cc} = \psi \Leftrightarrow V_i^\psi \neq \emptyset \ \forall 1 \le i \le N$), guaranteed by all the iterations of the KMT Newton algorithm \cite{DBLP:journals/corr/KitagawaMT16} (more on this below). 
From this observation, not only "ghost cells"
are no longer needed, but also we can much more accurately compute the boundary of the fluid, by computing the intersection between the Laguerre cells and a set of disks. Put differently, what is computed is the limit case where the number of ghost cells $M$ tends to infinity. Not only the result will be more precise, but the overall computational cost will be significantly reduced, since the number of vertices in the Laguerre diagram is not increased.


As shown in Figure \ref{fig:laguerre2}, to gain more intuition about this setting, it is also possible to take the "minimization diagram" point of view: the distance to $\mathcal{O}_0 = \Omega$
is equal to zero on $\Omega$, and $\psi_0=0$, thus the graph of the associated function $f_0$ is the horizontal plane $Z=0$. The other functions $f_i$
are shifted paraboloids just like before. Now, looking at the diagram from above, one will see the tip of the paraboloids intersected by the horizontal planes
(disks). The paraboloids can also intersect each other, forming polygonal cells with straight edges (that correspond to projected parabolas), just like in the previous configuration. 

\section{Semi-Discrete Partial Optimal Transport}
\label{sect:pot-algo}

We shall now explain how to design a numerical solution mechanism. The associated algorithm takes the following inputs and produces the following outputs:
$$
\begin{array}{ll}
   \mbox{\textbf Input:}    &- \mbox{ The domain } \Omega\ (\mbox{ can be } [0,1]^d \mbox{ or a simplicial mesh})\\
                        &- \mbox{ a set of } N \mbox{ points } \bx_i \in \Omega \\
                        &- \mbox{ prescribed volumes } (\nu_i)_{i=1}^N \mbox{ such that } \sum_i \nu_i \le |\Omega| \\[2mm]
   \mbox{\textbf Output:} &- \mbox{the (unique) vector } \psi^* \in \mathbb{R}^N \mbox{ that maximizes } K(.) \\
                      &- \mbox{the cells } (V_i^{\psi^*})_{i=1}^N
                         \mbox{ such that } |V_i^{\psi^*}| = \nu_i \\
\end{array}   
$$

The core of the algorithm solves the following optimization problem:
\begin{equation}
\begin{array}{l}
    \sup_\psi \left[ K(\psi) = 
       \sum\limits_{i=1}^N\ \int\limits_{V_i^\psi} (\| \bx - \bx_i \|^2 - \psi_i) d\bx + \sum\limits_{i=1}^N \nu_i \psi_i
    \right] \\[4mm]
    \mbox{where } V_i^\psi = Lag_i^\psi \cap \{ \bx \ | \ \| \bx - \bx_i \|^2 \le \psi_i \} \\[2mm]
    \mbox{subject to } V_i^\psi \neq \emptyset \quad \forall 1 \le i \le N
\end{array}
\label{eqn:KPOT}
\end{equation}

The KMT algorithm due to Kitagawa, Merigot and Thibert \cite{DBLP:journals/corr/KitagawaMT16} is a Newton algorithm that produces a series of iterates $\psi^{(k)}$ that provably converges to the solution of Eq. \ref{eqn:KPOT}. Each iterate satisfies the convexity constraint $\psi^{cc} = \psi$ (that is, $V_i^\psi \neq \emptyset \quad \forall k, \forall 1 \le i \le N$, see also \cite{Levy2018}). Adapting the KMT algorithm to our specific context has numerical and geometrical aspects detailed below.

\subsection{Numerical aspects}

The KMT algorithm follows the classical structure of a Newton algorithm, already seen in \Cref{sect:algo} but updated to this Partial Optimal Transport scheme:
$$
\begin{array}{ll}
   \ (1): & \psi \leftarrow [ 0 \ldots 0 ] \\
   \ (2): & \mbox{Loop} \\
   \ (3): & \quad \mbox{Compute the cells } (V_i)_{i=1}^N \\
   \ (4): & \quad \mbox{Compute the gradient } \nabla K(\psi) \\
   \ (5): & \quad \mbox{If } \| \nabla K(\psi) \|_\infty < \epsilon_K \mbox{ then Exit loop} \\
   \ (6): & \quad \mbox{Compute the Hessian matrix } \nabla^2 K (\psi) \\
   \ (7): & \quad \mbox{Solve for } \bp \in {\mathbb R}^n \mbox{ in  } \nabla^2 K(\psi) \bp = -\nabla K(\psi) \\
   \ (8): & \quad \mbox{Find the descent parameter } \alpha \\
   \ (9): & \quad \psi \leftarrow \psi + \alpha \bp \\
   \ (10): & \mbox{End loop}
\end{array}
$$

The algorithm above needs to evaluate the gradient and Hessian matrix of $K(.)$ at each main loop iteration. The coefficients of the gradient and Hessian matrix can be deduced from the cells $(V_i^\psi)_{i=1}^N$ that are computed at step (3). The algorithm that computes the cells is detailed later in the next subsection on the geometric aspects. 
The components ${\partial K}/{\partial \psi_i}$ of the gradient $\nabla K(.)$ are given by the following expression \cite{DBLP:journals/corr/KitagawaMT16, Levy2018}:
\begin{equation}
  \frac{\partial K}{\partial \psi_i}  =  \nu_i - |{V_i^\psi}|. 
\end{equation}
In other words, each component of the gradient corresponds to the prescribed volume $\nu_i$ associated with a point $\bx_i$ minus the volume of the cell $V_i^\psi$. For the vector $\psi^*$ that maximizes $K(.)$, all components of the gradient vanish, which means that each cell $V^{\psi^*}_i$ has exactly the prescribed volume $\nu_i$.  

This expression of the gradient leads also to a natural stopping criterion (line 5), the largest component of the gradient corresponds to the maximum volume error. We stop the algorithm as soon as it is smaller than a prescribed $\epsilon_K$ (typically one percent of $\nu_i$).

We now consider the Hessian matrix computed at step (6). Adapting the formulas in \cite{DBLP:journals/corr/KitagawaMT16, Levy2018} to our context, the coefficients of the Hessian are given by:
\begin{eqnarray}
  \nonumber
  \frac{\partial^2 K}{\partial \psi_i \partial \psi_j} & = & \frac{1}{2} \frac{|V^\psi_{ij}|}{\| \bx_j - \bx_i \|} \quad \mbox{if } i \neq j \\
  \frac{\partial^2 K}{\partial \psi_i^2} & = & - \left(\sum_{j \neq i} \frac{\partial^2 K}{\partial \psi_i \partial \psi_j}\right)
                                               - \frac{1}{2}\frac{|V^\psi_{i0}|}{\sqrt{\psi_i}}
\end{eqnarray}
where $|V_{ij}^\psi| = |\partial V_i^\psi \cap \partial V_j^\psi|$ denotes the area of the intersection between the border of the cell
$V_i^\psi$ and the border of the cell $V_j^\psi$, that is, a polygonal facet in the Laguerre diagram clipped by the two balls $B(\bx_i,\sqrt{\psi_i})$ and $B(\bx_j,\sqrt{\psi_j})$. The term $|V^\psi_{i0}|$ corresponds to the free surface area in $\partial V_i^\psi$, that is the portion of the border of $V_i^\psi$ included in the sphere $S(\bx_i, \sqrt{\psi_i})$, that touches the unoccupied portion of $\Omega$. 

Step (7) of the algorithm computes the Newton step vector $\bp$, by solving a linear system. Except the boundary term in $\partial^2{K}/\partial \psi_i^2$, this linear system is identical as the one solved in \cite{levy:hal-03081581}: it corresponds to a Poisson equation discretized with finite elements. The matrix is sparse, with a non-zero coefficient at coefficient $(i,j)$ if and only if the cells $V^\psi_i$ and $V^\psi_j$ touch each other. The additional term $(1/2) \left( | V_{i0}^\psi | / \sqrt{\psi_i} \right)$  in $\partial^2 K / \partial \psi_i^2$ does not change the sparsity pattern of the Hessian: since the Lagrange multiplier associated with $\psi_0$ is fixed and equal to 0, there is no partial derivative with respect to it except the diagonal term. The same linear solver as in \cite{levy:hal-03081581} can be used (Jacobi-preconditioned conjugate gradient \cite{Hestenes1952} with sparse matrix stored in CRS format, with optional GPU acceleration), as well as the same stopping criterion $\| H \bp - \bg \| / \| \bg \| \le 10^{-3}$.

Once the step vector $\bp$ is computed, a good descent parameter $\alpha$ needs to be found (step (8)). In the provably convergent KMT algorithm \cite{DBLP:journals/corr/KitagawaMT16}, the descent parameter $\alpha$ is determined as follows:
$$
\begin{array}{ll}
   \ (1): & \alpha \leftarrow 1 \\
   \ (2): & \mbox{Loop} \\
   \ (3): & \quad \mbox{If  } \inf_i | V^{\psi + \alpha \bp}_i | >a_0 \\
   \ (4): & \quad \mbox{and } \| \nabla K(\psi + \alpha \bp) \| \le (1 - \alpha / 2) \| \nabla K(\psi) \| \\
   \ (5): & \quad \quad \mbox{then Exit loop} \\
   \ (6): & \quad \alpha \leftarrow \alpha / 2 \\
   \ (7): & \quad \mbox{Compute the cells } (V^{\psi+\alpha \bp}_i)_{i=1}^N \\
   \ (8): & \mbox{End loop}
\end{array}
$$
where $a_0 = \frac{1}{2}\min\left( \inf_i  |V^{\psi^(0)}_i|  , \inf_i(\nu_i) \right)$.

The KMT algorithm iteratively halves the descent parameter $\alpha$ until two criteria are met: the volume of the smallest cell needs to be larger than a threshold $a_0$ (line 3), and the norm of the gradient needs to decrease sufficiently (line 4). The threshold $a_0$ for the minimum cell volume corresponds to (half) the minimum cell volume for $\psi = 0$ (also called Voronoi diagram) and minimum prescribed area $\nu_i$.

Equipped with the KMT algorithm above, one can compute the descent parameter $\alpha$, by plugging the algorithm above into line (8) of the Newton algorithm at the beginning of this section.


\subsection{Geometrical aspects}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{Figures/primal_dual.png}
    \caption{The Bowyer-Watson algorithm computes the Laguerre diagram
    (cell boundaries displayed in white and light grey). It uses internally
    a dual representation (triangulation in black). One also needs to
    compute the intersection with $\Omega$ (the square in this
    example).} \label{fig:triangulation}
\end{figure}

To evaluate the components of the gradient and coefficients of the Hessian of $K(.)$, the KMT algorithm needs to construct the cells $V_i^\psi = Lag_i^\psi \cap B(\bx_i,\sqrt{\psi_i})$, that is, the intersection between Laguerre cells and balls centered on the $\bx_i$'s for each $\psi^{(k)}$ iterate. To do so, we first compute the Laguerre diagram, using the classical Bowyer-Watson algorithm \cite{DBLP:journals/cj/Bowyer81,journals/cj/Watson81},
see Figure \ref{fig:triangulation}. A readily-available open-source implementation is available in our GEOGRAM library \cite{WEB:GEOGRAM} (and also in the CGAL library \cite{DBLP:journals/comgeo/BoissonnatDPTY02} for programmers who have a taste for C++ templates). The details of our highly optimized implementation in GEOGRAM are given in \cite{levy:hal-03081581}.  

Given the points $\bx_i$ and the vector $\psi$, the Bowyer-Watson algorithm computes the Laguerre diagram. Internally, it uses the dual triangulation, formed by the black triangles in Figure \ref{fig:triangulation} (in 3D, they become tetrahedra). Each vertex of the Laguerre diagram is deduced from one of the triangles (resp. tetrahedra), and the Laguerre cells are obtained by traversing the triangles (tetrahedra) incident to a given $\bx_i$. Note that the Boywer-Watson algorithm computes the Laguerre diagram in the entire $\mathbb{R}^d$ domain, with infinite Laguerre cells (grey straight lines in the Figure). To obtain the $V_i$ cells, one needs to compute intersections with the domain $\Omega$ (and also with the balls centered on the $\bx_i$'s):
$$
   V_i^\psi = Lag_i^\psi \cap \Omega \cap B(\bx_i,\sqrt{\psi_i})
$$
 
In 2D, constructing the cells $V_i^\psi$ means computing intersections between convex polygons and disks, which does not represent too much difficulty (see Figure \ref{fig:POT}). However, it becomes a challenging task in 3D:
as can be seen, a large number of intersections needs to be computed, and 3D mesh intersection is known to be a delicate operation, subject to numerical precision problems: the intersection algorithm depend on a small set of functions, called \emph{geometric predicates}, that take as input some points, and that returns a discrete value. For instance, such a geometric predicate can tell whether a point $\bp_1$ is above or below another point $\bp_2$. These geometric predicates are often polynomials in the coordinates of the points. The limited precision of floating point numbers can have catastrophic consequences, for instance when the algorithm estimates that a certain point $\bp_1$ is above another point $\bp_2$, and later 
inconsistently estimates that $\bp_2$ is above $\bp_1$. Such inconsistent behavior can result in invalid combinatorics (see e.g. \cite{DBLP:conf/compgeom/Shewchuk96}).

It is important to stress that a \emph{huge} number of intersections (billions) needs to be computed, hence these numerical issues that correspond to what may be initially thought of as very unlikely corner cases occur in fact \emph{thousands times} in a typical simulation. The goal now is to find a strategy to overcome this robustness issue while keeping reasonable performance. 
The strategy is to express the problem in terms of a small number of "atomic" operations, that can be used as solid foundations to build the rest of the algorithm. These "atomic" operations are simple and can be robustly implemented:
\begin{itemize}
    \item (1) the algorithm will solely manipulate \emph{convex polytopes} and their intersections. A convex polytope $C$ is defined 
    as the intersection of $N_v$ half-spaces $\Pi^+_i$:
$$
    C = \bigcap_{i=1}^{N_v} \Pi^+_i = \bigcap_{i=1}^{N_v} \{ \bx \ | \ a_i x + b_i y + c_i z + d_i \le 0 \}
$$
    \item (2) to avoid inconsistencies, all the used predicates will be computed 
      in arbitrary precision \cite{DBLP:conf/compgeom/Shewchuk96}, with arithmetic filters to speed-up
      computations in the easy cases \cite{meyer:inria-00344297}.  
\end{itemize}

A consequence of $(1)$ is that the balls $B(\bx_i,\sqrt{\psi_i})$ will be approximated by convex polytopes:
$$
  B(\bx,R) \simeq \hat{B}(\bx,R) = \bigcap_{i=1}^{N_u} \{ \by \ | \ \| \bx - \by \|^2 \le \| \bx + 2 R \bu_i - \by \|^2 \}  
$$
where $(\bu_i)_{i=1}^{N_u}$ is a predetermined set of unit vectors $\bu$ that uniformly sample the unit sphere. 
While this strategy shares with the "ghost cells" approach the fact that it approximates the balls as convex polytopes, it   
differs from it in two important aspects: 
\begin{itemize}
\item first, computations are \emph{local} to a given
$V_i^\psi$ cell, with a small number of vertices (typically 50 to 200, as compared to the total number of points $N$), whereas the "ghost cells" 
strategy needs inserting all additional vertices into the \emph{global} Laguerre diagram, with $N+M$ vertices;
\item second, each $\hat{B}(\bx,R)$ can be instanced from a \emph{precomputed} unit sphere template $\hat{B}\left({\mathbf 0},1\right)$ without needing to compute any intersection, by scaling and translating the unit sphere template (no combinatorics needs to be computed). 
\end{itemize}

This means a large number $N_u$ (typically 100 to 200) of unit directions $\bu_i$ can be used without significantly 
increasing the computation time: the "partial optimal transport" algorithm computes a solution that is equivalent to using $N_u \times N$ "ghost cells", without the associated computational cost. Not only does it significantly reduce the complexity, but also it makes it possible to use a large value for $N_u$, that gives a precise approximation of the balls $B(\bx_i, \sqrt{\psi_i})$. Alternatively, one can also use a data structure for convex polytopes with spherical faces \cite{DBLP:journals/siamnum/LeclercMSS20}, but this comes with longer execution times, and robustness is more difficult to ensure.

The domain $\Omega$ can be either a convex polytope, or a simplicial mesh $\Omega = \bigcup_{i=1}^{N_T} T_i$. In both cases the intersections between 
$\Omega$, $Lag_i^\psi$ and $\hat{B}(\bx_i, \sqrt{\psi_i})$ involved in the definition of $V_i^\psi$ can be computed as intersections between convex polytopes.

\subsubsection{Convex polytopes: H and V representation}


The "atomic" operation computes the intersection between a convex polytope $C$ and a half-space $\Pi^+$:
$$
   C \leftarrow C \cap \Pi^+.
$$

\begin{figure}
    \centerline{
    \includegraphics[width=.49\columnwidth]{Figures/cell.png}\hfill
    \includegraphics[width=.49\columnwidth]{Figures/cell_dual.png}
    }
    \caption{Representation of convex polytopes.}
    \label{fig:cells}
\end{figure}


In the end, to estimate the gradient and Hessian of $K(.)$, one needs
to compute the volume of the cells and the areas of the cell borders.
To do so, we need to convert from the initial representation of the
cell, that is, as an intersection of half-spaces" (H-representation),
$(\Pi_i^+)_{i=1}^{Nv}$ into an explicit representation of all the
vertices and facets of the cell (V-representation). As shown in Figure
\ref{fig:cells}, if the planes are in generic position (we will
discuss degenerate cases further on), each vertex (black) is shared
by three planes $\Pi_i, \Pi_j, \Pi_k$. It is then natural to represent
the polytope in dual form (see Figure \ref{fig:cells}), by a
triangulation $(T_l = \{ i_l, j_l, k_l \})_{l=1}^{N_t}$ shown in grey
in the Figure (the same representation is exploited by the
Boywer-Watson algorithm that we use to compute the Laguerre diagram).
The coordinates $q_x, q_y, q_z$ of the vertex $\bq_l$ associated with
triangle $T_l$ are the solution of:
$$
\left[ \begin{array}{ccc}
   a_i & b_i & c_i \\
   a_j & b_j & c_j \\
   a_k & b_k & c_k
\end{array} \right]
\left[ \begin{array}{c}
  q_x \\ q_y \\ q_z
\end{array} \right]
= -
\left[ \begin{array}{c}
  d_i \\ d_j \\ d_k
\end{array} \right]
$$

Using Cramer's formula, the coordinates $q_x, q_y, q_z$ are given by:
\begin{equation}
\left[ \begin{array}{c}
  q_x \\ q_y \\ q_z
\end{array} \right] = 
{
-\frac{1}{\left| 
\begin{array}{ccc}
   a_i & b_i & c_i \\
   a_j & b_j & c_j \\
   a_k & b_k & c_k
\end{array} \right|}
\left[
\begin{array}{c}
\\
\left| 
\begin{array}{ccc}
   d_i & b_i & c_i \\
   d_j & b_j & c_j \\
   d_k & b_k & c_k
\end{array} \right| \\ \\
\left| 
\begin{array}{ccc}
   a_i & d_i & c_i \\
   a_j & d_j & c_j \\
   a_k & d_k & c_k
\end{array} \right| \\ \\
\left| 
\begin{array}{ccc}
   a_i & b_i & d_i \\
   a_j & b_j & d_j \\
   a_k & b_k & d_k
\end{array} \right| \\ \\
\end{array}
\right]}
\label{eqn:q}
\end{equation}
(the i-th component of the solution corresponds to the determinant of the system with its i-th column replaced with the r.h.s divided by the determinant of the system, see \cite{DBLP:journals/siamcomp/Aurenhammer87,DBLP:journals/tog/EdelsbrunnerM90} for similar computations).

\begin{figure}
    \centerline{
    \includegraphics[width=.32\columnwidth]{Figures/clip1.png}\hfill
    \includegraphics[width=.32\columnwidth]{Figures/clip3.png}\hfill
    \includegraphics[width=.32\columnwidth]{Figures/clip2.png}
    }
    \caption{Clipping a convex polytope by a half-space.}
    \label{fig:clipping}
\end{figure}

The algorithm to compute the intersection between a cell $C$ and a half-space $\Pi^+(a,b,c,d)$ works in three phases (see Figure \ref{fig:clipping}):
$$
\begin{array}{ll}
   \mbox{\textbf Input:}    &- \mbox{ A convex polytope } C  \mbox{ and a half-space } \Pi^+  \\
   \mbox{\textbf Result:}   &- \mbox{ replaces } C \mbox { with } C \cap \Pi^+   \\
\ (1): & \mbox{ classify the vertices } \\
\ (2): & \mbox{ discard the triangles } \\
\ (3): & \mbox{ triangulate the hole } \\
\end{array}
$$
\textbf{The first phase} determines the triangles that correspond to
vertices $\bq_l$ that are on the negative side of $\Pi^+$ (highlighted
in Figure \ref{fig:clipping}-Left), that is, $a q_x + b q_y + c
q_z + d < 0$. This can be rewritten as:
$$
{ %\tiny
\left|
\begin{array}{ccc}
  a_i & b_i & c_i \\
  a_j & b_j & c_j \\
  a_k & b_k & c_k 
\end{array}
\right|
\times
\left|
\begin{array}{cccc}
  a_i & b_i & c_i & d_i \\
  a_j & b_j & c_j & d_j \\
  a_k & b_k & c_k & d_k \\
  a   & b   & c   & d
\end{array}
\right|} \le 0
$$
by replacing the coordinates of $\bq$ by their expression (one obtains
an expression of the determinant above, developed w.r.t. to its last
row). Note that the vertices of the triangles can be consistently
oriented in such a way that the sign of the first term (the $3\times3$
determinant) is always positive. In the end, to classify a vertex, one
only needs to compute the sign of the $4 \times 4$ determinant $|
\Pi_i, \Pi_j, \Pi_k, \Pi |$ formed by the coefficients of the planes.
This is the only predicate (besides {\texttt orient3d} and {\texttt
in\_weighted\_sphere} used by the Bowyer-Watson algorithm that
constructs the Laguerre diagram) used by the algorithm. Note also
that during the computation of the intersections, only the equations
of the $\Pi_i$ planes are used, the $\bq$ points do not need to be
explicitly computed. To compute the exact sign of these determinants,
we use expansion arithmetics \cite{DBLP:conf/compgeom/Shewchuk96}, and
arithmetic filters \cite{meyer:inria-00344297} to quickly determine
the sign in the easy cases. The implementation in
\cite{DBLP:journals/cad/Levy16} generates the code that computes the
predicate from its formula. This strategy ensures that all
combinatoric decisions taken by the algorithm will remain consistent.
The points $\bq$ are computed once all clipping operations are
computed, right before evaluating the volumes and areas involved in
the gradient and Hessian of $K(.)$.
 
\textbf{The second phase} discards the triangles $(i,j,k)$ such that
$| \Pi_i, \Pi_j, \Pi_k, \Pi |$ $<$ $0$ (highlighted in Figure
\ref{fig:clipping}-Center). This leaves a hole in the triangulation
(the border of the hole is shown as dashed lines). As often done in
implementations of the Bowyer-Watson algorithm, the discarded
triangles are kept in a linked list, so that they can be reused in the
subsequent steps;

\textbf{The third phase} constructs a new triangle for each edge on the border of the hole created at the previous phase
(see Figure \ref{fig:clipping}-Right).

The algorithm has the same structure as the Bowyer-Watson algorithm, with the difference that the $\texttt in\_sphere$ predicate, that determines which triangle to discard, is replaced with ${\texttt sign}|\Pi_i, \Pi_j, \Pi_k, \Pi|$. Another difference is the way the Bowyer-Watson determines the list of triangles to be discarded. In our case, we test all the triangles. We could instead adapt the Bowyer-Watson strategy, that starts from a random vertex, and walks along the edges of the dual graph until it finds a triangle to be discarded, but since the number of triangles per polytope remains small in our case (typically $ < 100$), we did not observe a speedup using this technique. 

\subsubsection{Computing the cells $V_i^\psi$}

We shall now explain how to use the algorithm of the previous paragraph to compute the cells. Each cell 
$V_i^\psi = Lag_i^\psi \cap \Omega \cap \hat{B}(\bx_i, \sqrt{\psi_i})$ is the intersection between three
objects: the Laguerre cell $Lag_i^\psi$, the domain $\Omega$, represented by a simplicial mesh, that is,
a set of tetrahedra, and the (approximated) ball $\hat{B}(\bx_i, \sqrt{\psi_i})$. 
Note that the domain $\Omega$ is not necessarily convex, but since it is decomposed into simplices, we only
need to compute intersections between convex objects. The algorithm that computes a cell $V_i^\psi$ is detailed below:
 
$$
\begin{array}{ll}
   \mbox{\textbf Input:}    &- \mbox{ a simplicial mesh } \Omega \\
                        &- \mbox{ the set of points } (\bx_i)_{i=1}^N \mbox{ and the vector } (\psi_i)_{i=1}^N \\
                        &- \mbox{ the index } i \mbox{ of a cell} \\
   \mbox{\textbf Output:}   &- \mbox{ A set of convex polytopes } \mathcal{\mathbf C} = \{ C_k \} \mbox{ such that } V_i^\psi = \bigcup C_k\\[2mm]
%   \hline \\[4mm}
\end{array}
$$

$$
\begin{array}{ll}
%\hline \\
   \ (1):  & \mathcal{T} \leftarrow \{ t \subset \Omega \ | \ B(\bx_i,\sqrt{\psi_i}) \cap t \neq \emptyset \} \\
   \ (2):  & \mbox{If } \exists t \in \mathcal{T} \ | \ \partial t \cap \partial \Omega \neq \emptyset \\
   \ (3):  & \quad \mbox{For each } t \in \mathcal{T} \\
   \ (4):  & \quad \quad C \leftarrow t \\
   \ (5):  & \quad \quad C \leftarrow C \cap Lag_i^\psi \\
   \ (6):  & \quad \quad C \leftarrow \hat{B}(\bx_i, \sqrt{\psi_i}) \cap C \\
   \ (7):  & \quad \quad {\mathbf C} \leftarrow {\mathbf C} \cup \{ C \} \\
   \ (8): & \quad \mbox{End for}  \\
   \ (9): & \mbox{Else}  \\
   \ (10): & \quad C \leftarrow Lag_i^\psi \\
   \ (11): & \quad C \leftarrow \hat{B}(\bx_i, \sqrt{\psi_i}) \cap C \\
   \ (12): & \quad \mathcal{\mathbf C} \leftarrow \mathcal{\mathbf C} \cup \{ C \} \\
   \ (13): & \mbox{End if}   \\
\end{array}
$$

In line (1), we first determine the set of simplices $\mathcal{T}$ that have a non-empty intersection with the ball 
$B(\bx_i, \sqrt{\psi_i)}$. It is done using an axis-aligned bounding box tree (AABB tree), see e.g \cite{WEB:OPCODE}. 

Then, we distinguish two different cases: \textbf{first case:} if one of the intersected tetrahedra touches the boundary $\partial \Omega$ of the domain, then $V_i^\psi$ is not necessarily convex, but can be easily decomposed into convex objects: to do so, we compute the intersections $Lag_i^\psi \cap \hat{B}(\bx_i, \sqrt{\psi_i}) \cap t$ for each $t$  in $\mathcal{T}$. \textbf{Second case:} the ball $B(\bx_i,\sqrt{\psi_i})$ is included in $\Omega$ (no tetrahedron in $\mathcal{T}$ touches $\partial \Omega$), then $V_i^\psi = Lag_i^\psi \cap B(\bx_i, \sqrt{\psi_i})$ is convex and can be directly computed (lines (9) to (13)). 

In lines (6) and (11), we compute the intersection between a convex polytope and the approximated ball $\hat{B}(\bx_i, \sqrt{\psi_i})$. If $C$'s bounding sphere centered on $\bx_i$ has a radius smaller than $\sqrt{\psi_i}$ (that is, $\max_{\bq \in C} \{ \| \bx_i - \bq \|$ $<$ $\sqrt{\psi_i}$), then $C$ is entirely contained in the ball and one can skip the intersection computation. Else, the approximated ball $\hat{B}(\bx_i, \sqrt{\psi_i})$ is 
instanced (copied from the unit ball, translated and scaled), and then clipped by all the $\Pi_i$'s from the H-representation of $C$ (in this order, because $\hat{B}(\bx_i, \sqrt{\psi_i})$ has a larger number of facets than $C$ in general).
 
In line (4), the convex polytope $C$ is initialized from a tetrahedron of $\Omega$. Combinatorics are initialized from a fixed tetrahedron template, and the plane equation of each facet $\bp_1, \bp_2, \bp_3$ is
given by $ax + by + cz + d = 0$, where $(a,b,c) = (\bp_2-\bp_1) \times (\bp_3-\bp_1)$ and 
$d = -(a,b,c) \cdot \bp_1$.

In line (10), $C$ is initialized from a Laguerre cell. Combinatorics are copied from the Laguerre diagram, also stored in dual representation (tetrahedra). The triangles of $C$ correspond to the facets of the tetrahedra incident to vertex $i$ opposite to $i$. The plane equations $\Pi_{i,j}$ are given by: 
$
\bx \in \Pi_{ij} \quad \Leftrightarrow \quad
\| \bx - \bx_i \|^2 - \psi_i = \| \bx - \bx_j \|^2 - \psi_j
$
or: $ax + by + cz + d = 0$ where $(a,b,c) = 2(\bx_i - \bx_j)$ and $d = \bx^2_j -  \bx_i^2 - \psi_j + \psi_i$. \\

Up to now, we have supposed that all planes were in generic position (that is, each vertex of the polytope is shared by exactly three planes). To handle the specific configurations where more than three planes meet at a single vertex, one would think of using the symbolic perturbations, as in the "simulation of simplicity" technique \cite{DBLP:journals/tog/EdelsbrunnerM90}. In our case, the symbolic perturbation is trivial: one can use the predicate $|\Pi_i, \Pi_j, \Pi_k, \Pi| \le 0$ as is (provided that it is implemented with exact arithmetics, see \cite{DBLP:journals/cad/Levy16} and references herein). It means that each time a vertex $\bq$ is exactly located on the boundary of a clipping half-space $\Pi^+$ the associated triangle will be discarded. To evaluate the efficiency of this technique, degenerate configurations are tested in the next section.
 
It is not necessary to store the set of convex polytopes $\mathcal{\mathbf C}$: one can process them in a "streaming" manner, when assembling the gradient and Hessian of $K(.)$. Then, the algorithm stores a single $C$ at any time, and has only a small memory overhead besides the storage taken by the Laguerre diagram (stored as a 3D triangulation). Moreover, the algorithm is very simple to parallelize (then it stores one convex polytope per thread). When assembling the gradient and Hessian, it just needs a synchronization primitive (spinlock) to ensure that two threads do not modify the same entry at the same time.
 

\section{Fluid Simulation}

We shall now demonstrate how the algorithm can be used to implement a Lagrangian mesh with controlled volume that can change topology, with a very simple incompressible Navier-Stokes simulator for fluids with free boundaries.
In the context of this article, this "toy" simulator is meant to demonstrate that our Lagrangian mesh can reproduce some of the typical behaviors of a fluid while accurately tracking collisions and changes of topology. A more extensive quantitative evaluation of the simulator and the calibration of its parameters will be the subject of another article.

The simulator is obtained by adding viscosity and surface tension terms to the Gallouet-Merigot scheme \cite{Gallouet2017}, that simulates an incompressible Euler fluid. Inspired by Brenier and Benamou's point of view on incompressible fluids \cite{BrenierPFMR91,DBLP:journals/nm/BenamouB00}, the Gallouet-Merigot scheme softly projects the motion onto the incompressibility constraint by adding a "spring" force, that moves each point $\bx_i$ towards the centroid of its cell $V_i$. This additional "pressure" force makes the motion of the fluid tangent to the manifold of incompressible motions, and the computed fluid motion provably converges to the solution of an incompressible Euler fluid, that is, a geodesic on the manifold of incompressible motions (see \cite{WEB:Tao,Marsden} on geodesic flows and the Euler-Arnold equation). The reader is referred to \cite{Gallouet2017} for all the details.

The fluid is modeled as a set of $V_i$ cells, parameterized by the $\bx_i$ points. The volume of each $V_i$ cell is controlled and remains constant throughout the whole simulation. All the cells have the same mass $m$. Each cell is subjected to the following forces:

$$
  \begin{array}{lclcl}
     F_p(\bx_i) & = & \frac{1}{\epsilon_p^2} (\bx_i - \bg_i) 
     & & \mbox {where } \bg_i = \frac{1}{|V_i|} \int_{V_i} \bx dx  
     \\[4mm]
     F_g(\bx_i) & = & -m g {\mathbf z} \\[4mm]
     F_v(\bx_i) & = & \mu \hat{\Delta} \bv (\bx_i) & = & \mu \sum_{j \in N_i} \frac{1}{2\| \bx_j - \bx_i \|} \ |V_{ij}|\ (\bv_j - \bv_i) \\[4mm]
     F_t(\bx_i) & = & \gamma \hat{\Delta} \bx (\bx_i) & = & \gamma \sum_{j \in N_i} \frac{1}{2\| \bx_j - \bx_i \|} \ |V_{ij}|\ (\bx_j - \bx_i).
  \end{array}
$$
Following \cite{Gallouet2017}, the "pressure" $F_p$ may be thought of as a "spring" that connects each point $\bx_i$ to the centroid $\bg_i$ of the associated cell $V_i$ with a stiffness $\epsilon_p$. The additional forces are viscosity and surface tension: viscosity $F_v$ is proportional to the Laplacian of the velocity. Surface tension $F_t$ is derived in a volumetric manner, from the mutual attraction in the fluid. Inside the fluid, the net result is zero. On the boundary, the net result is a force that pulls towards the interior of the fluid. Both viscosity and surface tension use the $\mathbb{P}_1$ Laplacian: 
 
$$
 \hat{\Delta} f(\bx) = \sum_{j \in N_i} \frac{1}{2 \| \bx_j - \bx_i \| } \ | V_{ij} |\ (f(\bx_j) - f(\bx_i))
$$
where $|V_ij| = |\partial V_i \cap \partial V_j|$ denotes the area of the facet common to the cells $V_i$  and $V_j$, and where $N_i$ denotes the set of cells touching $V_i$. In the expressions of $F_v$ and $F_t$, it is computed component-wise.
 
The simulation results shown below use a semi-implicit time-stepping, with implicit integration for the viscosity, as follows:

$$
  \begin{array}{lcl}
    \bx^{(k+1)} & = & \bx^{(k)} + \delta t \bv^{(k)} \\[2mm]
    \bv^{(k+1)} & = & \bv^{(k)} + \frac{\delta t}{m} \left( 
       F_p + F_g + F_t + \mu \hat{\Delta} \bv^{(k+1)} 
    \right) \\
    \mbox{or: }  
    \left( \mbox{Id} - \mu \frac{\delta t}{m} \hat{\Delta}\right) \bv^{(k+1)} & = &
       \bv^{(k)} + \frac{\delta t}{m} \left( F_p + F_g + F_t\right) \\[2mm]
    \mbox{or: }  
    \left( 
      \hat{\Delta} - \frac{m}{\mu \delta t} \mbox{Id}
    \right) \bv^{(k+1)} & = & - \frac{1}{\mu} \left( 
       \frac{m}{\delta t} \bv^{(k)} + F_p + F_g + F_t 
    \right).
  \end{array}
$$
The linear system is solved for each component $x,y,z$ separately. Note that the matrix of this linear system corresponds to $\hat{\Delta}$ plus a diagonal term. To assemble it, the same code that constructs the Hessian of $K()$ can be reused. Sub-stepping is used whenever the CFL condition is violated, that is whenever a displacement $\delta t\ \bv$ is such that more than one cell is crossed, that is whenever there exists an $i$ such that $\delta_t \ |\bv_i| > \sqrt{\psi_i}$.
 
The method is tested in four different configurations:
\begin{itemize}
    \item a simulated crown splash, shown in Figure \ref{fig:crown}, computed with 500000 cells. The simulated mesh captures all the changes of topology, and preserves the fine details of the motion throughout the simulation;
    % \item  (Figure \ref{fig:honey}), demonstrating the interactions between 
    % the fluid and the boundary of the domain;
    \item a fluid in a zigzag domain (Figure \ref{fig:mercury}) with a higher surface tension;
    \item a simulation of the experimental results in \cite{hurdles} shown in Figure \ref{fig:hurdlesreal}, where a droplet slides along a slope and bounces over obstacles. An example of a simulation mesh
    is shown in Figure \ref{fig:hurdles2}, as well as a cross-section, revealing the internal structure of the mesh, and the effect of surface tension that clusters the points $\bx_i$ near the surface (due to the volume constraint, the associated cells are more elongated than the others). Some frames of the full animation are shown in Figure \ref{fig:hurdles}. The fluid-boundary interactions and the typical oscillations of the droplet are reproduced in the simulation.
\end{itemize}
 
The table below indicates the parameters for the four simulations: 
$$
   \begin{array}{l|l|l|l|l|l|c}
     & g & m & \mu & \gamma & N & \mbox{time (per timestep)} \\
   \hline
   \mbox{Splash}   & 10  & 3 & 0.001 & 1.5 & 500000 & 98s \\
   \mbox{Zigzag 1} & 1   & 3 & 0.001 & 0.5 & \ 25000 & 12s \\
   \mbox{Zigzag 2} & 1   & 3 & 0.001 & 2   & \ \ 5000 & 2.8s \\
   \mbox{Hurdles}  & 0.3 & 3 & 0.001 & 3   & \ \ 5000 & 3.5s 
   \end{array}
$$
The table reports the parameters for gravity ($g$), viscosity ($\mu$), surface tension ($\gamma$), the number of cells ($N$) and the computation time for one timestep. The domain $\Omega$ is normalized in the $[0,1]^3$ box. For fine-scale tension-dominated simulations (zigzag and hurdles), the simulation space is globally scaled-down by applying a scaling factor to $g$. In the four simulations, the threshold for Newton iterations $\epsilon_K$ is set to $1\%$. During the whole simulation, each cell volume $|V_i|$ remains within a $1\%$ tolerance from the specified volume. Note that unlike some methods that use an Eulerian mesh, volume errors do not accumulate, the $1\%$ bound is met by all cells at all timesteps. The parameter $\epsilon_p$ for the "pressure" is set to $0.004$, and the timestep is set to $0.004$. The spheres are approximated using $N_u = 162$ half-spaces.

We shall now report timing breakdown, measured with the crown splash example (Intel Core i9-9880H, 2.3 GHz, multithreaded implementation):
$$
\begin{array}{l|l}
  \mbox{phase} & \mbox{time} \\
  \hline
  \mbox{Laguerre (Bowyer-Watson)}   & 9.7\% \\
  \mbox{Linear solve} & 2.2\% \\
  \mbox{Evaluate gradient} & 32.8\% \\
  \mbox{Evaluate and assemble Hessian}  & 41.5\% \\
  \mbox{Euler update and implicit viscosity lin. solve} & 13.8\%
\end{array}
$$
The phases that take the largest amount of time are evaluating the gradient and Hessian. It is because these phases involve all the geometric computations (query the Axis-Aligned Bounding Box tree, compute intersections between convex polytopes and compute volumes and areas).

Concerning the geometric predicate used by the convex polytope clipping algorithm (the sign of a 4x4 determinant), during 100 timesteps of the crown splash example, it was called 15752996161 times (15.7 billions), and the arithmetic filter triggered the arbitrary precision mode 70034 times (that is, 0.0004 \%). Clearly, degenerate configurations represent a tiny proportion, but given the huge number of computed intersections, they appear in each simulation. Our robust intersection algorithm properly handles them, with negligible overhead, thanks to the arithmetic filters (see also \cite{DBLP:journals/cad/Levy16} and references herein).

% \cph{Conclusion and introduction for the practicals?}
 

\begin{figure}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/crown0}    
    \includegraphics[width=.495\columnwidth]{Figures/crown1}
    }
    \vspace{2mm}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/crown2}
    \includegraphics[width=.495\columnwidth]{Figures/crown3}
    }
    \vspace{2mm}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/crown4}
    \includegraphics[width=.495\columnwidth]{Figures/crown5}
    }
    \caption{Testing surface tension and topology changes: crown splash}
    \label{fig:crown}
\end{figure}

\begin{figure}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/mercury_1}
    \hspace{2mm}
    \includegraphics[width=.495\columnwidth]{Figures/mercury_2}
    }
    \vspace{2mm}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/mercury_3}
    \hspace{2mm}
    \includegraphics[width=.495\columnwidth]{Figures/mercury_4}
    }
    \caption{Free-surface flow and interactions with boundaries: higher surface tension.}
    \label{fig:mercury}
\end{figure} 

\begin{figure}
    \centerline{
    \includegraphics[width=\columnwidth]{Figures/hurdles}    
    }
    \caption{Droplet hurdles race, courtesy Hélène de Maleprade, Rachid Bendimerad, Christophe Clanet and David Quéré.}
    \label{fig:hurdlesreal}
\end{figure}

\begin{figure}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/voro_cut_3}    
    \includegraphics[width=.495\columnwidth]{Figures/voro_cut_1}
    }
    \caption{Simulated droplet hurdles race: Mesh details and cross-section.}
    \label{fig:hurdles2}
\end{figure}
 

\begin{figure}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/hurdles1}    
    \includegraphics[width=.495\columnwidth]{Figures/hurdles2}
    }
    \vspace{2mm}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/hurdles3}    
    \includegraphics[width=.495\columnwidth]{Figures/hurdles4}
    }
    \vspace{2mm}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/hurdles5}    
    \includegraphics[width=.495\columnwidth]{Figures/hurdles6}
    }
    \vspace{2mm}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/hurdles7}    
    \includegraphics[width=.495\columnwidth]{Figures/hurdles8}
    }
    \vspace{2mm}
    \centerline{
    \includegraphics[width=.495\columnwidth]{Figures/hurdles9}    
    \includegraphics[width=.495\columnwidth]{Figures/hurdles10}
    }
    \caption{Simulated droplet hurdles race}
    \label{fig:hurdles}
\end{figure}

\bibliographystyle{eg-alpha-doi} 
\bibliography{bibliography}

\end{document}
