%~Mouliné par MaN_auto v.0.27.3 2023-09-27 10:10:46
\documentclass[CRMECA, Unicode, XML, screen, Thematic, HideArticleType]{cedram}
%
%\TopicEN{Section information missing}
%\TopicFR{Information section manquante}
\usepackage{xcolor}
\usepackage[all]{xy}
\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bbE}{\mathbb{E}}
\newcommand{\bbP}{\mathbb{P}}
\newcommand{\setR}{\mathbb{R}}
\newcommand{\setN}{\mathbb{N}}
\renewcommand*\O {\mathcal{O}}
\newcommand{\clM}{\mathcal{M}}
\newcommand{\clL}{\mathcal{L}}
\newcommand{\diag}{\mathrm{diag}}
\newcommand{\flux}{\mathrm{flux}}
\newcommand{\corrl}{\mathrm{co}}
\newcommand{\OO}{\Omega}
\newcommand{\p}{\partial}
\newcommand{\ps}{\partial^\star}
\newcommand*\intV {\mathring{V}}
%\newcommand*\<[1],#2>{{\left<#1\mathbin,#2\right>}}
%\newcommand*\s<[1],#2>{{\left\langle #1\,|\,#2\right\rangle}}
%\newcommand*\<[1],#2>{{\left<#1\,,\,#2\right>}}
\newcommand*\set[1] {\left\{#1\right\} }
\newcommand*\abs[1]{\left|#1\right|}
\newcommand*\virg { \,, \; }
\newcommand*\dsp {\displaystyle}
\newcommand*\vseq {\vspace{4pt} }
\newcommand*\oB {\mathring{B}}
\newcommand*\op {\mathring{p}}
\newcommand*\intq {\mathring{q}}
\newcommand*\intp {\mathring{p}}
\makeatletter
\def\cdr@preabstracthook{%
\vspace*{-5pt}
\vskip2pt plus 3pt minus 1pt
}
\makeatother
\newcounter{casecount}
\renewcommand*\thecasecount{\roman{casecount}}
\newenvironment{case}{\refstepcounter{casecount}\begin{proof}[{Case~\thecasecount}]}{\let\qed\relax\end{proof}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}
\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}
%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}
\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\RubricEF{The scientific legacy of Roland Glowinski}{L'héritage scientifique de Roland Glowinski}
\title{Defective Laplacians and paradoxical phenomena in crowd motion modeling}
\alttitle{Laplaciens déficients et phénomènes paradoxaux dans la modélisation des mouvements de foule}
\author{\firstname{Bertrand} \lastname{Maury}}
\address{LMO, Université Paris-Saclay, F-91405 Orsay cedex\\
\& DMA, Ecole Normale Supérieure, PSL University, Paris}
\email{bertrand.maury@universite-paris-saclay.fr}
\keywords{Discrete Laplace operator, maximum principle, crowd motion, faster-is-slower effect}
\altkeywords{Laplacien discret, principe du maximum, mouvements de foules, effet faster-is-slower}
\begin{abstract}
In both continuous and discrete settings, Laplace operators appear quite commonly in the modeling of natural phenomena, in several context: diffusion, heat propagation, porous media, fluid flows through pipes, electricity\dots. In these contexts, the discrete Laplace operator enjoys all the properties of its continuous counterpart, in particular: self-adjointness, variational formulation, stochastic interpretation, mean value property, maximum principle, \dots In a first part, we give a detailed description of the correspondence between these mathematical properties and modeling considerations, in contexts where the continuous and the discrete settings perfectly match. In a second part, we describe a pathological situation, in the context of granular crowd motion models. Accounting for the non-overlapping constraint between hard discs leads to a particular operator acting on a field of Lagrange multipliers, defined on the dual graph of the contact network. This operator is defective in a certain sense: although it is the microscopic counterpart of the macroscopic Laplace operator, this discrete operator indeed lacks some properties, in particular the maximum principle. We investigate here how this very \emph{defectivity} may explain some \emph{paradoxical} phenomena that are observed in crowd motions and granular materials, phenomena that are not reproduced by macroscopic models.
\end{abstract}
\begin{altabstract}
Aux niveaux continu et discret, l'opérateur de Laplace intervient de façon très courante dans la modélisation de phénomènes naturels, dans de nombreux contextes: diffusion, propagation de la chaleur, milieux poreux, écoulement de fluides dans des conduits, électricité\dots. Dans ces contextes, le laplacien discret possède toutes les propriétés de son pendant continu: caractère auto-adjoint, structure variationnelle, interprétation stochastique, propriété de la valeur moyenne, principe du maximum, \dots. Dans une première partie, nous proposons une description détaillée des liens entre ces propriétés et les aspects de modélisations, dans des contextes où les notion continues et discrètes se correspondent parfaitement. Dans une seconde partie, nous décrivons une situation pathologique, dans le contexte de la modélisation de foules d'un point de vue granulaire. La prise en compte de la contrainte de non recouvrement entre grains rigides conduit à un opérateur particulier qui agit sur les champs de multiplicateurs de Lagrange, définis sur le graphe dual du réseau de contacts. Cet opérateur est déficient dans un certain sens~: bien qu'il apparaisse comme le pendant discret du laplacien continu, il ne vérifie pas certaines des propriétés usuelles du laplacien, en particulier le principe du maximum. Nous explorons comment cette déficience permet d'expliquer certains effets paradoxaux observés en mouvements de foules, que les modèles continus ne reproduisent pas.
\end{altabstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\dateposted{2023-12-21}
\begin{document}
\maketitle
\section{Introduction, continuous and discrete Laplacians}\label{sec:intro}
The Laplace operator occupies a central position in the modeling of physical phenomena. In many situations, it combines a phenomenological phenomenon together with a conservation principle. In the context of diffusion of some substance represented by its density $p= p~(x,t)$, Fick's Law states that the flux vector $J$ is proportional to the gradient of $p$:
\[
J = - D \nabla p,
\]
where $D>0$ is the diffusion coefficient. The flux vector $J$ is such that, for any oriented surface $\Sigma$, with normal vector $n$, $\int_\Sigma J\cdot n$ is the flux of substance through $\Sigma$, counted positively in the\linebreak $n$-direction. Now considering a domain $\omega$, assuming that there is a source term $f$, mass balance writes
\[
\frac {d}{dt} \int_\omega p = \int_\omega {\partial _t}p = \int_{\omega } f - \int_{\partial \omega } J \cdot n = \int_{\omega } f - \int_{\omega } \nabla \cdot J
\]
which leads to
\[
\int_\omega \left ({\partial _t}p + \nabla \cdot J - f \right) = 0.
\]
Since $\omega $ is arbitrary, and $J = - D \nabla p$, it yields the heat equation
\[
{\partial _t}p - D \Delta p = f,
\]
the stationary version of which is the Poisson equation
\[
- D \Delta p = f,
\]
set in a domain $\Omega$ occupied by the substance. Let us mention another context in which $-\Delta$ appears as a composition of $\nabla$ (gradient) and $-\nabla\cdot $ (minus the divergence operator). We consider the flow of an incompressible viscous fluid in a porous medium. Darcy's law writes
\[
u = - k \nabla p
\]
where $p$ is the pressure, $k$ the permeability of the medium, and $u$ is the Darcy velocity (which should be thought of as a volume flux, expressed in cubic meter per square meter per second, which is indeed homogeneous to a velocity). Volume conservation writes
\[
-\nabla \cdot u = - \nabla \cdot k \nabla p = 0,
\]
which is again the Laplace equation in a possibly non homogeneous medium.
\begin{figure}[!htbp]
\includegraphics[width=0.75\linewidth]{DiffCell.pdf}
\caption{Diffusion process between cells.}\label{fig:DiffCell}
\end{figure}
To obtain a discrete Laplacian, one may e.g. consider a set of cells filled with some medium (see Figure~\ref{fig:DiffCell}). A substance freely diffuses in each cell, so that the concentration within a cell can be considered as uniform. Cells are separated by a membrane, and the flux of substance between two cells is assumed to be proportional to the difference between concentrations. We define the abstract non-oriented graph $(V,E)$ as follows : the set of vertices $V$ is the set of cells (symbolized by black dots in the figure), and the set of edges $E$ contains $(x,y)$ if and only if cells $x$ and $y$ share a piece of membrane. Note that $E$ is symmetric by construction: the graph is said to be \emph{non oriented}. We shall write $x\sim y$ whenever $(x,y)\in E$. If we denote by $p_x$ the concentration (or partial pressure) on the cell $x$, and $c_{xy}$ the permeability of the membrane ($c_{xy}$ plays the role of a \emph{conductance} in electrical / fluid networks), the flux from $y$ to $x$ writes
\[
J_{y \to x} = c_{xy} (p_y-p_x).
\]
Note that the flux is symbolically associated to the edge connecting $x$ and $y$, but it corresponds here to a transfer that is spread over the interface which separates $x$ and $y$ cells (grey arrows in the figure).
Denoting by $M_x$ the volume of cell $x$, the product $M_x \, p_x$ is the total quantity of substance at $x$, and mass balance at each cell $x$ writes
\[
M_x \frac {dp_x}{dt} = - \sum_{y\sim x} c_{xy} \left(p_x - p_y\right).
\]
It can be written globally
\[
M \frac {dp}{dt} + Lp = 0,
\]
where $L$ is a \emph{discrete Laplacian}, i.e. the discrete counterpart of $-\Delta$, defined as
\[
L\; : \; p \in \bbR^V \longmapsto Lp = \left (\sum_{y\,\sim\,x} c_{xy} \left(p_x - p_y\right)\right) _{x\,\in\,V},
\]
and $M = \diag(M_x)$.
In another context, natively static, one may also consider a network of interconnected pipes, in which a viscous fluid is flowing. If we denote by $V$ the set of bifurcation points, by $E \subset V\times V$ the (symmetric) set of pipes between vertices, the flow through a pipe $(x,y)$ follows the so-called Poiseuille's law
\[
u_{xy} = c_{xy} \left(p_x - p_y\right),
\]
where $c_{xy}$ is the conductance of the edge $(x,y)$, which depends on geometrical properties of the actual pipe that it represents. Now consider that the flow in conservative at some vertices, the set of which is denoted by $\intV$, and possibly non conservative on the set of remaining points $\Gamma = V \setminus \intV$ (see Figure~\ref{fig:network}). Conservation on $\intV$ writes
\[
\sum _{y\,\sim\,x } u_{xy} = \sum _{y\,\sim\,x } c_{xy} \left(p_x - p_y\right) = 0 \quad \forall x \in \intV.
\]
As a consequence, if the value of the pressure is prescribed at boundary points (field $P\in \bbR^\Gamma$), the pressure field is the solution to a Dirichlet problem
\begin{equation*}
\left |
\begin{array}{rcl}
Lp& = &0 \quad \text{ in } \intV \\
p &=& P \quad \text{ on } \Gamma.
\end{array}
\right.
\end{equation*}
This setting has been intensively used to model the airflow in the respiratory tree (see e.g.~\cite{RSEMaury13} or~\cite{SGT20}). It can be transposed to resistive electrical networks, where the potential plays the role of the pressure, and electric intensity the role of fluid fluxes.
\begin{rema}\label{rem:droniou}
Let us make it clear that the discrete equations written above do not aim here at approximating the solution to continuous Poisson or Laplace equations, they model intrinsically discrete phenomena, and so will the discrete Laplace operator introduced in Section~\ref{sec:crowd}. Let us nevertheless mention the deep link between the diffusion between cells which we considered (Figure~\ref{fig:DiffCell}) and the elaboration of numerical methods of the Finite Volume type, for which cells are virtual, since they do not correspond in general to a specific physical zone of the underlying domain. We refer to~\cite{droniou14} for a very detailed and documented account of discretization strategies of the Finite Volume type. Requiring that the discrete operator verify the maximum is essential to establish $\ell^\infty$ stability of the numerical scheme. Note that, in the context of approximating PDE's, convergence also relies on consistency with respect to the continuous operator, which calls for strong requirements on the mesh that is used. In the present setting, as detailed in the next section, the maximum principle is automatically verified since the models are compliant with the principles of thermodynamics, whereas consistency with respect to the continuous Laplacian is not an issue.
\end{rema}
\begin{figure}[tbp]
\includegraphics[width=0.8\linewidth]{network.pdf}
\caption{Network $(V,E,c,\Gamma)$}\label{fig:network}
\end{figure}
\section{Mathematical properties and their physical interpretations}\label{sec:propslap}
We recall in this section the main properties of the standard Laplace operator on euclidean domains, and formulate their discrete counterparts for resistive networks. Some of the results at the continuous level require the domain to be \emph{smooth} in some sense. We shall simply call smooth a domain with Lipschitz boundary.
\subsection{Self-adjointness, positivity, variational setting}\label{sec:savs}
Consider the case of Dirichlet boundary conditions, i.e.\ $p$ is set to $0$ on $\partial \Omega$. From a functional analysis standpoint, $ - \Delta$ is an unbounded operator in $L^2(\Omega)$, which is self-adjoint, thanks to Green's formula:
\[
\left\langle - \Delta p\,\middle|\,q \right\rangle = - \int_\Omega q \Delta p =
\int_\Omega \nabla p \cdot \nabla q = \left\langle - \Delta q\,\middle|\, p\right\rangle.
\]
Its self adjoint character comes from the fact that ``$-\nabla \cdot $'' (opposite of the divergence) and $\nabla $ (gradient) are mutually adjoint for the $L^2$-- duality. Note that, although it is intensively used in the theoretical study of linear elliptic equations, this fact is in some way \emph{fortuitous}, and makes a connections between considerations of very different natures. On the one hand, the divergence operator expresses a universal conservation principle, sometimes called Lavoisier' Principle. It would not make any sense to imagine an alternative way, e.g. nonlinear, to express this conservation principle. On the other hand Fick's law is phenomenological, it was inferred from experimental measurements, and in some contexts it takes a nonlinear form (see e.g.~\cite{nonFick15}).
This makes it possible to formulate the Dirichlet problem as a minimization problem, the Poisson equation corresponding to optimality conditions. It can be expressed in the form of a very classical existence and uniqueness property, which we shall write extensively anyway to emphasize the similarity with the discrete context (see Proposition~\ref{prop:FVdirdisc} thereafter).
\begin{prop}\label{prop:FVdir}
Let $\Omega$ be a bounded smooth domain, and let $P\in H^{1/2}(\Gamma)$ be given. Then the Dirichlet problem
\begin{equation}\label{eq:dir}
\left |
\begin{array}{rcl}
-\Delta p & = &0 \quad \text{ in } \Omega\\
p &=& P \quad \text{on } \Gamma.
\end{array}
\right.
\end{equation}
admits a unique weak solution in $H^1(\OO)$, in the sense that there exists a unique $p \in H^1_P(\OO)$ such that
\[
\int _\OO \nabla p\cdot \nabla q = 0 \quad \forall q\in H^1_0(\OO),
\]
with
\[
H^1_P(\OO) = \set{q\in H^1(\OO)\virg q_{|\Gamma} = P },
\]
where $q_{|\Gamma} $ denotes the trace of $q$ on $\Gamma$. Furthermore, $p$ is the unique minimizer of
\[
v \longmapsto J(v) = \frac 1 2 \int_\OO \abs {\nabla q} ^2
\]
in $H^1_P$.
\end{prop}
\begin{proof}
See e.g.~\cite{LionsNH72}.
\end{proof}
In the discrete setting, a similar framework can be set. Let us consider a resistive network $(V,E,c,\Gamma)$, where $c\in]0,+\infty[^E$ is the collection of conductances, and $\Gamma$ is the subset of vertices which may exchange fluid (or intensity in the electrical context) with the outside world. We shall consider pressure fields $p\in \setR^V$ and flux fields $u \in \setR^E$, with the convention that fluxes are skew symmetric, i.e. $u_{xy} = - u_{yx}$. As detailed in~\cite{RSEMaury13}, we can define a discrete divergence operator (it actually stands for the \emph{opposite} of the divergence):
\begin{align*}
\partial \;: \; u \in \bbR^E& \longmapsto \partial u \in \setR^V \\
(\partial\, u)_x& = \sum_{y\,\sim\,x} u_{yx}.
\end{align*}
The conservation at interior vertices simply writes $(\partial\, u)_x = 0 $ for any $x\in \intV$.
In the same spirit, we define a discrete gradient (or ``drop'') operator as
\begin{align*}
\ps \;: \; p \in \setR^V &\longmapsto \ps p \in \setR^E \\
\left(\ps p\right)_{xy} &= p_y - p_x.
\end{align*}
Poiseuille / Ohm's law then writes $u = - c \ps p$.
The two operators $ \p \;: \; \setR^E\longrightarrow \setR^V $ and $ \ps \;: \; \setR^V\longrightarrow \setR^E $ defined above can be straightforwardly checked to be mutually adjoint, as stated by the next proposition.
\begin{prop}\label{prop:divgrad}
We have that
\[
\left\langle \p u\,\middle|\, p\right\rangle _V = \left\langle u\,\middle|\,\ps p\right\rangle_E\quad \forall u \in \setR^E,\, p \in \setR^V,
\]
where $\langle \cdot\,|\,\cdot\rangle_V$ and $\langle \cdot\,|\,\cdot \rangle_E$ are the canonic scalar products on $\setR^V$ and $\setR^E$, respectively.
\end{prop}
\begin{proof}
For $u \in \setR^E$, $p\in \setR^V$, it holds that
\[
\left\langle\p u\,\middle|\, p\right\rangle_V =
\sum _{x\,\in\,V} p_{x} (\p u)_x =
\sum _{x\,\in\,V} p_{x} \sum_{y\,\sim\,x} u_{yx} =
\sum_{(x,y)\,\in\,E} u_{xy} \underbrace{\left(p_y - p_x\right)}_{\ps p~(x,y)} =
\left\langle u\,\middle|\,\ps p\right\rangle_E,
\]
which ends the proof.
\end{proof}
\begin{figure}[tbp]
\[
\xymatrix{
u \in\bbR^E (1-\text{ chains}) \ar[r]^\partial & \partial u \in\bbR^V (0-\text{ chains})\\
\partial^\star p \in\bbR^E (1-\text{ cochains}) & \ar[l]^{\quad\partial^\star} p \in\bbR^V (0\text{ cochains})}
\]
%\includegraphics[width=0.72\linewidth]{VV1.pdf}
\caption{Abstract setting: chains and cochains.}\label{fig:VVchain}
\end{figure}
\begin{rema}
In the language of algebraic topology, the operators $\p$ and $\ps$ can be considered as boundary and coboundary operators, respectively, whereas fluxes are $1$-chains, mass conservation defects are $0$-chains, pressure fields are $0$-cochains, and pressure drops are $1$-cochains (see Figure~\ref{fig:VVchain}). The duality between $1$-chains and $1$-cochains (on the left-hand side of the figure) corresponds to the power dissipated within edges, whereas the duality between $0$-chains and $0$-cochains corresponds to the power of external pressures (i.e. pressure at non-conservative vertices where some fluid is exchanged with the outside world). In the spirit of Remark~\ref{rem:droniou}, it is fruitful to establish a link between this abstract setting and discretization strategies for Partial Differential Equations. We refer to~\cite{rapetti09} for a detailed description of how principles borrowed from exterior calculus can be integrated in the elaboration of Finite Element discretization strategies which are respectful of the underlying equations, from this standpoint.
\end{rema}
The problem takes the form of a discrete Darcy problem, which combines a phenomenological law (Poiseuille, Fick, or Ohm), together with a conservation principle:
\begin{equation}\label{eq:darcydisc}
\left |
\begin{array}{lcl}
u + c\ps p & = &0 \quad \text{ over } E \\
\p u &=& 0 \quad \text{ in } V.
\end{array}
\right.
\end{equation}
Like in the continuous setting, the discrete Laplace operator appears through the elimination of the velocity in the Darcy system. It can be written $L = \p c \ps $ where, for any $p\in \setR^V$, $c \ps p$ simply stands for the term-by-term product of $c$ and $\ps p$, i.e. $c \ps p\in \setR ^E$, with $(c \ps p)_e = c_e \, (\ps p) _e$. Its self-adjointness is a straight consequence of Proposition~\ref{prop:divgrad}:
\[
\left\langle Lp\,\middle|\,q\right\rangle = \left\langle \p c \ps p\,\middle|\,q\right\rangle = \left\langle c \ps p\,\middle|\,\ps q\right\rangle = \left\langle \ps p\,\middle|\,c\ps q\right\rangle = \left\langle p\,\middle|\,\p c\ps q\right\rangle =\left\langle p\,\middle|\,Lq\right\rangle.
\]
Like in the continuous setting again, this property can be seen as fortuitous in terms of modeling, since it relies on a strong relation between a universal principle (mass conservation) and a phenomenological law (Poiseuille's in a fluid context, or Ohm's law for a electrical network).
Proposition~\ref{prop:FVdir} has a discrete counterpart: the Dirichlet problem is well-posed under general assumptions, and the problem can be formulated in a variational way.
\begin{prop}\label{prop:FVdirdisc}
Let $(V,E,c,\Gamma)$ be a connected resistive network, with boundary $\Gamma\neq \emptyset$, and $P\in \setR^\Gamma$ a collection of boundary value. Then the Dirichlet problem
\begin{equation}\label{eq:dirdisc}
\left |
\begin{array}{rcll}
\left(L p\right)_x = \dsp \sum_{y\,\sim\,x}c_{xy} \left(p_x-p_y\right) & = &0 & \forall x \in \intV\vseq\\
p_x &=& P _x& \forall x \in \Gamma.
\end{array}
\right.
\end{equation}
admits a unique solution in $\setR^V$. It is the unique field $p\in H_P$ such that
\[
\sum_{(x,y)\,\in\,E }c_{xy} \left(p_x-p_y\right) \left(q_x-q_y\right) = 0 \quad \forall q \in H_0,
\]
with
\[
H_P = \set{q\in \setR^V\virg q_{|\Gamma} = P }.
\]
Furthermore, $p$ is the unique minimizer of
\[
q \longmapsto J(q) = \frac 1 2\sum_{(x,y)\,\in\,E }c_{xy} \left(q_x-q_y\right)^2
\]
in $H^1_P$.
\end{prop}
\begin{proof}
Any field $q$ of $\setR^V$ can be written $(\intq,q_\Gamma)$ so that minimizing $J$ over $H_P$ amounts to minimizing
\[
\intq \longmapsto \Phi\left(\intq\right) = J\left(\intq,P_\Gamma\right)
\]
over $\setR^{\intV}$. Let us show that $\Phi$ is coercive. Let $q = (\intq,P_\Gamma)$ be given in $H_P$. Since the graph is connected and $\Gamma\neq\emptyset$, any $x\in \intV$ can be connected to $\Gamma$ by a path
\[
x_0\in \Gamma \sim x_1 \sim\dots \sim x_n = x.
\]
It holds that
\begin{align*}
\abs { q_x} &= \abs {q_{x_0} + \sum_{j= 1}^{n} \left(q_{x_{j-1}} - q_{x_{j}}\right) }
\leq \abs{q_{x_0} } + \sum_{j = 1}^{n} \abs { q_{x_{j-1}} - q_{x_{j}} }
\\
&\leq \abs{P_{x_0} } + \underbrace{\left (\sum_{j = 1}^{p} c_{x_{j-1}x_j} \abs { q_{x_{j-1}} - q_{x_{j}} }^2 \right) ^{1/2}}_{\leq\,J(\intp) ^{1/2}}
\left (\sum_{j = 1}^{p} \frac 1 {c_{x_{j-1}x_j}} \right) ^{1/2}.
\end{align*}
Boundedness of $\Phi(\intq)$ therefore implies boundedness of $q$ (for any norm in this finite dimensional context). The functional $\Phi$ is therefore coercive: $\Phi(\intq)\rightarrow +\infty$ when $\abs {q}\rightarrow +\infty$. Since it is continuous, it admits a minimizer $\intp$, at which $\nabla \Phi $ vanishes, so that $p = (\intp,P_\Gamma)$ solves~\eqref{eq:dirdisc}, which ensures the existence of a solution. Let us now prove that $\Phi $ is strictly convex. Considering $\intp $ and $\intp'$ two different fields, we consider $p = (\intp,P_\Gamma)$ and $p' = (\intp',P_\Gamma)$. Since $p\neq p'$, while they identify at least at one vertex, and since the graph is connected, there exist at least an edge $(x,y)$ on which $p_x - p_y \neq p'_x - p'_y $. The associated quadratic term of $\Phi$ therefore ensures strict convexity. The minimizer of $\Phi$ is therefore unique. By convexity of $\Phi$, any field $\intp$ that vanishes the gradient (i.e. such that $(\intp,P_\Gamma)$ is a solution to~\eqref{eq:dirdisc}) is a minimizer, which ensures that the solution to~\eqref{eq:dirdisc} is unique.
\end{proof}
Let us insist on the fact that, in the discrete setting, the notion of boundary is somewhat arbitrary, whereas for a euclidean domain its designs the topological boundary $\overline \Omega \setminus \Omega$. The latter expression makes no sense for a network, and $\Gamma$ is simply the set of non conservative vertices. We may even consider the caricatural case $\Gamma = V$, $\intV = \emptyset$, which leads to a trivial Dirichlet problem $p\equiv P$.
In both the continuous and the discrete settings, the positivity of the Laplace operator that results from the Darcy system expresses dissipation of energy. For the continuous Darcy problem, $c \abs{\nabla p}^2$ is the dissipated power per unit volume (or area in 2d), due to viscous stress in the fluid. At the discrete level, $(c \abs { \ps p} ^2)_e$ is the dissipated power within edge $e$, through dissipative stress in the fluid context, or through Joule's effect in the electrical context ($p$ then represents an electrical potential).
\subsection{Mean value property and maximum principle}
The Laplacian of a function at some point measures the discrepancy between the value at this point and the mean value around it, as stated by the following proposition.
\begin{prop}\label{prop:mvp}
Let $\OO\in \setR^d$ be a domain (nonempty open set). A scalar field $p$ is an harmonic function in $\OO$, i.e. $-\Delta p \equiv 0$ in $\Omega$, if and only if for any $x\in \OO$, any $r>0$ such that $\overline B(x,r)\subset \OO $,
\[
p(x) = \frac 1{\abs{B(x,r)}}\int_{B(x,r)} p(y)\, dy =
\frac 1 {\abs{S(x,r)}}\int_{S(x,r)} p(y)\, dy,
\]
where $\abs{B(x,r)}$ is the $d$-volume of the ball $B(x,r)$, and $\abs {S(x,r)}$ the $d-1$ dimensional measure of the sphere $S(x,r) = \partial B(x,r)\in \setR^d$.
\end{prop}
\begin{proof}
See~\cite[Chapter~I]{Doob84}.
\end{proof}
Beyond this characterization of harmonic functions, the Laplacian operator estimates the local gap between the value of a function at some point $x$ and the mean value in the neighborhood of $x$. This general remark can be formulated as follow :
\begin{prop}\label{prop:lapmean}
Let $p$ be a $C^4$ scalar field in an open set containing $x\in \setR^d$. Then
\[
- \Delta p~(x) =
\frac 1 {\epsilon ^2}\left (p(x) -
\frac 1 {\abs{S(x,\epsilon)}} \int_{S(x,\epsilon)} p(y)\, dy \right) + \O\left(\epsilon^2\right).
\]
\end{prop}
\begin{proof}
For $\epsilon >0$ sufficiently small, for any $\sigma\in S(0,1)$, it holds that
\[
f(x + \epsilon \sigma)- f(x) =
\sum_{\abs \alpha = 1}^3 \frac {\partial ^\alpha f }{\partial x^\alpha } (x) \, \epsilon^{\abs \alpha}\, \sigma ^\alpha + \O\left(\epsilon^4\right),
\]
with $\abs{\alpha} = \alpha_1+\dots+ \alpha _d$ for $\alpha = (\alpha_1,\dots,\alpha_d)\in \setN^d$, and $\sigma^\alpha = \sigma_1^{\alpha_1} \times \dots \sigma_d^{\alpha_d}$. We integrate over $S(0,1)$, and we use the fact that the terms associated to $\abs{\alpha} = 1$, $\abs{\alpha} = 3$, and $\abs{\alpha} = 2$ except for the diagonal terms (trace of the Hessian matrix), all vanish by symmetry of $S(0,1)$. The proposed expansion of $- \Delta p~(x) $ is obtained by dividing the Taylor expansion by the measure of $S(x,\epsilon)$.
\end{proof}
\begin{prop}
Let $\Omega$ be a smooth domain, and $p$ harmonic in $\Omega$. Then the maximum of $p$ is attained on the boundary of $\Omega$.
\end{prop}
\begin{proof}
This is a straight consequence of Proposition~\ref{prop:mvp}.
\end{proof}
Another type of maximum principle, namely Hopf maximum principle, will be important to interpret the behaviour of macroscopic crowd motion models
\begin{prop}
Let $\Omega$ be a smooth domain, and $p$ a smooth scalar field such that $-\Delta p \geq 0$. At any point of $\partial \Omega$ which minimizes $p$, it holds that $\partial p / \partial n < 0$.
\end{prop}
\begin{coro}\label{coro:hopf}
Let $\Omega$ be a smooth domain, and $p$ such that $-\Delta p = f $, with $f$ a smooth nonnegative function that is non identically zero. Then $p$ is positive in $\Omega$, and $\partial p / \partial n <0$ on $\partial \Omega$.
\end{coro}
\begin{proof}
This a direct consequence of~\cite[Theorem~3.5]{GT}.
\end{proof}
In the discrete setting, the value of an harmonic function at a vertex is the mean value on the set of direct neighbors.
\begin{prop}\label{prop:mvpdisc}
Let $(V,E,c,\Gamma)$ be a resistive network, and $p$ a field that is harmonic on $\intV$. Then, for any $x\in \intV$, $p_x$ is an explicit convex combination of values at neighboring vertices:
\[
p_x = \sum_{y\sim x} \pi_{xy} p_y\virg \text{ with } \pi_{xy} = \frac { c_{xy}}{C_x} \virg C_x =
\sum_{y\sim x} c_{xy}.
\]
\end{prop}
\begin{proof}
This is a straight consequence of the harmonicity at $x$.
\end{proof}
Note that the counterpart of Proposition~\ref{prop:lapmean} for the continuous Laplacian is simply the definition of $L$, introduced as the difference between the vertex value and the mean value over neighbors.
\begin{rema}
The explicit formula given by the latter proposition is only valid for the strict neighborhood of a vertex. One can actually express a more general property, and express the value at some point $x$ as a convex combination of values taken in a set $\gamma$ of vertices that \emph{surrounds} $x$ in some sense. Yet, unlike in the continuous setting, this combination is not fully explicit, it relies on the \emph{harmonic measure} over $\gamma$, relative to $x$ (see Proposition~\ref{prop:MVPgamma}).
\end{rema}
\subsection{Stochastic setting}
The Dirichlet problem~\eqref{eq:dir} has the following stochastic interpretations (see e.g.~\cite{Doob84}).
\begin{prop}\label{prop:DirStoch}
Consider a smooth and bounded domain $\Omega$ in $\setR^d$, and $P$ a smooth function defined on the boudary $\Gamma = \partial \Omega$. For any $x\in \Omega$, let $w_x^t$ be the brownian motion from $x$. Let $\tau_x$ be the hitting time of $\Gamma = \partial \Omega$ by $w_x$, i.e.
\[
\tau _x = \inf \set{ t\virg w_x^t \in \Gamma}.
\]
Define $p$ as the expected value of $P$ at the hitting point:
\[
p(x) = \bbE\left(P\left(w_x^{\tau_x}\right)\right).
\]
Then $p$ is the solution to Dirichlet Problem~\eqref{eq:dir}.
\end{prop}
\begin{proof}
See~\cite[Chapter~IX]{Doob84}, or~\cite[Chapter~4]{karatzas91}.
\end{proof}
Another property involves the so-called \emph{harmonic measure} relative to a point.
\begin{prop}\label{prop:mesharm}
Consider a smooth and bounded domain $\Omega$ in $\setR^d$, and $x\in \Omega$. Let $w_{x_0}(t)$ be the brownian motion from ${x_0}$ and $\tau_{x_0}$ the hitting time of $\Gamma = \partial \Omega$ by $w_{x_0}$ like in the previous proposition. Then $Y_{x_0} = w_{x_0}^{\tau_{x_0}} \in \Gamma$ is a random variable valued in $\Gamma$, the law of which is the measure $\mu_{x_0}$ on $\Gamma$. This measure admits a density with respect to the Lebesgue measure on $\Gamma$, which is $-\partial q / \partial n$, where $q$ is the solution to the homogeneous Dirichlet problem with singular right-hand side
\begin{equation}\label{eq:dirdirac}
\left |
\begin{array}{rcll}
-\Delta q & = &\delta_{x_0} & \text{ in } \; \Omega\\
q &=& 0 & \text{ on } \Gamma.
\end{array}
\right.
\end{equation}
In other words, for any measurable subset $A\subset \Gamma$
\[
\bbP(Y_x \in A) = -\int_A \frac {\partial q}{\partial n}.
\]
The mesure $\mu_{x_0}$ with density $-\partial q/\partial n$ is called the harmonic measure relative to $x_0$.
\end{prop}
\begin{proof}
See~\cite{joensuu}.
\end{proof}
As a direct consequence of the two previous propositions, one can express the solution to the Dirichlet problem
\begin{prop}\label{prop:pmuy}
Let $\Omega$ be a smooth and bounded domain, and $p$ the solution to Dirichlet problem~\eqref{eq:dir} associated to the boundary condition $P$. Then, at any $x\in \Omega$, $p$ expresses
\[
p(x) = \int _{\partial \Omega} P(y) \, d\mu_x (y)
\]
where $\mu_x$ is the harmonic measure relative to $x$.
\end{prop}
This property is the conceptual core of Monte Carlo methods to solve non-homogeneous Dirichlet problems, we refer to the pioneering paper~\cite{montecarlo56} for more details.
%\medskip
Proposition~\ref{prop:DirStoch} admits a discrete counterpart, where the brownian motion is replaced by the random walk canonically associated to the resistive network (see e.g.~\cite{LyonsPeres17}).
\begin{prop}\label{prop:DirStochdisc}
Let $(V,E,c,\Gamma)$ be a connected resistive network, with boundary $\Gamma\neq \emptyset$, and $P\in \setR^\Gamma$ a collection of boundary values. We consider the Markov process on $V$ defined by $x\to y$ transition probabilities
\[
\pi_{xy} = \frac {c_{xy}}{C_x}\virg C_x = \sum_{y\,\sim\,x} c_{xy},
\]
with $c_{xy} = 0$ whenever $(x,y) \notin E$. For any ${x_0}\in V$, we denote by $(w_{x_0}^k)_k$ the random walk from ${x_0}$, according to the transition probabilities indicated above. We define the hitting time $k_{x_0}$ as
\[
k_{x_0} = \inf \set{ k\virg w_{x_0}^k \in \Gamma}
\]
and $Y_{x_0} = w_{x_0}^{k_{x_0}}$ the hitting point, that is a random variable valued in $\Gamma$. We finally define $p\in \setR^V$ by
\[
p_{x_0} = \bbE\left(P_{Y_{x_0}}\right).
\]
Then $p$ is the solution to the discrete Dirichlet problem~\eqref{eq:dirdisc}.
\end{prop}
\begin{proof}
If $x_0\in \Gamma$, then $ k_{x_0} = 0$ and $p_{x_0} = P_{x_0}$, so that $p$ verifies the boundary condition. If $x_0\notin \Gamma$, the property is a straight consequence of the law of total expectation. Indeed,
\[
p_{x_0} = \sum_{y\,\sim\,x_0} \pi_{ {x_0} y } \bbE\left(P_{Y_{y}}\right) =\sum_{y\,\sim\,x_0} \pi_{ {x_0} y } p_y =
\sum_{y\,\sim\,x_0} \frac {c_{ {x_0} y }}{C_ {x_0} } p_y,
\]
so that
\[
\sum_{y\,\sim\,x_0} {c_{ {x_0} y }}\left(p_{x_0} - p_y\right) = 0,
\]
i.e. $p$ is harmonic, which ends the proof.
\end{proof}
Proposition~\ref{prop:mesharm} also admits a discrete counterpart, which characterizes the discrete harmonic measure.
\begin{prop}[{Harmonic measure in the discrete setting}]\label{prop:mesharmdisc}
Let $(V,E,c,\Gamma)$ be a connected resistive network, with boundary $\Gamma\neq \emptyset$. Like in the previous proposition, for any $x_0\in \intV$, we define $Y_{x_0} = w_{x_0}^{k_{x_0}}$ as the hitting point in $\Gamma$ of the random walk from $x_0$. For any $y \in \Gamma$, the probability that $Y_{x_0} $ is $y_0$ writes
\[
\bbP\left(Y_{x_0} = y_0\right) = \p u(y_0) = - \p c \ps p (y_0),
\]
where $p$ is the solution to the Poisson problem
\begin{equation}\label{eq:dirdelta}
\left |
\begin{array}{rcll}
(\p c \ps p) _x& = &\delta_{x_0,x} & \forall x\in \intV\\
p &=& 0 & \text{ on } \Gamma.
\end{array}
\right.
\end{equation}
where $\delta_{x_0,x}$ is the Kronecker delta. The field $- \p c \ps p $ on $\Gamma$ is the harmonic measure with respect to~$x_0$.
\end{prop}
\begin{proof}
The proof relies on an auxiliary Laplace problem with Dirichlet boundary conditions. Let $q$ be the solution to
\begin{equation}
\left |
\begin{array}{rcll}
\p c \ps q (x)& = &0& \forall x\in \intV\vseq\\
q _y&=& \delta_{y,y_0} & \forall y \in \Gamma.
\end{array}
\right.
\end{equation}
Since $q$ is harmonic on $\intV$, and $p$ vanishes on $\Gamma$, it holds that
\[
\sum_{x\,\in\,V} p_x \left(\p c \ps q\right)_x = 0.
\]
Since the Laplacian $\p c\ps$ is selfadjoint, the latter expression also writes
\[
\sum_{x\,\in\,V} p_x \left(\p c \ps q\right)_x =
\sum_{x\,\in\,V} q_x \left(\p c \ps p\right)_x = q_{x_0} + \left(\p c \ps p\right)_{y_0},
\]
so that $q_{x_0} = - (\p c \ps p)_{y_0}$. Now, thanks to Proposition~\ref{prop:DirStochdisc}, $q_{x_0}$ is the probability that $Y_{x_0} $ is $y_0$, which ends the proof of Proposition~\ref{prop:mesharmdisc}.
\end{proof}
\begin{rema}
In a euclidean domain, if one extends the pressure field of Proposition~\ref{prop:mesharm} by $0$, $-\partial p / \partial n$ expresses mass unbalance on the boundary $\Gamma$. It is also $u\cdot n$, with $u = - \nabla p$. It identifies to the singular part of $-\Delta p $, which is a single layer distribution on $\Gamma$. In the discrete setting, there is no such notion as ``normal derivative'', but the mass unbalance at a boundary vertex $y_0$ is $(\p u)_{y_0}$ (that is the outflow through $y_0$), with $u = - c \ps p$, so that the mass unbalance writes $- (\p c \ps p)_{y_0}$, which is therefore the discrete counterpart of $-\partial p / \partial n$ at $y_0$, as expressed by the previous proposition.
\end{rema}
The notion of harmonic measure makes it possible to establish a generalization of the Mean Value Property in the discrete setting. The following proposition is a straight discrete counterpart of Proposition~\ref{prop:pmuy}.
\begin{prop}\label{prop:MVPgamma}
Let $(V,E,c,\Gamma)$ be a connected resistive network, with boundary $\Gamma\neq \emptyset$, and $p\in \setR^V$ a field that is harmonic in $\intV$. Let $x_0\in \intV$ be given, and $\gamma \subset \intV$ a set of interior vertices which \emph{isolates} $x_0$ from $\Gamma$, in the sense that all the paths from $x_0$ to any vertex of\ $\Gamma$ necessarily meet $\gamma$. We denote by $V_{x_0}$ the set of vertices which can be connected to $x_0$ without meeting $\gamma$. Then $p_x$ writes as a convex combination of values taken by $p$ on $\gamma$:
\begin{equation}\label{eq:pxeq}
p_x = \sum _{y\,\in\,\gamma } \theta _y p_y\virg \theta _y \geq 0\quad \forall y \in \gamma\virg
\sum _{y\,\in\,\gamma } \theta _y = 1.
\end{equation}
The field of weights $(\theta_y)_\gamma$ is the harmonic measure on $\gamma $ relative to $x_0$, in the sense of Proposition~\ref{prop:mesharmdisc}, i.e.
\[
\theta_y = - \left(\p c \ps q\right) _y,
\]
where $q$ is the solution to the Dirichlet problem on $V_{x_0}\cup \gamma$
\begin{equation}\label{eq:dirdeltagamma}
\left |
\begin{array}{rcll}
\left(\p c \ps q\right) _x& = &\delta_{x_0,x} & \forall x\in \intV_{x_0}\vseq\\
q_y &=& p_y & \forall y \in \gamma.
\end{array}
\right.
\end{equation}
\end{prop}
\begin{proof}
Let us consider the subnetwork $(V',E',c',\Gamma')$, where $V' = V_{x_0} \cup \gamma$, $E'$ contains the edges of $E$ that connect two vertices in $V'$, $c'$ is simply the restriction of $c$ to $E'$, and $\Gamma' = \gamma$. We apply Proposition~\ref{prop:mesharmdisc} to the subnetwork $(V',E',\Gamma')$. It ensures that $\theta _y = - (\p c \ps q) _y$ is the probability that a random walk from $x_0$ meets $\gamma $ at $y$, where $q$ solves~\eqref{eq:dirdeltagamma}.
%~\ref{prop:mesharmdisc}.
Now since $\gamma $ isolates $x_0$ from $\Gamma$, any random walk starting from $x_0$ and ending at first hit with $\gamma$ does not meet $\Gamma$, thus it only visits vertices in $\intV$, at which the field $p$ is harmonic, which makes it possible to use Proposition~\ref{prop:DirStochdisc} applied to the new network $(V',E',\Gamma')$, with $\Gamma' = \gamma$. As a consequence, $p_{x_0}$ is the expected value of $p$ at the first hit with $\gamma$, which is expressed by~\eqref{eq:pxeq}.
\end{proof}
\section{Macroscopic and microscopic congestion models}\label{sec:crowd}
The approach which we present in this section was initially motivated by the handling of collisions in fluid-particle simulations (see e.g.~\cite{GloMaury97,Glo99,MauryGran06}). We shall present a common strategy to handle congestion at microscopic and macroscopic scales. This strategy can be carried out for second order in time evolution problems, see~\cite{MauryGran06} for microscopic granular flows, and~\cite{MauryPreux17} for the macroscopic version. We also refer to~\cite{degond18}, where the authors propose a numerical scheme to solve Euler equations with a singular pressure meant to account for congestion, or~\cite{perrin21} for the asymptotic analysis of a similar system. We shall present here the overdamped versions of these inertial problems, with a hard handling of the maximal density constraint, in the context of crowd motion modeling.
We describe in this section the macroscopic and the microscopic settings, and we show how they both lead to a Poisson-like problem. The next section proposes a deeper look on the discrepancy between both settings, and a detailed analysis of the pathologies of the discrete Laplacian.
\subsubsection*{Macroscopic setting}
The macroscopic crowd motion model proposed in~\cite{MRS10} reads as follows. Given a domain $\O$ (which covers the zone possibly occupied by the crowd), the crowd is represented by a density $\rho = \rho(x,t)$ subject to remain below a certain maximal value, which we set to $1$. The set of feasible density is thus
\[
K = \set{ \rho \in L^1(\O)\virg 0\leq \rho \leq 1\text{ a.e.}}.
\]
We denote by $U = U(x)$ the desired velocity, that is given: $U(x)$ is the velocity that a person at $x$ would take, would they be alone. We write that $\rho$ is transported by the actual velocity $u$ :
\begin{equation}\label{eq:transport}
\partial _t \rho + \nabla \cdot (\rho u) = 0,
\end{equation}
where $u$ is is the closest to $U$ among all those velocities which do not lead to a violation of the constraint, i.e. which are expansive in the zone which are already saturated. More formally, it amounts to prescribe a nonnegative divergence in the saturated zone. This constraint can be expressed in a dual way by subjecting $u$ to belong to
\[
C_\rho = \set{ v \in L^2(\O)\virg \int_\O v \cdot \nabla q \leq 0 \quad \forall q \in H_\rho^1(\O)\virg q \geq 0 \text{ a.e.} }
\]
were
\[
H_\rho^1(\O) = \set{ q \in H^1(\O)\virg q (1-\rho) = 0 \text{ a.e.} }.
\]
Note that the condition $q (1-\rho)=0$ contrains $q$ to vanish in zones where $\rho<1$, i.e. where the constraint is \emph{not} saturated. We may now write the full model by defining the actual velocity as
\begin{equation}\label{eq:uPU}
u = P_{C_\rho} (U),
\end{equation}
where the projection is defined in the $L^2$ sense.
\begin{rema}
Model~\eqref{eq:uPU}, like its microscopic counterpart~\eqref{eq:evolmic} below, is not meant to account for sophisticated individual tendencies, it merely encodes a very basic strategy in terms of behavior: individual strive to achieve their desired velocity, regardless of others. The resulting interactions are physical in nature (contact forces), and they may induce nonlocal effects over connected saturated zones.
\end{rema}
Boundary conditions can be accounted for through the definition of the pressure field. As an archetypal example, consider the situation represented in Figure~\ref{fig:room} where $\O$ is a room with boundary $\Gamma$, union of $ \Gamma_{w}$ (walls) and $ \Gamma_{out}$ (exit). The saturated zone $[\rho=1]$ corresponds to the grey domain, which we denote by $\Omega$. In this situation, we prescribe $q = 0$ on $\Gamma_{out}$, whereas $q$ is ``free'' on $\Gamma_w$. Note that, due to the constraint $q(1-\rho) = 0$, the pressure vanishes on the upstream side of the saturated zone $\Omega$. Let us now comment on the definition of $C_\rho$, in this particular situation. Suppose $\Omega = [\rho=1]$ is a smooth domain, and suppose $u\in C_\rho $ is smooth, considering smooth and nonnegative test functions $q$ in $H_\rho^1$, compactly supported in $\OO$, the dual condition yields
\[
\int_\O u \cdot \nabla q = - \int_\Omega q\nabla \cdot u \leq 0,
\]
so that $\nabla \cdot u\geq 0$ a.e.\ in $\OO$. Now consider test functions which are no longer compactly supported in $\Omega$, which can take positive values on $\Gamma_w$. The same integration by part yields a boundary term, so that we obtain
\[
\int _{\Gamma_{w}} q u\cdot n \leq 0,
\]
which imposes $u\cdot n \leq 0$ on $\Gamma_w$: people may not exit through the wall.
\begin{figure}[tbp]
\includegraphics[width=0.65\linewidth]{room.pdf}
\caption{Evacuation from a room.}\label{fig:room}
\end{figure}
We refer to~\cite{MRS10} for a full analysis of the evolution problem~\eqref{eq:transport}\eqref{eq:uPU}. This analysis is delicate due to the nonlinearity and lack of smoothness of $\rho \longmapsto u = P_{C_\rho} U$, which rules out standard PDE approaches. The natural framework is the Wasserstein space, built upon a distance between densities that is based on optimal transportation. We shall disregard here these aspects, and rather focus on the projection step~\eqref{eq:uPU}, which defines $u$. In a paradoxical way, we shall see that this step is much more delicate and rich at the microscopic scale than at the macroscopic one. Indeed, implementing~\eqref{eq:uPU} amounts to minimize the functional
\[
v \longmapsto J(v) = \frac{1}{2}\int _{\O} \abs{ v - U }^2,
\]
on the feasible set $C_\rho$. Since $C_\rho$ is closed and convex in $L^2(\O)$, and since $J$ is continuous and strongly convex, the projection problem admits a unique solution $u$.
\subsubsection*{Saddle-point formulation and Poisson problem}
Let us introduce
\[
B\; : \; v \in L^2(\O) \longmapsto -\nabla \cdot v \in H^{-1}_\rho = \left (H^1_\rho\right) '
\]
defined by
\[
\left\langle-\nabla \cdot v,q\right\rangle = \int _{\O} v\cdot \nabla q.
\]
The projection problem~\eqref{eq:uPU} can be formulated as a saddle-point problem (see e.g.~\cite{VarIneqGlo}) for the Lagrangian
\begin{equation}\label{eq:Ldef}
L(v,q) =J(v) + \left\langle Bv\,\middle|\,q\right\rangle = \frac{1}{2}\int _{\O} \abs{ v - U }^2 + \int _{\O} v\cdot \nabla q.
\end{equation}
\begin{prop}\label{prop:sp}
The Lagrangian defined by~\eqref{eq:Ldef} admits a unique saddle point $(u,p) $ in $ V \times \Lambda_+$, with
\[
V = L^2(\O)\virg \Lambda_+ =\set{ q\in H_\rho^1(\O) \virg q \geq 0\quad \text{ a.e. } }.
\]
This saddle point is characterized by $u \in L^2(\O)$, $p\in H^1_\rho(\O)$
\begin{equation}\label{eq:uBU}
\left |
\begin{array}{rlll}
u + \nabla p& = &U& \text{ in } \O\vseq \\
\qquad -\nabla \cdot u &\leq & 0 & \text{ on }\, [\rho = 1] \vseq \\
p &\geq & 0 & a.e. \text{ in }\, \O\vseq \\
\dsp \int_{\O} u \cdot \nabla p&=&0 &.
\end{array}
\right.
\end{equation}
where the second equation is meant in a weak sense, i.e.
\[
\int _{\O} u\cdot \nabla q\leq 0
\]
for any non-negative test function which vanishes outside $ [\rho = 1]$. The primal part of this saddle-point is the projection of $U$ on $C_\rho$.
\end{prop}
\begin{proof}
This problem fits into the abstract framework of Proposition~\ref{prop:absspp} in Section~\ref{sec:abs}, in the infinite dimensional situation. Well-posedness comes from that fact that $B$ is surjective, which is a straight consequence of Poincaré inequality (see~\cite[Proposition~3]{BourdinMaury2020} for details).
\end{proof}
Now consider a situation where the saturated zone is a smooth domain (like in the evacuation situation represented in Figure~\ref{fig:room}). As detailed in~\cite{BourdinMaury2020}, if the desired velocity field in concentrating, i.e. $\nabla \cdot U < 0$, then $p$ cannot vanish on an open subset on $\Omega$ (it would rule out the constraint in this subdomain). Thus, the complementarity slackness condition ($4^{\rm th}$ line in~\eqref{eq:uBU}) imposes that the constraint is saturated, i.e. $ \nabla \cdot u = 0$, in $\Omega$. As a consequence, one may eliminate the velocity to obtain a Poisson problem on the pressure with mixed boundary conditions:
\begin{equation}\label{eq:dirneu}
\left |
\begin{array}{rclll}
-\Delta p & = &-\nabla \cdot U & \text{ in } &\Omega\vseq\\
p &=& 0 & \text{ on } &\Gamma_{out} \cup \Gamma_{up} \vseq \\
\dsp \frac {\partial p }{\partial n } &= &0 & \text{ on } &\Gamma_{w}.
\end{array}
\right.
\end{equation}
As we shall see in the next section, this Poisson problem is central in the behavior of this congestion model.
\subsubsection*{Microscopic setting}
The microscopic counterpart of the previous approach was introduced in~\cite{MVm2an}. It is Lagrangian in nature, based on the position of $N$ individuals $x = (x_1,\dots, x_N)\in \setR^{2N}$. The individuals are identified to rigid discs of radius $r$, they are subject to a non-overlapping constraint, i.e. $g$ must remain in the feasible set
\[
K = \set{x\in \setR^{2N}\virg D_{ij}(x) = \abs{x_j - x_i} - 2r \geq 0 \quad \forall i\neq j }.
\]
We consider a collection of desired velocities $U = (U_1,\,\dots,\,U_N) \in \setR^{2N}$. The constraint on the divergence at the macroscopic level becomes here a non-overlapping constraint: Feasible velocities are those fields which do not decrease distances that are already $0$. For any $x\in K$, the feasible set for velocities writes
\[
C_x = \set{ v\in \setR^{2N}\virg D_{ij}(x) = 0 \Longrightarrow e_{ij}\cdot \left(v_i - v_j\right) \leq 0 },
\]
where $e_{ij} = (x_j-x_i) / \abs{x_j-x_i}$ is the unit vector pointing from $x_i$ to $x_j$. The evolution problem, which is the microscopic counterpart of~\eqref{eq:transport}\eqref{eq:uPU}, writes
\begin{equation}\label{eq:evolmic}
\left |
\begin{array}{rcl}
\dsp \frac {dx}{dt}& = &u\vseq\\
u &=& P_{C_x} U.
\end{array}
\right.
\end{equation}
Note that, in this purely Lagrangian setting, the Eulerian transport equation~\eqref{eq:transport} is replaced here by a simple definition of the velocity as the time derivative of the position vector.
Like in the macroscopic context, the set of feasible velocities is a closed convex cone, which ensures existence and uniqueness of a projection. The analysis of the time evolution problem, which is not the topic of the present paper, is detailed in~\cite{MVm2an}, it relies on the so-called prox-regularity of the position feasible set $K$, which essentially says that the projection on this set is locally well-defined. We shall focus here on the projection step which defines $u$. Let us first introduce a matrix formulation of the constraints: the set of feasible velocities can be written
\[
C_x = \set{ v\virg Bv \leq 0}
\]
where each row of $B$ stands for an active contact, i.e. a couple $(i,j)$ such that $D_{ij} = 0$, and thus writes
\[
B_{ij} = \left(\underset{\substack{\uparrow\\1}}{0},\underset{\substack{\uparrow\\2}}{0},\,\dots,\,\underset{\substack{\uparrow\\i}}{e_{ij}}, 0,\,\dots,\,0, \underset{\substack{\uparrow\\j}}{-e_{ij}}, 0,\,\dots,\, \underset{\substack{\uparrow\\N}}{0}\right)\in \setR^{2N}.
\]
The matrix $B$ belongs to $\clM_{N_c\times 2N} (\setR)$, where $N_c$ is the current number of contacts. This number may vary in time, but we shall here focus on the instantaneous problem.
The constrained minimisation problem admits a saddle-point formulation for the Lagrangian
\begin{equation}\label{eq:Ldefdisc}
L\left(v,q\right) =J(v) + \left\langle Bv\,\middle|\,q\right\rangle = \frac 1 2 \sum_{i=1}^N \abs{ v_i - U_i }^2 + \sum_{i\,<\,j,D_{ij} = 0} p_{ij}\times e_{ij}\cdot \left(v_i - v_j\right),
\end{equation}
as expressed in the following proposition (discrete counterpart of Proposition~\ref{prop:sp}).
\begin{prop}\label{prop:spdisc}
The Lagrangian defined by~\eqref{eq:Ldefdisc} admits a saddle point $(u,p) $ in $ V \times \Lambda_+$, with
\[
V = \setR^{2N}\virg \Lambda_+ =\setR_+^{N_c}.
\]
This saddle point is characterized by
\begin{equation}\label{eq:uBUdisc}
\left |
\begin{array}{rcl}
u + B^\star p& = &U \\
B u &\leq & 0 \vseq \\
p &\geq & 0 \vseq \\
\left\langle p\,|\,Bu\right\rangle &=& 0.
\end{array}
\right.
\end{equation}
The primal part of this saddle-point is the projection of $U$ on $C_x$.
\end{prop}
\begin{proof}
This elementary property relies on classical arguments: like in the macroscopic setting (Proposition~\ref{prop:sp}), this problem fits into the abstract framework of Proposition~\ref{prop:absspp}, in the case where the number of constraints is finite. Note that no assumption is made on $B$, in particular surjectivity is not required, so that uniqueness may be ruled out for the dual component $p$.
\end{proof}
\begin{figure}[tbp]
\includegraphics[width=0.95\linewidth]{row.pdf}
\caption{One-dimensional cluster}\label{fig:row}
\end{figure}
All matrices take a simple form in the one-dimensional setting. Consider a cluster of $N$ discs in a row (see Figure~\ref{fig:row}). Matrix $B$, discrete opposite of the divergence, and $BB^\star$, write
\begin{equation}\label{eq:1dlap}
B =
\left (
\begin{array}{ccccc}
1& -1 & 0 & \dots&\dots \\
0 & 1 & -1 & \dots&\dots \\
0 & 0 & \ddots & \ddots&\dots \\
0 & 0 & \dots&1 & -1 \\
\end{array}
\right) \virg B B^\star =
\left(
\begin{array}{ccccccc}
2& -1& 0 &\cdot&\cdot& 0 \vseq\\
-1 & 2& -1 &0&\cdot& \cdot \vseq\\
\cdot& & \cdot &\cdot& \cdot& \cdot \vseq\\
\cdot& & &-1& 2& -1 \vseq\\
0& \cdot & \cdot &0& -1& 2
\end{array}
\right),
\end{equation}
which is the discrete Laplacian in dimension one.
\begin{rema}
The operator $B$ plays here the role of $\p$ for resistive networks (see Section~\ref{sec:savs} and Figure~\ref{fig:VVchain}), which is the discrete counterpart of the opposite of the divergence. Mass balance is replaced here by overlap estimate. Indeed, by definition, each line of $B$ expresses the first order term of the expansion of the distance, e.g
\[
D_{ij}(x + \epsilon U) = \underbrace{D_{ij}(x)}_{=0} + \epsilon (BU)_{ij} + o(\epsilon).
\]
Accordingly, $B^\star$ plays the role of $\ps$, which is a discrete gradient. For a resistive fluid network, a pressure field induces a collection of fluxes proportional to $\ps p$, which write $u_{xy} = c_{xy} (p_x - p_y)$ in each edge $(x,y)$. In the present context, a pressure field $p\in \setR^{N_c}$ induces a force field $B^\star p$ over the primal degree of freedom, which are the positions of the grains. Note that the fact that $B$ and $B^\star$ are mutually adjoint is not fortuitous like in the case of resistive networks, this property is inherent to the saddle-point formulation which defines $p$.
\end{rema}
\begin{rema}\label{rem:dir}
Note that this matrix corresponds to the Finite Difference discretization of a Laplace problem \emph{with Dirichlet boundary condition}. Since the grain $1$ (see Figure~\ref{fig:row}) has no neighbor on its left-hand side, it is as if a zero pressure were imposed at the grey circle on its left. In the Finite Difference context, these grey circles correspond to boundary points at which a zero value is prescribed, so that they do not correspond to unknowns, and they are in general removed from the matrix formulation. This will play an important role in what follows: we shall define Laplace-like operators on clusters of discs, and those operators will act on pressures which are defined on \emph{edges} of the contact network, which is the dual graph. The resulting discrete Laplacian corresponds to a situation where a zero pressure is prescribed on the border.
\end{rema}
\begin{figure}[tbp]
\includegraphics[width=0.95\linewidth]{primaldual.pdf}
\caption{Primal and dual graphs}\label{fig:primaldual}
\end{figure}
For discs in $\setR^2$, following the idea that $B$ is the discrete counterpart of the opposite of the divergence, it can be expected that the matrix $BB^\star$ can be interpreted as a discrete Laplacian. Figure~\ref{fig:primaldual} represents such a cluster of discs, the associated primal graph, and the associated dual graph (in red), the vertices of which are the points of application of the interactions forces. The operator $BB^\star$ acts on the \emph{dual graph} of the contact network. As illustrated in the figure, the dual graph is not planar as soon as some discs have more than $3$ neighbors, whereas the primal graph is planar. The thick red lines in the figure correspond to the \emph{stencil} associated to a ``vertex'' of the dual graph, that is the contact $(i,k)$. This dual vertex has $7$ neighbors, which are the other contacts that implicate either $i$ or $k$. One may imagine that some substance circulates among discs, with a flux between two discs given by a Fick-like law. This would lead to a standard Laplacian of the type $\p c\ps$, defined on the primal graph (in black in the figure), like in Section~\ref{sec:intro}. The maximal degree would be $5$ (grain $k$ has 5 neighbors), whereas the maximal degree of the dual graph is $7$ here. As we shall see, the operator $BB^\star$ which emerges from the saddle-point formulation, which is dual in some way to $\p c\ps$, presents very different characteristics.
Like in the macroscopic setting, one may deduce from Proposition~\ref{eq:uBUdisc} that the pressure $p$ verifies some sort of Poisson problem. Yet, an extra difficulty comes from the fact that the constraint $Bu \leq 0 $ is not necessarily saturated, even in the case of a concentrating velocity field $U$, i.e. such that $BU >0$, which is the first major difference with the macroscopic setting. Let us explain this first discrepancy in a very simple situation, illustrated by Figure~\ref{fig:coin} (right), which is a discrete version of the macroscopic situation represented in the left-hand side of the figure. The situation is represented in the upper right corner of the figure : 4 discs form a cluster, and the desired velocities push them toward the center. For any two discs in contacts, moving the discs along desired velocities would lead to a violation of the constraint, which writes $BU> 0$.
A macroscopic counterpart of this situation is represented on the left-hand part of the figure: the domain is an ellipse, and the desired velocity points toward the center, in such a way that all parts of the domain are compressed by $U$. It writes $-\nabla \cdot U > 0$. In the macroscopic situation, the pressure field is the solution to the Poisson problem $-\Delta p = - \nabla \cdot U>0$, so that $p>0$ in $\Omega$ by Corollary~\ref{coro:hopf}. In this situation, minimizing $\abs {v- U}$ under the constraint $-\nabla\cdot v \leq 0$ thus leads to a pressure field that is positive, so that minimizing over divergence-free fields would lead to the very same solution. In the microscopic setting the situation is fully different. Minimizing over $\ker B$ amounts to consider that grains are glued to each other (they actually can slip on each other, but distances of grains in contact have to remain equal to $0$). For symmetry reasons, the projection of $U$ over $\ker B$ is obviously $0$. In the saddle-point formulation, the interaction force between the two grains that are on the vertical axis is \emph{negative}, whereas the other pressures are positive (see the `+' and `-'signs on the figure), which is in full opposition to the macroscopic situation. As we shall see, it comes from the fact the Laplace operator that we are going to exhibit does not verify any maximum principle. If one now considers the inequality constraint $Bu\leq 0$, the solution will be different: the two grains on the vertical axis shall be pushed apart, due to the horizontal squeezing action on the grains on the horizontal axis (\emph{wedge} effect).
\begin{figure}[tbp]
\includegraphics[width=0.95\linewidth]{coin.pdf}
\caption{Wedge principle.}\label{fig:coin}
\end{figure}
The example that we just considered shows that supposing $BU>0$ (grains are pushed agains each other) does not guarantee that the constraint is saturated (i.e. $Bu = 0$). To obtain a Laplace operator by eliminating the velocity thus requires that we first remove from $B$ the rows which correspond to non-active constraints, i.e. such that the actual velocity fields creates a \emph{positive} relative velocity: grains are strictly pulled apart each other.
\begin{prop}\label{prop:oB}
For a given desired velocity field $U$, $(u,p)$ being a solution to the saddle-point formulation~\eqref{eq:uBUdisc}, we denote by $\oB$ the matrix obtained from $B$ by removing all the rows associated to a non active constraint, i.e. for couple $(i,j)$ such that $e_{ij} \cdot (v_i - v_j) < 0$. Then $p$ verifies the discrete Poisson problem
\begin{equation}\label{eq:Poissondisc}
\oB \oB^\star \op = \oB U,
\end{equation}
where $\op$ is obtained from $p$ by removing all the entries that correspond to removed constraints.
\end{prop}
\begin{proof}
By construction, we have that $\oB u =0$ (all the rows with a strict inequality have been removed). Now, thanks to the complementarity condition, for all those raws $(i,j)$ which have been removed, it holds that $p_{ij} = 0$, so that $\oB^\star \op = B^\star p$. We thus obtain the discrete Darcy problem
\begin{equation}\label{eq:darcydiscB}
\left |
\begin{array}{lcl}
u + \oB^\star \op& = &U \\
\qquad \oB u &=& 0
\end{array}
\right.
\end{equation}
which yields, after elimination of the velocity, to the discrete Poisson equation~\eqref{eq:Poissondisc}.
\end{proof}
Let us insist on the fact that the discrete Laplacian $\oB \oB^\star$ does not only depend on the configuration. For a given $x$ in $K$ (collection of non-overlapping discs), $BB^\star$ is the basic discrete Laplacian but, given a compressing field $U$ (i.e. such that $BU>0$), the induced pressure will verify a Poisson problem $\oB \oB^\star\op = \oB U$, where $\oB$ is obtained by removing some rows from $B$, in a way \emph{that depends on} $U$. In other words, there is no canonical Laplacian associated to the dual graph, but a (finite) collection of Laplacian operators which are defined on subgraphs of the dual graph.
\begin{rema}\label{rem:coins}
Let us give a very concrete illustration of the previous considerations. Consider a collection of many coins on a table, with no overlapping but possibly touching. If one exerts centripetal forces from the outside of the cluster, the induced motion may push some coins away from each other, and the collection of couples that are split depends on the way the squeezing force is prescribed.
\end{rema}
To better understand the origin of this pathology, let us describe in more detail the matrix $B$, in the form of a proposition which expresses a negative property.
\begin{prop}
Let $x\in K$ represent a collection of $N$ non-overlapping discs. If there exists at least a disc $i$ which is in contact with $j$ and $k$ such that $e_{ij}$ and $e_{jk}$ form an \emph{acute} angle, then the matrix $BB^\star$ is \emph{not} a $Z$--matrix, i.e. it contains extra-diagonal positive entries.
\end{prop}
\begin{proof}
This comes from the expression of the matrix $BB^\star$. The simple example of such a situation is represented in Figure~\ref{fig:eijeik}, which is an extraction of a $3$-cluster from the configuration represented in Figure~\ref{fig:primaldual}. Let us compute the $2\times 2$ contribution to $BB^\star$ that comes from the interaction between $(i,j)$ and $(i,k)$ (red dots in the figure). The corresponding matrices are
\[
B_{(i,j),(i,k)} =
\left (
\begin{array}{ccc}
e_{ij}& -e_{ij} &0 \\
e_{ik}& 0 &-e_{ik}
\end{array}
\right) \in\clM_{2,6}(\setR)\virg
\]
so that the contribution of this interaction to the global Laplacian is
\[
B_{(i,j),(i,k)} B^\star_{(i,j),(i,k)} =
\left (
\begin{array}{ccc}
2 \abs{e_{ij}}&e_{ij}\cdot e_{ik} \\
e_{ij}\cdot e_{ik}& 2 \abs{e_{ik}}
\end{array}
\right) =
\left (
\begin{array}{ccc}
2 &\dsp1/ 2 \vseq\\
\dsp 1/ 2 & 2
\end{array}
\right)
\in\clM_{2,2}(\setR).
\]
By gathering all those contributions, the extra diagonal term which corresponds to the $(i,j),(i,k)$ interaction will remain unchanged, thus positive. As a consequence, the resulting matrix $BB^\star$ has \emph{positive} extra-diagonal entries.
\end{proof}
\begin{figure}[tbp]
\includegraphics[width=0.5\linewidth]{eijeik.pdf}
\caption{Acute angle condition}\label{fig:eijeik}
\end{figure}
In the continuous setting, a smooth function defined in a smooth connected domain is constant if and only if its gradient is identically $0$. In the discrete setting, where $B^\star$ plays the role of the gradient, this equivalence is no longer true, and this negative property is not due to boundary effects.
If one considers the situation represented in Figure~\ref{fig:primaldual}, the disc $k$ is in contact with $5$ discs, and the contact points are not distributed in a symmetric way. As a consequence, a unit pressure field induces a non-zero force on grain $k$, i.e. $B^\star 1 \neq 0$.
Conversely, in highly congested situations, the discrete gradient $B^\star$ may have a rich kernel. Consider for example a triangular lattice of discs with a common radius. If the number of discs is $N$, there are $2N$ primal degrees of freedom, whereas the number of contacts $N_c$ scales like $3N$. The kernel of $B\in \clM_{N_c, 2N}$ has therefore a dimension that scales like $N$. As a direct consequence, the Poisson problem
\[
BB^\star p = b
\]
which is the discrete counterpart of a standard Poisson problem with Dirichlet boundary conditions (see Remark~\ref{rem:dir}), is highly \emph{ill-posed} in terms of uniqueness. Let us nevertheless add that, even in situations where the problem is very degenerate, the unilateral homogenous problem is well-posed: the problem $ BB^\star p = 0$, with $p\in \setR_+ ^{N_c}$ admits $0$ as unique solution. This is a straight consequence of Hahn--Banach theorem, as detailed in~\cite{MauryGran06}. A simple illustration of this property is the following: consider a static collection of coins on a table, possibly in contact with each other. Assuming that coins are not glued to each other, the question is: is is possible to determine interaction forces ? The answer is yes. Consider an extreme coin (the center of which is not in the convex hull of other centers). Since forces are repulsive, if the forces acting upon this coin were not zero, they would push the coin away from the cluster. So they are zero, the coin can be removed without perturbing the equilibrium, and the process can be further implemented until exhaustion. Note that, in the same situation, if some coins are glued together (i.e. pressures can be positive or negative), it is impossible to determine interaction forces from the sole equilibrium property. It can also be proved that the solution set of the Poisson problem with positivity constraint is bounded. Indeed (see again~\cite[Lemma~1]{MauryGran06}), it holds that $\ker B^ \star \cap \setR_+^{N_c}$ is bounded, so that any Poisson problem $BB ^\star p = b$ set in $\ker B^ \star \cap \setR_+^{N_c}$ admits a bounded (possibly empty) solution set.
Before exploring the effect of the (negative) properties of $BB^\star$ upon the behavior of the evolution model~\eqref{eq:evolmic}, let us summarize these properties. Its form provides it with a variational setting. Indeed, as illustrated in Figure~\ref{fig:VV}, the pressure that is a solution to the saddle-point formulation~\eqref{eq:uBUdisc}, minimizes $\abs {U- B^\star q}^2$. As for the other properties listed in Section~\ref{sec:propslap}, most of them are \emph{not} verified by $BB^\star$. Since some extra diagonal terms can be positive, there is no mean value properties (Proposition~\ref{prop:mvpdisc} is not true), no stochastic interpretation (Propositions~\ref{prop:DirStochdisc}, \ref{prop:mesharmdisc}, and~\ref{prop:MVPgamma} do not make any sense, since there is no random walk canonically associated to $BB^\star$).
\section{Capacity drop and Faster is Slower effect at the macroscopic and macroscopic scales}
We explain in this section how the ``nice'' properties of the macroscopic (standard) Laplace operator prevent the corresponding model to reproduce some observed effect and, on the other side, how the pathologies of the microscopic operator make it possible to reproduce these effects.
\subsubsection*{Capacity drop}
In the context of evacuation, the \emph{Capacity Drop} (CD) effect designs a reduction of the flux though an exit door or a bottleneck when the number of evacuees upstream increases, see e.g.~\cite{CapaD05}. Let us first show that the macroscopic model~\eqref{eq:transport}~\eqref{eq:uPU} does \emph{not} reproduce the CD effect.
\begin{prop}\label{prop:FiF}
Let us consider the situation represented in Figure~\ref{fig:room}: pedestrians are represented at some time by a density $\rho$ upstream the exit $\Gamma_{out}$, the saturated zone ($\rho=1$) is a smooth domain $\Omega$, and the desired velocity field points to the door, in such a way that $\nabla \cdot U \leq 0$, and this quantity is not identically $0$ on $\Omega$. Then, according to the model~\eqref{eq:transport}~\eqref{eq:uPU}, people actually exit \emph{faster} than they would if there were no congestion.
\end{prop}
\begin{proof}
This is a straight consequence of the Maximum Principle. We give here a formal proof of this property, and we refer to~\cite{BourdinMaury2020} for a more detailed proof. Since $U$ is concentrating, the contraint $\nabla \cdot u\geq 0$ is saturated over $\Omega$, so that $\nabla \cdot u$ is identically $0$. The unilateral Darcy problem~\eqref{eq:uBU} thus becomes a standard Darcy problem. Elimination of the velocity yields the Poisson problem
\[
- \Delta p = - \nabla \cdot U,
\]
with homogeneous Dirichlet boundary conditions on $\Gamma_{out}$ (the exit), $\Gamma_{up}$ (the upstream front), and homogeneous Neumann condition on $\Gamma_{w}$ (walls). The field $p$ is such that $-\Delta p \geq 0$, and not identically $0$ in $\omega$, so that $p>0$ in $\Omega$. As a consequence, for any $x\in \Gamma_{out}$ the value of $p$ is strictly less that $p(y)$ for any $y\in \Omega$ and, by Corollary~\ref{coro:hopf}, $\partial p/\partial n (x) <0$. As a consequence
\[
u\cdot n = (U - \nabla p) \cdot n = U\cdot n - \frac {\partial p}{\partial n} (x) >U\cdot n,
\]
which exactly states that people exit faster than they would if they were alone (i.e. if they adopted their desired velocity).
\end{proof}
The previous property shows that the macroscopic model in its basic form is not relevant to recover the Capacity Drop phenomenon. Note that some attempts have been proposed to incorporate this phenomenon in macroscopic models. In~\cite{Andreianov14,Andreianov16}, the authors prescribe a limitation of the flux at the exit, in a way that non-locally depends on the number of people accumulated upstream the exit, which makes it possible to enforce the CD phenomenon. Yet, up to our knowledge, no macroscopic model is currently able to natively reproduce this phenomenon, in a way that would explain it.
We shall investigate now how the microscopic setting, because of the very defects of the discrete Laplacian $BB^ \star$, reproduce the Capacity Drop phenomenon. The discrete model actually exhibits an extreme version of the Capacity Drop phenomenon, that is the spontaneous emergence of static jams.
\begin{figure}[tbp]
\includegraphics[width=0.95\linewidth]{arche.pdf}
\caption{Principle of a static jam}\label{fig:arche}
\end{figure}
\begin{defi}[{Static jam}]
\label{def:staticjam}
In an evacuation situation like the one represented in Figure~\ref{fig:room}, a static jam is an equilibrium point of the evolution problem~\eqref{eq:evolmic}, i.e. a configuration $x\in K$ such that $P_{C_x} U = 0$. By Proposition~\ref{eq:uBUdisc}, it can also be defined as any configuration $x\in K$ such that there exists a pressure field $p\in \setR_+^{N_c}$ which balances $U$, in the sense that $B^\star p = U$.
\end{defi}
We may express in an informal way the ability of the microscopic model to reproduce the Capacity Drop phenomenon:
\emph{In an evacuation situation like the one represented in Figure~\ref{fig:room}, with desired velocities pointing to the outside of the room through the exit door, the microscopic evolution model~\eqref{eq:evolmic} leads in some situations to static jams in the sense of Definition~\ref{def:staticjam}, even when the width of the exit is significantly larger that the size of an individual.}
The mechanism of jams is illustrated in Figure~\ref{fig:arche}. The bottom corresponds to the macroscopic situation, it illustrates Proposition~\ref{prop:FiF}: the red arrows correspond to the variations of velocities due to congestion, i.e. $- \nabla p$, which points \emph{outward} the domain, meaning that people exit \emph{faster} than they would if they were alone (capacity rise phenomenon). This property is a direct consequence of the Maximum Principle for $-\Delta$, as detailed in Proposition~\ref{prop:FiF}. The top part of the figure corresponds to the microscopic situation. According to Remark~\ref{rem:dir}, one can consider that the grey dots, located at the zone of grains facing the exit, correspond to points where a zero pressure is prescribed. Now consider positive pressures between grain in contact, the corresponding correction on the velocity, that is $-B^\star p$ (discrete counterpart of $-\nabla p$), points \emph{inward} the domain, meaning that people exit \emph{slower} than they would if they were alone (CD phenomenon). The previous example explains the mechanism through which the desired velocity can be balanced by interaction forces. In can be checked that such jams spontaneously appear in numerical simulations, as illustrated by Figure~\ref{fig:StaticJam}, which represents a jammed situation with the desired velocities in black, and the forces induced by pressure in blue.
Figure~\ref{fig:PressureNetwork} represents the pressure field during an evacuation process. Pressures are associated to edges between two centers, and the width of edges are proportional to pressure values. Pressure arches are quite visible, in particular at the very exit, which instantiates the general principle described in Figure~\ref{fig:arche}. Note also that the highest pressures implicate discs upfront, to be compared to the macroscopic situation, where maximal pressures are attained on a zone significantly upstream the exit, so that people at the exit are pushed ahead.
\begin{figure}[tbp]
\includegraphics[width=0.95\linewidth]{StaticJam.pdf}
\caption{Stable static jam}\label{fig:StaticJam}
\end{figure}
\begin{figure}[tbp]
\includegraphics[width=0.95\linewidth]{HD_IsoPresMac.pdf}
\includegraphics[width=0.95\linewidth]{PressureNetwork2.pdf}
\caption{Macroscopic (top) and microscopic (bottom) pressure fields.}\label{fig:PressureNetwork}
\end{figure}
\subsubsection*{Faster is Slower effect}
The Faster is Slower (FiS) effect pertains to a wide class of systems which may perform globally worse as their individual components strive to do better (see e.g.~\cite{FiS15}). In the context of evacuation, it corresponds to a situation where the eagerness of some individuals to exit a room leads to a decrease of the evacuation efficiency. We explain here why the macroscopic model~\eqref{eq:transport}~\eqref{eq:uPU} does \emph{not} make it possible to recover the FiS effect, whereas its microscopic counterpart~\eqref{eq:evolmic} reproduces this effect. We propose a partial explanation of this paradoxical phenomenon, which relies on the \emph{defects} of the underlying Laplacian.
We consider again the situation represented by Figure~\ref{fig:room}, with a desired velocity field $U$ which points toward the exit. In~\cite{CinEqs18} it was proposed to change $U$ by introducing a scalar speed correction field $\beta(x)$, $\beta >1$ corresponding to a speed increase, and to investigate the dependence of the flux through the exit upon $\beta$. We extend here this approach by allowing variations of $U$ in any direction, and we shall keep $\beta$ to denote the variation around $U$, so that $U+\beta $ is the modified velocity field, and we shall investigate the dependence of the flux upon $\beta$. The FiS issue formalizes in this context:
\begin{quote}
Are there any zones such that a variation of $U$ in its own direction (so that $\beta \cdot U >0$) induces a \emph{decrease} of the flux at the exit ?
\end{quote}
The pressure $p_\beta$ associated to this desired velocity $U + \beta$ is the solution to the Poisson problem
\begin{equation}\label{eq:pdef}
- \Delta p_\beta = - \nabla \cdot U - \nabla \cdot \beta,
\end{equation}
with homogeneous Dirichlet boundary conditions $p_\beta=0$ on $\Gamma_{out}$ and $\Gamma_{up}$, and homogeneous Neuman boundary conditions $\partial p_\beta/\partial n =0$ on $\Gamma_{w}$. The associated flux through the exit is thus
\begin{equation}\label{eq:Jdef}
J(\beta) = \int_{\Gamma_{out}} u _\beta \cdot n =
\int_{\Gamma_{out}} U \cdot n -
\int_{\Gamma_{out}} \frac {\partial p_\beta}{\partial n}.
\end{equation}
\begin{prop}\label{prop30}
Let us assume that the desired velocity $U$ is tangent to the wall. Let $\beta \mapsto J(\beta)$ the flux functional defined by~\eqref{eq:Jdef}. We assume that desired velocities are not modified at the exit, i.e. $\beta \equiv 0$ on $\Gamma_{out}$. The gradient of $J$ at $\beta$, that is the field $g$ such that
\[
J\left(\beta + \tilde \beta\right) = J(\beta) + \int_\OO g \cdot \tilde \beta + o\left(\tilde \beta\right),
\]
for all variations $\tilde \beta$ vanishing on $\Gamma_{out}$ is $-\nabla q$, where $q$ is the solution to the adjoint problem
\begin{equation}\label{eq:adjoint}
\left |
\begin{array}{rcll}
-\Delta q &=& 0 & \text{ in }\, \Omega,
\\
q &= & -1 & \text{ on }\, \Gamma_{out} \\
q &= & 0 & \text{ on }\, \Gamma_{up} \vseq\\
\dsp \frac {\partial q }{\partial n } &=& 0 & \text{ on }\, \Gamma_{w}.
\end{array}
\right.
\end{equation}
\end{prop}
\begin{proof}
The gradient is determined by the adjoint problem method. We introduce a dual variable $q$, that is a function defined in $\OO$, and which makes it possible to express the relation between the control variable $\beta$ and the state variable $p$. For any such $q$ (assumed to be smooth), from~\eqref{eq:pdef}, we have that
\begin{equation}\label{eq:wstat}
b(p,q) = \int_\OO \nabla p\cdot \nabla q - \int_\Gamma \frac{\partial p}{\partial n } q +
\int_\OO q \nabla \cdot U + \int_\OO q \nabla \cdot \beta = 0,
\end{equation}
where the boundary integral over $\Gamma$ can be replaced by an integral over $\Gamma_0 = \Gamma_{out} \cup \Gamma_{up}$, thanks to the boundary condition on $\Gamma_{w}$. Now, since we assume that desired velocities are not modified at the exit, and since we are interested in variations of the flux only, we may drop the first term of~\eqref{eq:Jdef} in the definition of $J$ :
\[
J(\beta) = -
\int_{\Gamma_{out}} \frac {\partial p_\beta}{\partial n}.
\]
The Lagrangian
\[
L\left(p,\beta,q\right) = \underbrace{ -
\int_{\Gamma_{out}} \frac {\partial p}{\partial n} }_{\flux} +\underbrace{ \int_\OO \nabla p\cdot \nabla q - \int_{\Gamma_0} \frac{\partial p}{\partial n } q +
\int_\OO q \nabla \cdot U + \int_\OO q \nabla \cdot \beta }_{b\left(p,q\right) }
\]
is such that, for any $p_\beta $ coupled with $\beta$ by the state equation, $b(p_\beta,q) = 0$ for any $q$. Thus, for any $\beta$ and any $q$,
\[
L\left(p_\beta, \beta, q\right) = J(\beta).
\]
We differentiate with respect to $\beta$ both sides of this identity:
\begin{equation}\label{eq:Dbeta}
D_\beta J= D_pL \circ D_\beta p_\beta + D_\beta L.
\end{equation}
The approach consists in choosing an adjoint variable $q$ such that $D_p L = 0$, which circumvents the difficulty to explicitly estimate $D_\beta p_\beta$. Let $\tilde p $ denote a variation in the variable $p$. It holds that
\begin{align*}
D_pL \, \tilde p & = -
\int_{\Gamma_{out}} \frac {\partial \tilde p}{\partial n} + \int_\OO \nabla \tilde p\cdot \nabla q - \int_{\Gamma_{out}} \frac{\partial \tilde p}{\partial n } q - \int_{\Gamma_{up}} \frac{\partial \tilde p}{\partial n } q\\
&= - \int_{\Gamma_{out}} \frac {\partial \tilde p}{\partial n} \left(q+1\right) + \int_\OO \left(-\Delta q\right) \tilde p +
\int_{\Gamma} \frac{\partial q}{\partial n } \tilde p - \int_{\Gamma_{up}} \frac{\partial \tilde p}{\partial n } q.
\end{align*}
Since $\tilde p$ is a variation of $p$, which identically vanishes on $\Gamma_0$, the boundary integral over $\Gamma = \partial \Omega$ reduces to an integral over $\Gamma_w$. The adjoint problem is designed in order to make the previous expression vanish:
\[
- \Delta q = 0\text{ in } \Omega
\virg q = - 1 \text{ on } \Gamma_{out}\virg \frac{\partial q}{\partial n } = 0 \text{ on } \Gamma_{w}\virg q = 0 \text{ on } \Gamma_{up},
\]
which is exactly~\eqref{eq:adjoint}, up to a change in sign in the Dirichlet boundary condition. Let $q$ be the solution to this adjoint problem, so that $D_pL(p_\beta,\beta,q) = 0$. From~\eqref{eq:Dbeta}, it comes
\begin{align*}
D_\beta J &= 0 + D_\beta L \left(p_\beta,\beta,q\right),
\\
\intertext{with}
D_\beta J\, \tilde \beta &= \int_\OO q \nabla \cdot \tilde \beta = - \int _\Omega \tilde \beta \cdot \nabla q + \int_\Gamma q\tilde \beta \cdot n.
\end{align*}
Since $\tilde \beta = 0$ on $\Gamma_{out}$, $q = 0 $ on $\Gamma_{up}$, and $U\cdot n = 0$ on $\Gamma_{w}$, the boundary term vanishes, which yields
\[
\nabla _\beta J = - \nabla q,
\]
which ends the proof of Proposition~\ref{prop30}.
\end{proof}
Figure~\ref{fig:HD_grad} represents $- \nabla q$. The fact that it is the gradient of $J(\beta)$, that is the outflow through the exit, means that, at any point $x$ of the saturated domain, if a person located around $x$ change their desired velocity by a vector $\beta$, the effect upon the flux will be proportional to $- \nabla q \cdot \beta$. As shown in the figure, the gradient of $J(\beta)$ is strikingly similar to a desired velocity field, which means that any local variation of the desired velocity in its own direction will lead to an increase of the flux through the door, which corresponds to a \emph{Faster is Faster} effect.
Let us now show that the situation is fully different in the microscopic setting. As previously, we shall investigate the dependence of some observable quantifying the evacuation upon variations of individual desired velocities of people upstream the exit. We aim here at determining whether it might be possible that some agents would improve the evacuation efficiency by reducing their pace. In the discrete situation, the instantaneous flux through the exit is not a relevant marker of egress efficiency. We shall rather consider the velocity of a particular agent close to the exit. We consider a highly congested situation like the one represented in Figure~\ref{fig:ni}. We propose to formulate the FiS issue in the following terms: consider a person (i.e. a disk) of index $i$ on the front of the cluster, i.e. close to the exit (see Figure~\ref{fig:ni}).
\begin{figure}[tbp]
\includegraphics[width=0.65\linewidth]{ni.pdf}
\caption{Definition of $i$ and $n_i$.}\label{fig:ni}
\end{figure}
We denote by $n_i$ the direction of their desired velocity. Now consider that any other individual $j$ has the ability to slightly change their speed, i.e. by changing their desired velocity to $ U_j + \beta_j$. The FiS issue amounts to investigate whether taking $\beta _j$ in the direction of $U_j$ for some individuals might \emph{decrease} the velocity of $i$ in the desired direction, thereby decreasing the flow near the exit. To formalize this question, we first consider the actual velocity field associated to the desired velocity field $U$ at some instant, i.e. the couple $(u,p)$ solves System~\eqref{eq:uBUdisc}. As in Proposition~\ref{prop:oB}, we eliminate all those rows of $B$ that do not correspond to active constraints. We keep the same notation for the reduced $B$ and the associated reduced vector $p$ of Lagrange multipliers. The system writes
\begin{equation}\label{eq:PSBred}
\left |
\begin{array}{lcl}
\dsp u + B^\star p &=& U,\vseq\\
\qquad B u & = & 0,
\end{array}
\right.
\end{equation}
where all pressures are positive. The velocity can be eliminated, which yields
\[
BB ^\star p = BU,
\]
which is the discrete counterpart of the macroscopic Poisson problem ($-\Delta p = -\nabla \cdot U$).
We now consider the velocity as a function $u_\beta$ of $\beta \in \setR^{2N}$, i.e. the desired velocity of individual $j$ is $ U_j + \beta_j$. Admitting that the network of active contacts remains the same for small variations of $\beta$, it holds that
\[
BB ^\star p_\beta = B U + B \beta \virg u_\beta= U + \beta - B^\star p_\beta.
\]
The objective function we are interested in is the actual velocity of $i$ in its desired direction:
\begin{equation}\label{eq:obj}
u_\beta \cdot n_i = U_i \cdot n_i + \beta \cdot n_i - B^\star p_\beta \cdot n_i
\virg n_i = \left(0,\,\dots,\,0,U_i / \abs { U_i}, 0,\,\dots,\,0\right).
\end{equation}
Since person $i$ is the closest to the exit, with no one in front, any change of $i$'s desired velocity shall affect their actual velocity in the same direction. We shall therefore disregard variations of speed for $i$, so that the first two terms above are constant, and we define the objective function as
\begin{equation}\label{eq:micJdef}
J(\beta) = - B^\star p_\beta \cdot n_i.
\end{equation}
\begin{figure}[tbp]
\includegraphics[width=0.75\linewidth]{HD_gradRed.pdf}
\includegraphics[width=0.70\linewidth]{toto.pdf}
\caption{Solution to the adjoint problem, in the macroscopic (top) and microscopic (bottom) settings.}\label{fig:HD_grad}
\end{figure}
\begin{prop}\label{rop:adjointmicro}
Let $\beta \mapsto J(\beta)$ be the flux functional defined by~\eqref{eq:micJdef}. The gradient of $J$ with respect to $\beta$ is $-B^\star q$, where $q$ is the solution to the adjoint problem
\begin{equation}\label{eq:micadjoint}
BB ^\star q = Bn_i.
\end{equation}
\end{prop}
\begin{proof}
Like in the macroscopic setting, the core of the approach relies on a so-called \emph{Lagrangian}, that is a function of the state variable $p$, the control variable $\beta$, and a dual variable $q$. It is defined as the sum of the objective function (velocity of $i$ along the desired direction) expressed in the state variable $p$ (uncoupled from the control variable $\beta$), and a weak expression of the state equation $BB ^\star p_\beta = BU + B \beta$:
\[
L(p,\beta,q) = - B^\star p\cdot n_i +\left(BB^\star p - BU - B \beta\right) \cdot q.
\]
As a consequence, for any $\beta$ and any $q$,
\[
L\left(p_\beta, \beta, q\right) = J(\beta).
\]
Differentiation with respect to $\beta$ both sides of this identity yields Eq.~\eqref{eq:Dbeta}. The approach consists again in choosing an adjoint variable $q$ such that $D_p L = 0$. Let $\tilde p $ denote a variation in the variable $p$. It holds that
\[
D_pL \, \tilde p = - B^\star \tilde p\cdot n_i +\left(BB ^\star \tilde p\right) \cdot q =
\left(-Bn_i +BB ^\star q\right) \cdot \tilde p.
\]
The adjoint problem is designed in order to make the previous expression vanish, it therefore reads
\begin{equation}\label{eq:BBsq}
BB ^\star q = Bn_i.
\end{equation}
On the other hand, it holds that
\[
D_\beta L \, \tilde \beta = - B \tilde \beta \cdot q = - B \star q \cdot \tilde \beta.
\]
We finally obtain
\[
\nabla J(\beta) = - B^\star q
\]
where $q$ is the solution to the adjoint problem~\eqref{eq:BBsq}.
\end{proof}
This expression of the gradient makes it possible to investigate the Faster is Slower effect. Figure~\ref{fig:HD_grad} represents a congested situation, the black disk corresponds to $i$, the agent upfront. The black arrows represent to $\nabla J $, i.e. $- B^\star q$, where $q$ is the solution to the adjoint problem. This field is the straight discrete counterpart to the continuous field represented on the top of the figure, yet it presents very different characteristics. Contrary to the macroscopic field, some arrows indeed point \emph{upstream} the exit. For all such agents (represented in blue in the figure), changing their velocity against their desired direction~(i.e. slowing down) will \emph{increase} the speed of $i$, i.e. it will speed up the evacuation process.
\section{Abstract considerations}\label{sec:abs}
We propose here an abstract setting to handle the macroscopic and microscopic saddle-point formulations (systems~\eqref{eq:uBU} and~\eqref{eq:uBUdisc}, respectively), and we give an abstract criterium to distinguish those situations which lead to a maximum principle from those which do not.
\subsubsection*{Abstract saddle-point formulations}
The abstract setting is represented in Figure~\ref{fig:VV} (middle). The primal space $V$ is $L^2(\O)$ in the macroscopic setting, and $\setR^{2N}$ in the microscopic one. The figure represents the different cones which can be introduced. The left-hand side of the figure represents $C$, the feasible cone, and its polar
\[
C^\circ = \set{ v \in V \virg \left\langle v\,|\,u\right\rangle \leq 0 \quad \forall u \in C},
\]
in $(\ker B)^\bot$, so that $\ker B$ is reduced to single point, which is $C\cap C^\circ$, in this representation.
\begin{prop}\label{prop:absspp}
Let $V$ and $\Lambda$ be Hilbert spaces, $B\in \clL(V,\Lambda)$, $\Lambda_+$ a closed convex cone in $\Lambda$. We define its polar cone as
\[
\Lambda _+^ \circ = \set{q\in \Lambda \virg \left\langle q\,\middle|\,p\right\rangle\leq 0 \quad \forall p \in \Lambda_+}.
\]
For a given $U\in V$, we consider the problem which consists in minimizing
\[
v \longmapsto J(v) = \frac 1 2 \abs { v- U}^2,
\]
over $C \subset V$ defined by
\[
C = \set{ v\in V\virg \left\langle Bv\,|\,q\right\rangle \leq 0\quad \forall q \in \Lambda_+} = B^{-1}\left(\Lambda_+^\circ\right).
\]
There exists a unique minimizer for $J$ over $C$ and, in any of the following situations:
\begin{enumerate}\romanenumi
\item $\Lambda_+$ is spanned by a finite number of vectors, i.e. $\Lambda_+ = \corrl(g_1,\dots,g_p)$ (conic convex hull).
\item $B$ is surjective.
\end{enumerate}
the problem admits a variational formulation in the following sense: $u$ is the constrained minimizer if and only if there exists $p$ such that
\begin{equation}\label{eq:uBUabs}
\left |
\begin{array}{rcl}
u + B^\star p& = &U \\
B u &\in & \Lambda_+^\circ \vseq \\
p &\in& \Lambda _+ \vseq \\
\left\langle p\,\middle|\,Bu\right\rangle&=& 0.
\end{array}
\right.
\end{equation}
\end{prop}
\begin{proof}
The Hilbert space $V$ is decomposed into the sum of $C$ and its polar cone $C^\circ$, i.e. any $U\in V$ decomposes uniquely as the sum
\[
U = u + u^\circ\virg u \in C\virg u^o \in C^o\virg \left\langle u\,\middle|\,u^o\right\rangle = 0,
\]
which can be seen as a unilateral version of the direct sum of two orthogonal subspaces (see~\cite{moreau62}). Let us prove now that, in both cases~\ref{case1} and~\ref{case2}, $u^o$ can be written $B^\star p$ with $p\in\Lambda_+$.
%Case $(i)$:
\begin{case}\label{case1}
If $\Lambda_+$ is the closed convex conic hull of a finite number of vectors, i.e. if
\[
\Lambda_+ = \corrl\left(g_1,\,\dots,\,g_p\right) = \Big\{ \sum \lambda_i g_i\virg \text{ with } \lambda_i \geq 0\quad \forall i = 1,\,\dots,\,p\Big\},
\]
then $C$ can be written
\[
C = \set{ v \in V\virg \left\langle v\,\middle|\,B^\star g_i\right\rangle \leq 0\quad \forall i = 1,\,\dots,\,p},
\]
and the identity
\[
C^\circ = \corrl\left(B^\star g_1,\,\dots,\,B^\star g_p\right) = \Big\{ \sum \lambda_i B^\star g_i\virg \text{ with } \lambda_i \geq 0\quad \forall i = 1,\,\dots,\,p\Big\}
\]
is a straight consequence of Farkas' Lemma. Note that this lemma is usually expressed in a finite dimensional setting, but it can be easily checked that it generalizes to infinite dimensional spaces as soon as the number of unilateral constraints is finite (the cone has a finite number of facets). The previous identity implies that $u^\circ \in C^\circ$ can be written $u^\circ = B^ \star p$, with $p\in \lambda_+$, which settles Case~\ref{case1}.
\end{case}
%Case $(ii)$:
\begin{case}\label{case2}
In all generality it holds that
\[
C^\circ = \overline{B^\star \Lambda_+}
\]
Indeed, $C^\circ \subset \overline{B^\star \Lambda_+}$ is trivial, and if the inclusion is strict, there exists $v\in \overline{B^\star \Lambda_+}$, with $v\notin C^\circ$. By Hahn-Banach Theorem, there exist $w\in V$ and $\alpha\in \setR$ such that
\[
\left\langle B^\star q, w\right\rangle \leq \alpha < \left\langle v\,|\,w\right\rangle\quad \forall q\in \Lambda_+.
\]
Since $\Lambda _+$ is a cone, the first inequality implies that $\langle B^\star q\,|\, w\rangle\leq 0$ for any $q\in \Lambda_+$, so that $w\in C$, and $\alpha \geq 0$. We then have $\langle v\,|\,w\rangle0$, with $v\in C^\circ$ and $w\in C$, which is impossible.
\begin{figure}[tbp]
\includegraphics[width=0.90\linewidth]{VV2.pdf}
\caption{Abstract setting.}\label{fig:VV}
\end{figure}
As a consequence, if $B^\star \Lambda_+$ is closed, then $C^\circ =B^\star \Lambda_+$. Let us prove that $B^\star \Lambda_+$ is indeed closed. If $B$ is surjective between the Hilbert spaces $V$ and $\Lambda$, then there exists $\gamma >0$ such that (continuous inf-sup condition)
\[
\inf_\Lambda \sup_V \frac {\left\langle Bv\,\middle|\,q\right\rangle}{\abs v \abs q} \geq \gamma
\Longrightarrow
\abs{B^\star q} \geq \gamma \abs {q}\quad \forall q \in \Lambda.
\]
As a consequence, if a sequence $B^\star q_n$ converges in $V$, with $q_n\in \Lambda_+$, then $q_n$ is a Cauchy sequence that converges to $q\in \Lambda_+$, and the previous limit is $B^\star q\in B^\star \Lambda_+$, which proves that $B^\star \Lambda_+$ is closed, and thereby ends the proof of Proposition~\ref{prop:absspp}.\qedhere
\end{case}\let\qed\relax
\end{proof}
\begin{rema}
In the case of equality affine constraints, assuming that $B(V)$ is closed is sufficient to ensure existence of a Lagrange multiplier, whereas surjectivity ensures its uniqueness. In the present situation, with inequality constraints, assuming surjectivity, as we did, is essential to ensure existence. Indeed, an operator with closed range maps a closed convex cone to a convex cone which is not necessarily closed: if $B^\star (\Lambda)$ is closed (which is equivalent to say that $B$ has a closed range), it does not imply that $B^\star (\Lambda_+)$ is closed, even in the finite dimensional setting, we refer to~\cite{BourdinMaury2020} for a counter-example.
\end{rema}
\subsubsection*{Abstract maximal principle}
The difference between the macroscopic and the microscopic models can be described within the abstract setting presented above. In the macroscopic case, if one considers a desired velocity that is concentrating, i.e. $-\nabla \cdot U \geq 0$, it writes $BU\in - C$ (the operator $B $ is $-\nabla \cdot $). The fact that $-\Delta p = -\nabla\cdot U\geq 0 $ (which writes $BB^\star p = B U$), admits a solution that is non-negative means that there exists $p\in \Lambda_+$ such that $U - B^\star p$ is in $\ker B$. In other words, the maximal principle takes the abstract form
\[
- C \subset C^ \circ + \ker B,
\]
which is equivalent to $ - C \subset B^\star (\Lambda_+) + \ker B$. Figure~\ref{fig:VV} (middle), which gives a representation in $(\ker B)^\bot$, correspond to this situation : $-C$ appears as included in $C^ \circ$, and the vector $U'$ (in green) corresponds to a concentrating field (i.e. with nonpositive divergence).
The microscopic situation corresponds to the situation represented in the bottom. The inclusion above no longer holds :
\[
- C \not \subset B^\star (\Lambda_+) + \ker B.
\]
There exists vectors $U'$ in $-C$, i.e. concentrating fields like the one represented in Figure~\ref{fig:coin} (top-right), that are \emph{not} in $C^\circ + \ker B$ so that there exists no pressure fields in $\Lambda _+ = \setR_+^{N_c}$ such that $BB^\star p = B U$, which expresses the pathology of this Laplacian $BB^\star$.
\section{Conclusion}
We have recalled in this paper the tight links between continuous and discrete Laplacians, in contexts where the discrete Laplacian combines a phenomenological law of the Ohm's type together with a conservation principle. In the context of crowd motions (or granular flows), accounting for the non-overlapping constraint between solid grains can be done by duality, which leads to a unilateral discrete Darcy problem. Like in the continuous setting, eliminating the velocity leads to a Poisson-like problem for a pressure field defined on the dual graph of the contact network. The discrete Laplacian $BB^\star$ is the straight counterpart of $-\nabla \cdot \nabla = - \Delta$ in the one-dimensional situation (people / grains in a row), but the situation is fully different in higher dimensions. In particular, the matrix $BB^\star$ does not verify any maximum principle in general, and it does vanish on constant fields. We have investigated how these very defects make it possible to recover effects that are observed in reality, whereas in the continuous setting, with a standard Laplacian verifying the maximum principle, these effects are not reproduced.\linebreak These considerations also rule out the possibility to properly obtain the continuous model as a many-body limit of the discrete one, as illustrated for example by Figures~\ref{fig:PressureNetwork} and~\ref{fig:HD_grad}, which enlight the deep discrepancies between both scales, although the corresponding models are built upon the very same principles.
\section*{Conflicts of interest}
The authors declare no competing financial interest.
\section*{Acknowledgments}
The author would like to thank S. Faure for his contributions in numerical computations, and R. Glowinski for stimulating discussions over the last decades, which lead to these developments.
%\nocite{*}
\bibliographystyle{crunsrt}
\bibliography{crmeca20221259}
\end{document}