# Compositional maps for registration in complex geometries

###### Abstract

We develop and analyze a parametric registration procedure for manifolds
associated with the solutions to parametric partial differential equations in
two-dimensional domains. Given the domain $\Omega\subset\mathbb{R}^{2}$ and
the manifold $\mathcal{M}=\\{u_{\mu}:\mu\in\mathcal{P}\\}$ associated with the
parameter domain $\mathcal{P}\subset\mathbb{R}^{P}$ and the parametric field
$\mu\mapsto u_{\mu}\in L^{2}(\Omega)$, our approach takes as input a set of
snapshots from $\mathcal{M}$ and returns a parameter-dependent mapping
$\Phi:\Omega\times\mathcal{P}\to\Omega$, which tracks coherent features (e.g.,
shocks, shear layers) of the solution field and ultimately simplifies the task
of model reduction. We consider mappings of the form
$\Phi=\texttt{N}(\mathbf{a})$ where $\texttt{N}:\mathbb{R}^{M}\to{\rm
Lip}(\Omega;\mathbb{R}^{2})$ is a suitable linear or nonlinear operator; then,
we state the registration problem as an unconstrained optimization statement
for the coefficients $\mathbf{a}$. We identify minimal requirements for the
operator N to ensure the satisfaction of the bijectivity constraint; we
propose a class of compositional maps that satisfy the desired requirements
and enable non-trivial deformations over curved (non-straight) boundaries of
$\Omega$; we develop a thorough analysis of the proposed ansatz for polytopal
domains and we discuss the approximation properties for general curved
domains. We perform numerical experiments for a parametric inviscid transonic
compressible flow past a cascade of turbine blades to illustrate the many
features of the method.

Tommaso Taddei1

1 Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France
Inria Bordeaux Sud-Ouest, Team MEMPHIS, 33400 Talence, France,
<EMAIL_ADDRESS>

_Keywords:_ parameterized partial differential equations; registration in
bounded domains; model order reduction.

##  1 Introduction

###  1.1 Registration in bounded domains

Numerical methods based on Lagrangian (registration-based) approximations have
proven to be a promising research direction in scientific computing, both for
the numerical discretization of partial differential equations (PDEs) [29, 39]
and also for model order reduction of parametric systems (MOR, [11, 20, 21,
28, 32]). Given the field of interest $u$ defined over the domain
$\Omega\subset\mathbb{R}^{2}$, Lagrangian methods seek approximations of the
solution $u$ of the form $\tilde{u}\circ\Phi^{-1}$ where $\tilde{u}$ belongs
to a linear approximation space and $\Phi$ is a bijection $\Omega$. This class
of methods is designed for problems with sharp features such as shocks,
contact discontinuities or shear layers in fluid mechanics and fractures in
solid mechanics. The problem of finding the mapping $\Phi$ is dubbed as
registration problem and shares important features with geometry registration
[18, 23] mesh morphing [31, 36] techniques, and also optimal transportation
[5, 13, 26]. The goal of this paper is to develop a general registration
procedure for parametric problems with emphasis on parametric MOR
applications, and to provide a rigorous mathematical analysis of the problem
of registration in bounded two-dimensional domains.

We pursue an optimization-based approach to the problem of registration. We
denote by $\mu$ a vector of $P$ parameters in the parameter domain
$\mathcal{P}\subset\mathbb{R}^{P}$; given the domain
$\Omega\subset\mathbb{R}^{2}$, we denote by $\mathfrak{B}$ the space of
Lipschitz bijections from $\Omega$ in itself, and by $\mathfrak{D}$ the space
of diffeomorphisms from $\Omega$ in itself. Given the mapping $\Phi$, we
denote by $J(\Phi)$ the Jacobian determinant, and we denote by
$\texttt{id}:\mathbb{R}^{2}\to\mathbb{R}^{2}$ the identity map in
$\mathbb{R}^{2}$, $\texttt{id}(x)=x$ for all $x\in\mathbb{R}^{2}$. If we fix
the value of $\mu\in\mathcal{P}$, our goal is to minimize a target function
$\mathfrak{f}_{\mu}^{\rm tg}$ over all possible diffeomorphisms of $\Omega$,

$\min_{\Phi\in\mathfrak{D}}\mathfrak{f}_{\mu}^{\rm tg}(\Phi).$ (1)

Problem (1) is computationally intractable due to the fact that $\mathfrak{D}$
is an highly non-convex subset of $C^{1}(\Omega;\mathbb{R}^{2})$: research on
registration should hence focus on the development of computational strategies
to devise tractable counterparts of the statement (1).

###  1.2 Optimization-based registration

To devise a tractable registration procedure, we introduce an operator
$\texttt{N}:\mathbb{R}^{M}\to{\rm Lip}(\Omega;\mathbb{R}^{2})$ and a penalty
function $\mathfrak{f}_{\rm pen}:\mathbb{R}^{M}\to\mathbb{R}_{+}$ such that

$\left\\{\begin{array}[]{l}\displaystyle{\mathcal{B}_{\texttt{N}}=\left\\{\texttt{N}(\mathbf{a})\;:\;\mathfrak{f}_{\rm
pen}(\mathbf{a})\leq C\right\\}\subset\mathcal{B},\quad{\rm
for\;some}\;C>0;}\\\\[8.53581pt]
\displaystyle{\texttt{N}(\mathbf{a}=0)=\texttt{id},\;\;\mathfrak{f}_{\rm
pen}(\mathbf{a}=0)<C;}\\\\[8.53581pt]
\displaystyle{\texttt{N},\;\mathfrak{f}_{\rm pen}\;{\rm
are\;Lipschitz\;continuous}.}\\\ \end{array}\right.$ (2)

Equations (2)1 and (2)2 imply that the set $\mathcal{B}_{\texttt{N}}$ is not
empty; furthermore, exploiting (2)3 we find that the interior of
$\mathcal{B}_{\texttt{N}}$ is not empty — that is, bijectivity of the mapping
$\texttt{N}(\mathbf{a})$ is preserved for small perturbations of the mapping
coefficients $\mathbf{a}$. We further observe that, given the full rank matrix
$\mathbf{W}\in\mathbb{R}^{M\times m}$, the pair
$(\widetilde{\texttt{N}},\widetilde{\mathfrak{f}_{\rm pen}})$ such that
$\widetilde{\texttt{N}}(\cdot)=\texttt{N}(\mathbf{W}\cdot)$ and
$\widetilde{\mathfrak{f}_{\rm pen}}(\cdot)=\mathfrak{f}_{\rm
pen}(\mathbf{W}\cdot)$ satisfies (2): we can hence apply linear compression
methods such as proper orthogonal decomposition (POD, [30, 37]) to the mapping
coefficients $\mathbf{a}$ without fundamentally changing the properties of our
ansatz; as discussed in sections 2 and 4, this feature greatly simplifies the
task of dimensionality reduction for parametric problems.

Exploiting the previous definitions, we introduce the surrogate of (1):

$\min_{\mathbf{a}\in\mathbb{R}^{M}}\mathfrak{f}_{\mu}^{\rm
obj}(\mathbf{a}):=\mathfrak{f}_{\mu}^{\rm
tg}(\texttt{N}(\mathbf{a}))\;+\xi\;\mathfrak{f}_{\rm pen}(\mathbf{a}).$ (3)

Provided that $\xi$ is sufficiently large, solutions $\mathbf{a}^{\star}$ to
(3) satisfy the condition $\mathfrak{f}^{\rm pen}(\mathbf{a}^{\star})$ $\leq
C$. Note that (3) reads as a nonlinear non-convex unconstrained optimization
problem that can be tackled using standard gradient-descent optimization
algorithms.

The previous discussion highlights the two major questions in registration
methods: (i) how to construct $(\texttt{N},\mathfrak{f}_{\rm pen})$ that
satisfy (2) for a given domain $\Omega$; (ii) how to establish a rigorous
relation between the solutions to (1) and to (3) in the limit $M\to\infty$.
Note that the second question concerns the ability of approximating arbitrary
elements of $\mathfrak{D}$ using operators N that satisfy (2). We notice that
even if we focus on the approximation of diffeomorphisms we allow ourselves to
consider approximations in a less regular space: this choice is justified by
the particular discretization method (the finite element (FE) method) employed
in this work to represent the operator N and, more fundamentally, by the
strategy proposed here to define N.

###  1.3 Compositional maps for registration

The objective of this paper is to devise pairs $(\texttt{N},\mathfrak{f}_{\rm
pen})$ that satisfy (2) and study the approximation properties in
$\mathfrak{D}$ for arbitrary curved Lipschitz domains $\Omega$. In this work,
we say that the boundary of $\Omega$ is curved if it cannot be represented as
the union of a finite number of straight lines; if $\partial\Omega$ is curved,
we refer to $\Omega$ as to “curved domain”.

Towards this end, we first show that for polytopal domains affine maps
$\texttt{N}_{\rm
p}(\mathbf{a})=\texttt{id}+\sum_{i=1}^{M}(\mathbf{a})_{i}\varphi_{i}$ for
suitably chosen functions
$\\{\varphi_{i}\colon\mathbb{R}^{2}\to\mathbb{R}^{2}\\}_{i}$ and a suitable
penalty function $\mathfrak{f}_{\rm pen}$ satisfy (2) (cf. Proposition 2.1 and
Corollary 2.3); the function $\\{\varphi_{i}\\}_{i}$ should be chosen so that
their normal components vanish on $\partial\Omega$,
$\varphi_{i}\cdot\mathbf{n}|_{\partial\Omega}=0$ for $i=1,\ldots,M$, where
$\mathbf{n}:\partial\Omega\to\mathbb{S}^{1}=\\{x\in\mathbb{R}^{2}:\|x\|_{2}=1\\}$
is the outward normal to $\Omega$. We also prove that affine maps of this form
are dense for $M\to\infty$ in a meaningful subspace of diffeomorphisms (cf.
Proposition 2.2 and Corollary 2.4).

As rigorously shown in Lemma 3.1, affine maps are fundamentally ill-suited to
approximate non-trivial diffeomorphisms in curved domains; to address this
issue, we propose and analyze compositional maps of the form

$\texttt{N}(\mathbf{a})=\Psi\circ\texttt{N}_{\rm
p}(\mathbf{a})\circ\Psi^{-1},$ (4)

where $\Psi:\Omega_{\rm p}\to\Omega$ is a bijection from the polytope
$\Omega_{\rm p}$ to $\Omega$ and $\texttt{N}_{\rm p}(\mathbf{a}):\Omega_{\rm
p}\to\Omega_{\rm p}$ satisfies $\texttt{N}_{\rm
p}(\mathbf{a})=\texttt{id}+\sum_{i=1}^{M}(\mathbf{a})_{i}\varphi_{i}$ with
$\varphi_{i}\cdot\mathbf{n}|_{\partial\Omega_{\rm p}}=0$ for $i=1,\ldots,M$.
We propose an actionable strategy to define the polytope $\Omega_{\rm p}$ and
the mapping $\Psi$: our approach relies on the definition of a coarse-grained
curved high-order FE mesh $\mathcal{T}$ of $\Omega$, and to a fully-automated
procedure to define the polytope $\Omega_{\rm p}$ and the map $\Psi$ based on
the curved mesh $\mathcal{T}$.

The bijection $\Psi$ is independent of the mapping coefficients $\mathbf{a}$
and is designed to recast the registration problem of interest into a polytope
$\Omega_{\rm p}$ where affine maps can be employed. Figure 1 provides a
graphic representation of the sought mapping $\Psi$. Consider the problem of
finding a bijection of the semicircular domain $\Omega$ such that $\Phi(a)=b$:
first, we introduce the bijection $\Psi$ from the polytope (triangle)
$\Omega_{\rm p}$ to the curved domain $\Omega$; then, we seek a bijection
$\Phi_{\rm p}$ of $\Omega_{\rm p}$ such that $\Phi_{\rm p}(a_{\rm p})=b_{\rm
p}$ with $a_{\rm p}=\Psi^{-1}(a)$ and $b_{\rm p}=\Psi^{-1}(b)$. By
construction, the resulting mapping $\Phi=\Psi\circ\Phi_{\rm p}\circ\Psi^{-1}$
is a Lipschitz bijection of $\Omega$ that maps $a$ into $b$.

$a$$b$$\Omega$$\Omega_{\rm p}$${\Psi}$$a_{\rm p}$$b_{\rm p}$ Figure 1:
compositional maps. The bijection $\Psi$ is designed to recast the
registration task from the curved domain $\Omega$ to the polytope $\Omega_{\rm
p}$.

We exploit the analysis of registration in polytopes to study the properties
of maps of the form (4) and we discuss the approximation power. We show (cf.
Lemma 3.2) that single-layer compositional maps cannot approximate arbitrary
diffeomorphisms in curved domains. To improve the approximation power, we
consider the multi-layer generalization of (4) such that

$\texttt{N}\left(\mathbf{a}=[\mathbf{a}_{1},\ldots,\mathbf{a}_{\ell}]\right)=\texttt{N}_{1}\left(\mathbf{a}_{1}\right)\circ\ldots\circ\texttt{N}_{\ell}\left(\mathbf{a}_{\ell}\right);$
(5)

here, $\texttt{N}_{i}\left(\mathbf{a}_{i}\right)=\Psi_{i}\circ\texttt{N}_{{\rm
p},i}(\mathbf{a}_{i})\circ\Psi_{i}^{-1}$ where $\Psi_{i}:\Omega_{{\rm
p},i}\to\Omega$ is a Lipschitz bijection and $\Omega_{{\rm p},i}$ is a
suitable polytope, for $i=1,\ldots,\ell$. In Lemma 3.3, we show the superior
properties of the ansatz (5).

Our method is related to several previous works. A first extension of the
registration procedure in [32] to arbitrary domains was proposed in [34]: the
work of [34] relies on the introduction of a coarse-grained partition of the
domain $\Omega$ and on Gordon-Hall maps to morph each element of the partition
into the unit square. The approach in [34] requires that each element of the
partition is mapped in itself (_local bijectivity_) and relies on Gordon-Hall
maps: it is hence very sensitive to the choice of the coarse-grained
partition, it is not dense in any meaningful subset of diffeomorphisms and it
cannot be extended to three-dimensional domains. We also recall the work by
Zahr and Persson [40] for high-order implicit shock tracking methods: the
approach in [40] relies on a local parameterization of boundary degrees of
freedom; for this reason, it cannot be readily combined with linear
dimensionality reduction techniques to identify low-rank mapping spaces for
parametric systems. We finally remark that the multi-layer ansatz (5) is
closely related to registration methods appeared in the image processing
literature [7, 8].

We here rely on a standard $H^{1}$-conforming FE discretization to represent
the mapping; on the other hand, as in our previous works (e.g., [32, 34]), we
aim to include an $H^{2}$ penalization term in the objective
$\mathfrak{f}^{\rm obj}$. Towards this end, we propose a discrete $H^{2}$
broken norm that is inspired by the work by Mozolevski and coauthors [22] on
discontinuous Galerkin (DG) discretizations of the biharmonic equation.
Following the seminal work by Argyris and coauthors [1], several researchers
have considered $H^{2}$-conforming FE spaces for the discretization of fourth-
order operators: the use of standard FE discretizations simplifies the
implementation and also enables the application of polynomial bases of
arbitrary order.

The outline of the paper is as follows. Section 2 provides a complete analysis
of registration in polytopes; section 3 addresses the extension to curved
domains; section 4 discusses the construction of the polytope $\Omega_{\rm p}$
and the mapping $\Psi$ and reviews the parametric registration method proposed
in [34] and employed in the numerical experiments; section 5 contains
numerical investigations to illustrate the performance of the registration
method and its implications for model reduction. We here couple our
registration procedure with a non-intrusive (POD+regression) MOR procedure for
state estimation; the integration of registration in the offline-online
computational paradigm of projection-based MOR is the subject of ongoing
research (see [3]). Section 6 concludes the paper.

##  2 Affine maps in polytopes

In order to find N and $\mathfrak{f}^{\rm pen}$ that satisfy (2), we introduce
the sets

$A:=\\{\mathbf{a}\in\mathbb{R}^{M}\,:\,\texttt{N}(\mathbf{a})\in\mathfrak{B}\\},\quad
A_{\rm
jac}:=\\{\mathbf{a}\in\mathbb{R}^{M}:\inf_{x\in\Omega}J(\texttt{N}(\mathbf{a}))>0\\}.$
(6)

Then, we require that the condition $\inf_{x\in\Omega}J(\Phi)>0$ implies
bijectivity in $\Omega$ for any mapping $\Phi$ spanned by N, that is $A_{\rm
jac}\subset A$. It is possible to construct maps $\Phi$ that are bijective in
$\Omega$ for which $\inf_{x\in\Omega}J(\Phi)=0$ (cf. [27, Theorem 1.1]);
however, bijections $\Phi$ which satisfy $\inf_{x\in\Omega}J(\Phi)=0$ are of
little practical interest for scientific computing applications and in
particular for (projection-based) MOR.

Note that the pointwise condition $J(\texttt{N}(\mathbf{a}))(x)>0$ for all
$x\in\Omega$ cannot be directly translated into a penalty term for (3): we
address the construction of the penalty term in section 2.1. In section 2.2,
we present two technical results for affine maps in polytopes; in sections 2.3
and 2.4, we construct operators N that satisfy $A_{\rm jac}\subset A$ and are
dense in a meaningful subspace of diffeomorphisms; finally, in section 2.5 we
comment on the approximation of Lipschitz maps.

###  2.1 Definition of the penalty function

We denote by $\|\cdot\|_{\rm F}$ Frobenius norm, and by $H(w)$ the Hessian of
the field $w$. Given the polytope $\Omega_{\rm p}$, the constants $\epsilon>0$
and $C_{\rm exp}$ such that $C_{\rm exp}\ll\epsilon$, we introduce the
function $\mathfrak{f}_{\rm pen}^{(1)}(\Phi)=\mathfrak{f}_{\rm
jac}(\Phi)\>+\;\|\nabla J(\Phi)\|_{L^{\infty}(\Omega_{\rm p})}$ with

$\mathfrak{f}_{\rm jac}(\Phi)=\frac{1}{|\Omega_{\rm p}|}\int_{\Omega_{\rm
p}}{\rm exp}\left(\frac{\epsilon-J(\Phi)}{C_{\rm exp}}\right)\,dx.$ (7)

We can show that there exist constants $C,C_{\rm exp}$ such that the condition
$\mathfrak{f}_{\rm pen}^{(1)}(\Phi)\leq C$ implies that
$J(\Phi)\geq\epsilon/2$ for all $x\in\Omega_{\rm p}$ (cf. [32, section 2.2]),
that is $\\{\mathbf{a}:\mathfrak{f}_{\rm
pen}^{(1)}(\texttt{N}(\mathbf{a}))\leq C\\}\subset A_{\rm jac}$. Recalling
Jacobi’s formula, we find that $\|\nabla J(\Phi)\|_{2}\leq
C\|\nabla\Phi\|_{\rm F}\|H(\Phi)\|_{\rm F}$, for some constant $C$ that is
independent of $\nabla\Phi$; therefore, we can replace $\|\nabla
J(\Phi)\|_{L^{\infty}(\Omega_{\rm p})}$ with
$\|H(\Phi)\|_{L^{\infty}(\Omega_{\rm p})}$, to obtain

$\mathfrak{f}_{\rm pen}^{\rm th}(\Phi)=\mathfrak{f}_{\rm
jac}(\Phi)+\|H(\Phi)\|_{L^{\infty}(\Omega_{\rm p})}^{2}.$ (8)

The penalty (8) is defined for elements of the Sobolev space
$W^{2,\infty}(\Omega)=\\{v\in L^{\infty}(\Omega):\nabla v,H(v)\in
L^{\infty}(\Omega)\\}$. We anticipate that the penalty (8) is still not fully
amenable for computations: we address the problem of constructing an
actionable penalty for (3) in section 4.

###  2.2 Mathematical background

We say that $\Omega_{\rm p}\subset\mathbb{R}^{2}$ is a polytope if the
boundary of $\Omega_{\rm p}$, $\partial\Omega_{\rm p}$, consists of a finite
number of flat sides (_faces_). The boundary of two-dimensional polytopes is
described by a finite number of straight segments (_edges_) connected to form
(possibly several disjoint) closed polygonal chains; the points of
$\partial\Omega_{\rm p}$ where two consecutive non-parallel edges meet are
dubbed _vertices_ ; the union of all vertices is here denoted by
$V=\\{x_{i}^{\rm v}\\}_{i=1}^{N_{\rm v}}$. We consider polytopes that satisfy
the condition below: note that Definition 2.1 implies that each vertex $x^{\rm
v}$ of a regular polytope is the intersection of two edges (cf. Figure 2).

######  Definition 2.1.

The bounded polytope $\Omega_{\rm p}$ is said to be regular if there exist
$N+1$ bounded polytopes $\Omega_{{\rm int},1},\ldots,$ $\Omega_{{\rm int},N}$,
$\Omega_{\rm ext}$ such that (i) $\Omega_{\rm p}=\Omega_{\rm
ext}\setminus\bigcup_{i=1}^{N}\Omega_{{\rm int},i}$, (ii) $\Omega_{{\rm
int},1},\ldots,\Omega_{{\rm int},N}$ are pairwise disjoint polytopes that are
compactly embedded in $\Omega_{\rm ext}$ and (iii) $\Omega_{\rm
ext},\Omega_{{\rm int},1},\ldots,\Omega_{{\rm int},N}$ are isomorphic to the
unit ball.

$\Omega_{\rm p}$ (a)

$\Omega_{\rm p}$ (b)

Figure 2: interpretation of Definition 2.1. (a) regular polytope. (b)
irregular polytope.

We denote by $\mathbf{n}(x)$ the outward normal to $\Omega_{\rm p}$ at
$x\in\partial\Omega_{\rm p}$; we also denote by $\overline{\Omega}_{\rm p}$
the closure of $\Omega_{\rm p}$ in $\mathbb{R}^{2}$. Propositions 2.1 and 2.2
contain two important results for mappings in polytopes. The proof of
Proposition 2.1 is technical and is postponed to Appendix A.

######  Proposition 2.1.

Let $\Omega_{\rm p}$ be a regular bounded polytope. Define the space
${\mathfrak{U}}_{0}={\mathfrak{U}}_{0}(\Omega_{\rm p})=\\{\varphi\in
C^{1}(\overline{\Omega}_{\rm
p};\mathbb{R}^{2})\,:\,\varphi\cdot\mathbf{n}|_{\partial\Omega_{\rm p}}=0\\}$
and consider the vector-valued function $\Phi=\texttt{id}+\varphi$ with
$\varphi\in{\mathfrak{U}}_{0}$. Then, $\Phi$ is a bijection in $\Omega_{\rm
p}$ if $\min_{x\in\overline{\Omega}_{\rm p}}J(\Phi)>0$.

######  Proposition 2.2.

Let $\Phi$ be a diffeomorphism in $\Omega_{\rm p}$. Then, $\Phi(V)=V$.
Furthermore, if $\Phi(x^{\rm v})=x^{\rm v}$ for all $x^{\rm v}\in V$, then
$\Phi=\texttt{id}+\varphi$ with $\varphi\in\mathfrak{U}_{0}$.

###### Proof.

Let $x^{\rm v}\in V$ and define $y=\Phi(x^{\rm v})$; recalling Nanson’s
formula [19, page 7], we find that the normal to $\Phi(\Omega_{\rm p})$
satisfies $\mathbf{n}_{\Phi}(y)\propto
J(\Phi(x))(\nabla\Phi(x))^{-T}\mathbf{n}(x);$ therefore, if
$x\in\partial\Omega_{\rm p}\mapsto\mathbf{n}(x)$ is discontinuous at $x^{\rm
v}$, we must have that $x\mapsto\mathbf{n}_{\Phi}(x)$ is discontinuous at $y$.
Recalling the definition of vertices, we conclude that $y\in V$ for any
$x^{\rm v}\in V$, that is $\Phi(V)\subset V$. Since $\Phi$ is a bijection in
$\Omega_{\rm p}$, the image points $\\{\Phi(x^{\rm v}):x^{\rm v}\in V\\}$
should all be distinct and thus ${\rm card}(\Phi(V))={\rm card}(V)$, which
implies $\Phi(V)=V$.

Let $\Phi$ satisfy the condition $x^{\rm v}=\Phi(x^{\rm v})$ for all $x^{\rm
v}\in V$ and define the displacement $\varphi=\Phi-\texttt{id}$. Consider the
edge $F\subset\partial\Omega_{\rm p}$ and consider the parameterization
$\gamma_{\rm f}:[0,1]\to F$ such that $\gamma_{\rm f}(s)=x_{1}^{\rm
v}+s\|x_{2}^{\rm v}-x_{1}^{\rm v}\|_{2}\mathbf{t}_{\rm f}$. Since
$\partial\Omega_{\rm p}$ is closed and $\Phi(\partial\Omega_{\rm
p})=\partial\Omega_{\rm p}$, we either have $\Phi(F)=\partial\Omega_{\rm
p}\setminus F$ or $\Phi(F)=F$. Recalling the expression for the normal
$\mathbf{n}_{\Phi}$, we find that $\mathbf{n}_{\Phi}$ is continuous in
$\Phi(F)$: since $\partial\Omega_{\rm p}\setminus F$ contains at least one
vertex in addition to $x_{1}^{\rm v},x_{2}^{\rm v}$ (we here exploit the fact
that two-dimensional polytopes have at least three vertices), we must have
$\Phi(F)=F$.

The condition $\Phi(F)=F$ implies that $\Phi(\gamma_{\rm f}(s))=x_{1}^{\rm
v}+\alpha(s)\mathbf{t}_{\rm f}$ for all $s\in[0,1]$ and for some injective
function $\alpha:[0,1]\to[0,\|x_{2}^{\rm v}-x_{1}^{\rm v}\|_{2}]$. Given $x\in
F$, since $\mathbf{n}(x)=\mathbf{n}_{\rm f}$ and $\mathbf{n}_{\rm
f}\perp\mathbf{t}_{\rm f}$, we find

$\mathbf{n}(x)\cdot\varphi(x)=\mathbf{n}_{\rm
f}\cdot\left(\Phi(x)-x\right)=\mathbf{n}_{\rm f}\cdot\left(x_{1}^{\rm
v}+\alpha(s)\mathbf{t}_{\rm f}-x_{1}^{\rm v}-s\|x_{2}^{\rm v}-x_{1}^{\rm
v}\|_{2}\mathbf{t}_{\rm f}\right)=0,$

which is the desired result. ∎

###  2.3 Construction of finite-dimensional operators for registration

We introduce the finite-dimensional space $\mathcal{U}$ spanned by
$\\{\varphi_{i}\\}_{i=1}^{M}\subset{\mathfrak{U}}_{0}$, and the affine space
$\texttt{id}+\mathcal{U}$. We introduce the operator $\texttt{N}_{\rm
p}:\mathbb{R}^{M}\to\texttt{id}+\mathcal{U}$ such that $\texttt{N}_{\rm
p}(\mathbf{a})=\texttt{id}+\sum_{i=1}^{M}(\mathbf{a})_{i}\varphi_{i}$.
Exploiting Proposition 2.1 and the discussion in section 2.1, we show in
Corollary 2.3 that the operator $\texttt{N}_{\rm p}$ and the penalty (8)
satisfy (2): the stronger regularity assumption is required by the choice of
the penalty. The proof is straightforward and is here omitted.

######  Corollary 2.3.

Let $\mathcal{U}={\rm span}\\{\varphi_{i}\\}_{i=1}^{M}$ be an $M$-dimensional
subspace of ${\mathfrak{U}}_{0}\cap W^{2,\infty}(\Omega_{\rm p})$. Then, the
affine operator $\texttt{N}_{\rm p}$ and the penalty
$\mathbf{a}\mapsto\mathfrak{f}_{\rm pen}^{\rm th}(\texttt{N}_{\rm
p}(\mathbf{a}))$ (cf. (8)) satisfy (2).

###  2.4 Approximation properties of affine mappings

Proposition 2.2 shows that bijections are of the form
$\Phi=\texttt{id}+\varphi$ with $\varphi\in{\mathfrak{U}}_{0}$ if $\Phi(x^{\rm
v})=x^{\rm v}$ for all $x^{\rm v}\in V$ (i.e., $\Phi|_{V}=\texttt{id}$). Since
$\varphi|_{V}=0$ for all $\varphi\in{\mathfrak{U}}_{0}$, diffeomorphisms
$\Phi$ that do not satisfy $\Phi|_{V}=\texttt{id}$ do not belong to
$\texttt{id}+{\mathfrak{U}}_{0}$. It is easy to construct bijections that do
not fulfill the requirement $\Phi|_{V}=\texttt{id}$: to provide a concrete
example, consider the map $\Phi(x)=-x$ for the polytope $\Omega_{\rm
p}=(-1,1)^{2}$. However, in the setting of MOR and also geometry reduction, we
are interested in parametric maps that are smooth deformations of the identity
map: since registration is applied with respect to an element of the solution
manifold $\mathcal{M}=\\{u_{\mu}:\mu\in\mathcal{P}\\}$, we can indeed assume
that there exists $\mu\in\mathcal{P}$ such that $\Phi_{\mu}=\texttt{id}$. We
hence have the following result, which follows from Proposition 2.2.

######  Corollary 2.4.

Let $\Phi:\Omega_{\rm p}\times\mathcal{P}\to\Omega_{\rm p}$ be a continuous
function of the parameter $\mu\in\mathcal{P}$ and let $\Phi_{\mu}$ be a
diffeomorphism for all $\mu\in\mathcal{P}$. Assume that
$\Phi_{\mu^{\prime}}=\texttt{id}$ for some $\mu^{\prime}\in\mathcal{P}$ and
that $\mathcal{P}$ is simply connected. Then,
$\Phi_{\mu}=\texttt{id}+\varphi_{\mu}$ with
$\varphi_{\mu}\in{\mathfrak{U}}_{0}$ for all $\mu\in\mathcal{P}$.

###### Proof.

Exploiting the first statement of Proposition 2.2, we find that
$\Phi_{\mu}(V)=V$ for all $\mu\in\mathcal{P}$. Since the set of vertices $V$
is a discrete closed set, if the function $\Phi:\Omega_{\rm
p}\times\mathcal{P}\to\Omega_{\rm p}$ is continuous with respect to the
parameter $\mu\in\mathcal{P}$ and is bijective in $\Omega_{\rm p}$ for all
$\mu\in\mathcal{P}$, then the hypothesis
$\Phi_{\mu^{\prime}}|_{V}=\texttt{id}$ _for some_ $\mu^{\prime}\in\mathcal{P}$
implies $\Phi_{\mu}|_{V}=\texttt{id}$ _for all_ $\mu\in\mathcal{P}$.
Therefore, exploiting the second statement of Proposition 2.2, we conclude
that $\Phi_{\mu}=\texttt{id}+\varphi_{\mu}$ with
$\varphi_{\mu}\in{\mathfrak{U}}_{0}$ for all $\mu\in\mathcal{P}$. ∎

Given the parametric registration problem (3) with $\texttt{N}=\texttt{N}_{\rm
p}$, we denote by $\mathbf{a}_{\mu}^{\rm opt}$ a solution to (3) for any
$\mu\in\mathcal{P}$ and we define the corresponding manifold
$\mathcal{M}_{\Phi}:=\\{\mathbf{a}_{\mu}^{\rm
opt}:\mu\in\mathcal{P}\\}\subset\mathbb{R}^{M}$. If $\mathcal{M}_{\Phi}$ is
reducible — that is, there exists an orthogonal matrix
$\mathbf{W}\in\mathbb{R}^{M\times m}$ whose columns approximate the elements
of $\mathcal{M}_{\Phi}$ for $m\ll M$ — we can restrict the search space in (3)
to the image of $\mathbf{W}$, ${\rm col}(\mathbf{W})$. As discussed in the
introduction, since the pair $(\texttt{N}(\cdot)=\texttt{N}_{\rm
p}(\mathbf{W}\cdot),\mathfrak{f}_{\rm pen}^{\rm th}(\texttt{N}_{\rm
p}(\mathbf{W}\cdot))$ satisfies (2) for any full-rank matrix $\mathbf{W}$, the
choice of the latter can be made solely based on approximation considerations.

###  2.5 Approximation of Lipschitz maps

The proof of Proposition 2.1 exploits the Hadamard’s global inverse function
theorem (cf. [15, Theorem 6.2.8]); it hence relies on the assumption that
$\Phi$ is of class $C^{1}$. Our numerical investigations suggest that a
similar result — possibly with further conditions — might hold for piecewise-
smooth maps such as FE fields.

On the other hand, Proposition 2.2 does not hold for general Lipschitz maps,
as shown in the next example. First, we introduce the square domain
$\Omega=(0,1)^{2}$, and we define the map
$\Phi_{1}:\Omega=(0,1)^{2}\to\mathbb{R}^{2}$ such that

$\Phi_{1}(x)=\left[\begin{array}[]{cc}1/2&0\\\ -1/2&1\\\
\end{array}\right]x\mathbbm{1}_{x_{1}<x_{2}}(x)+\left[\begin{array}[]{cc}1&-1/2\\\
0&1/2\\\ \end{array}\right]x\mathbbm{1}_{x_{1}\geq x_{2}}(x).$ (9)

It is easy to verify that $\Phi_{1}$ is a Lipschitz map from $\Omega$ in the
unit triangle $D=\\{x\in\Omega:x_{1}+x_{2}<1\\}$ with Lipschitz inverse (cf.
Figure 3(a)). Then, we introduce a smooth bijection $\Phi_{2}:D\to D$ that
maps $x^{\star}=[1/2,1/2]$ into $y^{\star}=[1/4,3/4]$ (cf. Figure 3(b)):
$\Phi_{2}$ can be constructed using an expansion of quadratic polynomials; we
omit the explicit expression. Finally, we define the map
$\Phi=\Phi_{1}^{-1}\circ\Phi_{2}\circ\Phi_{1}$: clearly, $\Phi$ is a bijection
in $\Omega$; furthermore,
$\Phi([1,1])=\Phi_{1}^{-1}\left(\Phi_{2}([1/2,1/2])\right)=\Phi_{1}^{-1}\left([1/4,3/4]\right)=[1,1/2].$
We hence found a bijection $\Phi$ in $\Omega$ such that $\Phi(V)\neq V$.

$\Omega$$D$${\Phi_{1}}$ (a)

$x^{\star}$$y^{\star}$$\Omega$$D$${\Phi_{2}}$ (b)

Figure 3: example of Lipschitz map with $\Phi(V)\neq V$. (a) Lipschitz map
from the unit square to the unit triangle (cf (9)). (b) smooth map in the unit
triangle $D$.

##  3 Nonlinear maps for general domains

We address the extension to non-polytopal domains. First, we justify the need
for nonlinear ansätze; then, we investigate the approximation properties of
the nonlinear ansätze (4) and (5).

###  3.1 Inadequacy of affine maps for domains with curved boundaries

For many problems of interest, including the one considered in the numerical
results of section 5, it is of paramount importance to deform points that lie
on curved edges or faces; we hence seek an operator N and a penalty
$\mathfrak{f}_{\rm pen}$ that satisfy (2) and enable non-trivial deformations
of points on curved boundaries. Next result, which generalizes [34, Lemma
2.1], shows that non-trivial affine mappings for curved domains lead to
admissible sets $A$ with empty interior: therefore, affine maps in combination
with the penalty (8) do not satisfy (2). This observation shows the inadequacy
of affine maps for domains with curved boundaries and ultimately justifies the
need for a nonlinear ansatz.

######  Lemma 3.1.

Let $\Omega\subset\mathbb{R}^{2}$ be a domain with curved boundary
$\Gamma\subset\partial\Omega$; let $\texttt{N}:\mathbb{R}^{M}\to{\rm
Lip}(\overline{\Omega};\mathbb{R}^{2})$ be an affine operator of the form
$\texttt{N}(\mathbf{a})=\texttt{id}+\sum_{i=1}^{M}(\mathbf{a})_{i}\varphi_{i}$
and let $A\subset\mathbb{R}^{M}$ be the set of admissible (i.e., bijective in
$\Omega$) maps. Suppose that there exists $i\in\\{1,\ldots,M\\}$ and
$x\in\Gamma$ such that $\varphi_{i}(x)\neq 0$. Then, $0\in\mathbb{R}^{M}$ does
not belong to the interior of $A$. Furthermore, if $\Gamma=\partial\Omega$,
then $A$ has empty interior.

###### Proof.

By contradiction, there exists $r>0$ such that $\mathcal{B}_{r}(0)\subset A$;
therefore, there exists $\delta>0$ such that
$\Phi_{t}=\texttt{id}+t\varphi_{i}$ is a bijection in $\Omega$ for all
$t\in(-\delta,\delta)$. Since $\Phi_{t}$ is bijective in $\Omega$, we must
have $\Phi_{t}(\partial\Omega)=\partial\Omega$, which implies that the linear
segment $S_{\delta}=\\{x+t\varphi_{i}(x):t\in(-\delta,\delta)\\}$ is contained
in $\partial\Omega$. We conclude that $\partial\Omega$ is flat in the
neighborhood of $x\in\Gamma$. Contradiction.

We remark that the argument of the proof solely relies on the fact that
$\Phi_{t=0}(x)\in\Gamma$: if $\partial\Omega$ is entirely curved (that is, if
$\Gamma=\partial\Omega$) we can replicate the same argument for any
$\mathbf{a}\in A$ to prove that $A$ has empty interior. ∎

###  3.2 Approximation properties of compositional maps

Given the domain $\Omega$, we introduce the polytope $\Omega_{\rm p}$ and the
isomorphism $\Psi$ from $\Omega_{\rm p}$ to $\Omega$; we denote by
$V=\\{x_{i}^{\rm v}\\}_{i=1}^{N_{\rm v}}$ the vertices of $\Omega_{\rm p}$ and
we denote by $\\{F_{j}^{\rm p}\\}_{j=1}^{N_{\rm f,p}}$ the facets of
$\Omega_{\rm p}$; we also define the mapped facets $\\{F_{j}=\Psi(F_{j}^{\rm
p})\\}_{j=1}^{N_{\rm f,p}}$. In the remainder, we assume that the pair
$(\Omega_{\rm p},\Psi)$ satisfies the following assumption; the set $V_{\rm
ang}$ of angular points of $\Omega$ is the subset of $\partial\Omega$ where
the normal $\mathbf{n}$ is discontinuous.

######  Hypothesis 1.

The pair $(\Omega_{\rm p},\Psi)$ satisfies the following conditions: (i) the
vertices $V$ of $\Omega_{\rm p}$ belong to $\partial\Omega$; (ii) the set
$V_{\rm ang}$ of angular points of $\Omega$ is contained in $V$; (iii) the
linear facets of $\partial\Omega$ belong to $\partial\Omega_{\rm p}$; (iv) the
field $\Psi$ is a Lipschitz bijection from $\Omega_{\rm p}$ to $\Omega$ that
satisfies $\Psi|_{V}=\texttt{id}$ and $J(\Psi)>0$ a.e..

By construction, any bijection $\Phi$ of the form (4) — that is,
$\Phi=\Psi\circ\Phi_{\rm p}\circ\Psi^{-1}$ with $\Phi_{\rm
p}\in\texttt{id}+{\mathfrak{U}}_{0}(\Omega_{\rm p})$— satisfies
$\Phi|_{V}=\texttt{id}$ and $\Phi(F_{j})=F_{j}$ for $j=1,\ldots,N_{\rm f}$.
Conversely, we have the following result, which signifies that nonlinear maps
of the form (4) are dense in the space of diffeomorphisms that satisfy
$\Phi|_{V}=\texttt{id}$.

######  Lemma 3.2.

Let $(\Omega_{\rm p},\Psi)$ satisfy Hypothesis 1. Let $\Phi$ be a
diffeomorphism in $\Omega$ such that $\Phi|_{V}=\texttt{id}$. Then,
$\widetilde{\Phi}=\Psi^{-1}\circ\Phi\circ\Psi=\texttt{id}+\varphi$ for some
$\varphi\in{\rm Lip}(\Omega_{\rm p};\mathbb{R}^{2})$ such that
$\varphi\cdot\mathbf{n}|_{\partial\Omega_{\rm p}}=0$.

###### Proof.

Since $\Phi$ is a diffeomorphism, we must have
$\Phi(\partial\Omega)=\partial\Omega$; furthermore, the condition
$\Phi|_{V}=\texttt{id}$ implies that $\Phi$ preserves the orientation of the
boundary. In conclusion, we find ${\Phi}(F_{j})=F_{j}$ for $j=1,\ldots,N_{\rm
f}$. Therefore, we have that $\widetilde{\Phi}=\Psi^{-1}\circ\Phi\circ\Psi$ is
a bijection in $\Omega_{\rm p}$ that satisfies $\widetilde{\Phi}(F_{j}^{\rm
p})=F_{j}^{\rm p}$ for $j=1,\ldots,N_{\rm f}$. Then, exploiting the same
argument of the proof of Proposition 2.2, we find that
$\widetilde{\Phi}=\texttt{id}+\varphi$ with
$\varphi\cdot\mathbf{n}|_{\partial\Omega_{\rm p}}=0$. We omit the details. ∎

###  3.3 Multi-layer compositional maps

Lemma 3.2 shows that single-layer compositional maps cannot approximate
arbitrary diffeomorphisms in curved domains. To address this issue, we might
consider more general multi-layer maps of the form (5). We here study the
approximation properties of (5) for $\ell=2$ layers

$\texttt{N}(\mathbf{a}_{1},\mathbf{a}_{2})=\texttt{N}_{1}(\mathbf{a}_{1})\circ\texttt{N}_{2}(\mathbf{a}_{2}),\quad\texttt{N}_{i}(\mathbf{a}_{i})=\Psi_{i}\circ\texttt{N}_{{\rm
p},i}(\mathbf{a}_{i})\circ\Psi_{i}^{-1},\;\;i=1,2.$

The extension of the analysis to arbitrary number of layers is the subject of
ongoing investigations.

We denote by $\Omega_{{\rm p},1}$ and $\Omega_{{\rm p},2}$ the polytopes
associated to the bijections $\Psi_{1},\Psi_{2}$; we assume that
$(\Omega_{{\rm p},1},\Psi_{1})$ and $(\Omega_{{\rm p},2},\Psi_{2})$ satisfy
Hypothesis 1. We also denote by $V_{i}^{\rm fict}=\\{x_{j,i}^{\rm
v}\\}_{j=1}^{N_{{\rm v},i}}$ the _fictitious_ vertices of $\Omega_{{\rm
p},i}$, that is, the vertices of $\Omega_{{\rm p},i}$ that do not correspond
to angular points of $\partial\Omega$. Finally, given $x,y\in\partial\Omega$,
we define the set of curves

$A_{\partial\Omega}(x,y)=\left\\{\gamma\in
C^{1}([0,1];\partial\Omega):\gamma(0)=x,\;\gamma(1)=y,\;\|\dot{\gamma}\|_{2}\equiv{\rm
const}\;\right\\}$ (10a) and the geodesic distance ${\rm
dist}_{\partial\Omega}(x,y)=\left\\{\begin{array}[]{ll}\displaystyle{\inf\left\\{\int_{0}^{1}\|\dot{\gamma}\|_{2}\,dt:\;\;\gamma\in
A_{\partial\Omega}(x,y)\right\\}}&{\rm
if}\;A_{\partial\Omega}(x,y)\neq\emptyset\\\\[8.53581pt] +\infty&{\rm
if}\;A_{\partial\Omega}(x,y)=\emptyset\\\ \end{array}\right.$ (10b)

Note that $C^{1}$ curves cannot pass through a non-smooth point of the
boundary (angular point): therefore, if $x,y\in\partial\Omega$ are separated
by angular points (e.g., points on two separate faces of a polytope), we have
$A_{\partial\Omega}(x,y)=\emptyset$ and ${\rm
dist}_{\partial\Omega}(x,y)=+\infty$.

Next Lemma shows that two-layer maps can approximate arbitrary diffeomorphisms
under the assumption of small deformations. Note that boundary deformations
should be smaller than (i) the minimum distance between the fictitious
vertices of $\Omega_{{\rm p},1}$ and $\Omega_{{\rm p},2}$ — $\min_{x\in
V_{1}^{\rm fict},y\in V_{2}^{\rm fict}}{\rm dist}_{\partial\Omega}(x,y)$; (ii)
the minimum distance between two vertices of $\Omega_{{\rm p},1}$ and
$\Omega_{{\rm p},2}$ — $\min_{x,y\in V_{i},x\neq y}{\rm
dist}_{\partial\Omega}(x,y)$ with $i=1,2$.

######  Lemma 3.3.

Let $(\Omega_{{\rm p},1},\Psi_{1})$ and $(\Omega_{{\rm p},2},\Psi_{2})$
satisfy Hypothesis 1. Let $\Phi$ be a diffeomorphism in $\Omega$ with positive
Jacobian determinant such that $\max_{x\in\partial\Omega}{\rm
dist}_{\partial\Omega}$ $(x,\Phi(x))$ $<C$ where $C=C(\Omega_{{\rm
p},1},\Omega_{{\rm p},2})$ is given by

$\begin{array}[]{rl}C:=\min\Big{\\{}&\displaystyle{\min_{x\in V_{1}^{\rm
fict},y\in V_{2}^{\rm fict}}{\rm dist}_{\partial\Omega}(x,y),\;\;\min_{x,y\in
V_{1},x\neq y}{\rm dist}_{\partial\Omega}(x,y),}\\\\[8.53581pt]
&\displaystyle{\min_{x,y\in V_{2},x\neq y}{\rm
dist}_{\partial\Omega}(x,y)\Big{\\}}.}\\\ \end{array}$ (11)

where $V_{1}^{\rm fict},V_{2}^{\rm fict}$ are the vertices of $\Omega_{{\rm
p},1},\Omega_{{\rm p},2}$ that do not correspond to angular points of
$\partial\Omega$. Then, there exist $\Phi_{1}=\texttt{id}+\varphi_{1}$ and
$\Phi_{2}=\texttt{id}+\varphi_{2}$ such that
$\Phi=\Psi_{1}\circ\Phi_{1}\circ\Psi_{1}^{-1}\circ\Psi_{2}\circ\Phi_{2}\circ\Psi_{2}^{-1}$
and $\varphi_{i}\cdot\mathbf{n}|_{\partial\Omega_{{\rm p},i}}=0$ for $i=1,2$.

###### Proof.

We assume that $\Omega$ has only one curved boundary
$\Gamma\subset\partial\Omega$ with parameterization
$\boldsymbol{\gamma}:[0,1]\to\Gamma$. We define $N_{\rm v}=N_{{\rm
v},1}+N_{{\rm v},2}$, $\\{y_{\ell}\\}_{\ell=1}^{N_{\rm v}}\subset\Gamma$ and
$\\{t_{\ell}\\}_{\ell=1}^{N_{\rm v}}\subset(0,1)$ such that

$\\{y_{\ell}=\boldsymbol{\gamma}(t_{\ell})\\}_{\ell=1}^{N_{\rm v}}=V_{1}^{\rm
fict}\cup V_{2}^{\rm fict},\quad{\rm with}\;\;0<t_{1}<\ldots<t_{N_{\rm v}}<1.$

Definition (11) implies that ${\rm
dist}_{\partial\Omega}(y_{\ell},y_{{\ell+1}})\geq C$, while the hypothesis
$\max_{x\in\partial\Omega}$ ${\rm dist}_{\partial\Omega}$ $(x,\Phi(x))$ $<C$
implies that

$\Phi(y_{\ell})\in\left\\{\boldsymbol{\gamma}(t):t\in(t_{\ell-1},t_{\ell+1})\right\\},\quad{\rm
for}\;\ell=1,\ldots,N_{\rm v},$ (12)

with $t_{\ell=0}=0$, $t_{\ell=N_{\rm v}+1}=1$.

Exploiting (12) and the fact that $\Phi$ does not deform angular points of
$\partial\Omega$ (cf. Proposition 2.2), we find that there exists a
diffeomorphism $\widetilde{\Phi}_{1}$ such that $\widetilde{\Phi}_{1}(y)=y$
for all $y\in V_{1}$ and $\widetilde{\Phi}_{1}(y)=\Phi(y)$ for all $y\in
V_{2}$. Exploiting Lemma 3.2, we find that
$\widetilde{\Phi}_{1}=\Psi_{1}\circ\Phi_{1}\circ\Psi_{1}^{-1}$ for some
$\Phi_{1}\in\texttt{id}+{\mathfrak{U}}_{0}(\Omega_{{\rm p},1})$. Similarly,
since $\widetilde{\Phi}_{1}^{-1}\circ\Phi(y)=y$ for all $y\in V_{2}$, we find
that
$\widetilde{\Phi}_{1}^{-1}\circ\Phi=:\widetilde{\Phi}_{2}=\Psi_{2}\circ\Phi_{2}\circ\Psi_{2}^{-1}$
for some $\Phi_{2}\in\texttt{id}+{\mathfrak{U}}_{0}(\Omega_{{\rm p},2})$. In
conclusion, we find that $\Phi=\widetilde{\Phi}_{1}\circ\widetilde{\Phi}_{2}$
which is the desired result. ∎

The extension to multiple layers enables the approximation of diffeomorphisms
that involve larger deformations over curved edges. Figure 4 illustrates the
approximation power of multi-layer maps for deformations over an airfoil:
Figure 4(a) shows the two polytopes $\Omega_{\rm p,1},\Omega_{\rm p,2}$ that
are used to define the nonlinear ansatz and two points $x,y$ on the profile
such that $\Phi(x)=y$; Figures 4(b)-(c)-(d) show the action of the maps
$\texttt{N}_{3}$, $\texttt{N}_{2}$ and $\texttt{N}_{1}$, respectively. In more
detail, the map $\texttt{N}_{3}$ — which is associated to $\Omega_{\rm p,1}$ —
maps the point $x$ into the point $x_{1}$ (cf. Figure 4(b)); the map
$\texttt{N}_{2}$ — which is associated to $\Omega_{\rm p,2}$ — deforms the
point $x_{1}$ into the point $x_{2}$ (cf. Figure 4(c)); the map
$\texttt{N}_{1}$ — which is associated to $\Omega_{\rm p,1}$ — deforms the
point $x_{2}$ into the point $y$ (cf. Figure 4(d)). In conclusion, the
considered map satisfies
$\texttt{N}_{1}(\texttt{N}_{2}(\texttt{N}_{3}(x)))=y$.

$x$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\Omega_{\rm
p,2}}$$y$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Omega_{\rm
p,1}}$ (a)

$x$$y$$x_{1}$ (b) $\texttt{N}_{3}$

$x$$y$$x_{1}$$x_{2}$ (c) $\texttt{N}_{2}$

$x$$y$$x_{2}$ (d) $\texttt{N}_{1}$

Figure 4: approximation of large deformations over curved boundaries using a
multi-layer ($\ell=3$) compositional map.

###  3.4 Discussion

The analysis of this section shows that compositional maps (4) and their
multi-layer generalization (5) can be employed for registration in curved
domains. Since $\Psi$ in (4) and $\Psi_{1},\ldots,\Psi_{\ell}$ in (5) are
independent of the coefficients $\mathbf{a}$, we can apply the penalty term
(8) to the polytopal map $\texttt{N}_{\rm p}(\mathbf{a})$ in (4), or
$\\{\texttt{N}_{{\rm p},i}(\mathbf{a}_{i}\\}_{i=1}^{\ell}$ in (5); then, we
can exploit Proposition 2.1 (see also the discussion in section 2.3) to prove
that (4) and (5) satisfy (2).

Lemmas 3.2 and 3.3 investigate the approximation properties of the ansätze (4)
and (5): the analysis shows that multi-layer maps provide much more
approximation power, even for moderate number of layers $\ell$. We note,
however, that multi-layer maps are considerably more challenging to implement
and might also be significantly more expensive to evaluate: the solution to
(3) requires indeed many evaluations of the mapping $\texttt{N}(\mathbf{a})$
and its gradient, which involve the evaluation of the bijections
$\\{\Psi_{i}\\}_{i=1}^{\ell}$ and their inverses. In this work, we focus on
the implementation of registration strategies based on the more elementary
ansatz (4) and we refer to a future work for the implementation of
registration methods based on the model class (5).

##  4 Methodology

Given the domain $\Omega\subset\mathbb{R}^{2}$, we present a registration
procedure based on compositional mappings of the form (4). We recall (cf. (3))
that we consider the minimization problem

$\min_{\mathbf{a}\in\mathbb{R}^{M}}\mathfrak{f}_{\mu}^{\rm
tg}(\texttt{N}(\mathbf{a}))\;+\xi\;\mathfrak{f}_{\rm pen}(\mathbf{a}),$

where $\texttt{N}(\mathbf{a})=\Psi\circ\texttt{N}_{\rm
p}(\mathbf{a})\circ\Psi^{-1}$, $\Psi:\Omega_{\rm p}\to\Omega$ is a bijection
from the polytope $\Omega_{\rm p}$ to $\Omega$ and $\texttt{N}_{\rm
p}(\mathbf{a}):\Omega_{\rm p}\to\Omega_{\rm p}$ satisfies $\texttt{N}_{\rm
p}(\mathbf{a})=\texttt{id}+\sum_{i=1}^{M}(\mathbf{a})_{i}\varphi_{i}$ with
$\varphi_{i}\cdot\mathbf{n}|_{\partial\Omega_{\rm p}}=0$ for $i=1,\ldots,M$.
In section 4.1, we discuss the construction of the polytope $\Omega_{\rm p}$
and the mapping $\Psi$ in (4) based on a curved mesh of $\Omega$; in section
4.2, we introduce the FE space $\mathcal{U}={\rm span}\\{\varphi_{i}\\}_{i}$
for the displacement field in the polytope $\Omega_{\rm p}$ and we discuss the
choice of the functional norm $\|\cdot\|$ for $\mathcal{U}$ that enters in the
penalty $\mathfrak{f}_{\rm pen}$; then, in section 4.3, we discuss the
definition of the curved mesh of $\Omega$. In sections 4.4, 4.5, and 4.6, we
review the choice of the penalty $\mathfrak{f}_{\rm pen}$ and target functions
$\mathfrak{f}_{\mu}^{\rm tg}$ proposed in [32, 34] and the extension to
parametric problems. To simplify the presentation, we here assume that the
domain $\Omega$ is parameter-independent.

We recall that our ultimate goal is to devise Lagrangian approximations of
parametric fields. Given the parametric field
$u:\Omega\times\mathcal{P}\to\mathbb{R}^{D_{\rm u}}$ , Lagrangian
approximations read as

$u_{\mu}\approx\widehat{u}_{\mu}:=\widetilde{u}_{\mu}\circ\Phi_{\mu}^{-1},$
(13)

where $\widetilde{u}_{\mu}$ is a low-rank approximation of the mapped field
${u}_{\mu}\circ\Phi_{\mu}$. If we denote by $\mathcal{T}_{\rm pb}$ the high-
fidelity (HF) mesh used to approximate the mapped field, it might be important
(see, e.g., [3]) to ensure that the deformed mesh $\Phi_{\mu}(\mathcal{T}_{\rm
pb})$, which shares the connectivity matrix with $\mathcal{T}_{\rm pb}$ but
has deformed nodes according to $\Phi_{\mu}$, is well-behaved for all
$\mu\in\mathcal{P}$: in section 4.4, we discuss how to enforce this constraint
in the registration procedure.

###  4.1 Definition of the Lipschitz map $\Psi$ based on a curved mesh of
$\Omega$

Our point of departure is the definition of a triangular curved mesh
$\mathcal{T}$ of degree $\kappa>0$ of the domain $\Omega$: the mesh
$\mathcal{T}$ is uniquely identified by the mesh nodes and the connectivity
matrix. We define the reference (or master) element
$\widehat{\texttt{D}}=\\{x\in(0,1)^{2}:\sum_{i=1}^{2}(x)_{i}<1\\}$ with
reference vertices and nodes $\\{\tilde{x}_{i}^{\rm
p}\\}_{i=1}^{3}\subset\\{\tilde{x}_{i}\\}_{i=1}^{n_{\rm
lp}}\subset\overline{\widehat{\texttt{D}}}$, and the associated Lagrangian
bases $\\{\ell_{i}^{\rm fe,p}\\}_{i=1}^{3}$ and $\\{\ell_{i}^{\rm
fe}\\}_{i=1}^{n_{\rm lp}}$ for polynomials of degree one, $\mathbb{P}_{1}$,
and degree $\kappa$, $\mathbb{P}_{\kappa}$. We denote by
$\\{\texttt{D}_{k}\\}_{k=1}^{N_{\rm e}}$ the elements of the mesh, by
$\\{\texttt{F}_{j}\\}_{j=1}^{N_{\rm f}}$ the facets of the mesh; for each
element $\texttt{D}_{k}$ of the mesh $\mathcal{T}$, we define the nodes
$\\{x_{i,k}^{\rm hf}\\}_{i,k}$ such that $x_{i,k}^{\rm hf}$ is the $i$-th node
of the $k$-th element of the mesh $\mathcal{T}$ for $i=1,\ldots,n_{\rm lp}$
and $k=1,\ldots,N_{\rm e}$; similarly, we define the vertices $\\{x_{j,k}^{\rm
hf,v}\\}_{j,k}$ of the curved triangles for $j=1,\ldots,3$ and
$k=1,\ldots,N_{\rm e}$.

We have now the elements to define the polytope $\Omega_{\rm p}$ and the
mapping $\Psi$. We introduce the FE mappings

$\Psi_{k}(\tilde{x})=\sum_{i=1}^{n_{\rm lp}}x_{i,k}^{\rm hf}\ell_{i}^{\rm
fe}(\tilde{x}),\quad\Psi_{k,\rm p}(\tilde{x})=\sum_{i=1}^{3}x_{i,k}^{\rm
hf,v}\ell_{i}^{\rm fe,p}(\tilde{x}),\quad k=1,\ldots,N_{\rm e};$ (14a) for
each $k=1,\ldots,N_{\rm e}$, $\Psi_{k}$ maps the reference element into the
$k$-th element of the mesh, while $\Psi_{k,\rm p}$ is a linear map whose image
defines the $k$-th linearized element of the mesh, $\texttt{D}_{k,\rm
p}:=\Psi_{k,\rm p}(\widehat{\texttt{D}})$. Then, we define the polytope
$\Omega_{\rm p}$ such that $\overline{\Omega}_{\rm p}:=\bigcup_{i=1}^{N_{\rm
e}}\overline{\texttt{D}}_{k,\rm p},$ (14b) the linear triangular mesh
$\mathcal{T}_{\rm p}$ with elements $\\{\texttt{D}_{k,\rm
p}=\Psi^{-1}(\texttt{D}_{k})\\}_{k=1}^{N_{\rm e}}$ and facets
$\\{\texttt{F}_{j,\rm p}=\Psi^{-1}(\texttt{F}_{j})\\}_{j=1}^{N_{\rm f}}$, and
the FE space $\mathcal{X}_{\rm hf,p}=\left\\{\phi\in C(\Omega_{\rm
p}):\phi|_{\texttt{D}_{k,\rm p}}\in\mathbb{P}_{\kappa}\right\\}.$ (14c)
Finally, we introduce the field $\Psi$ from $\Omega_{\rm p}$ to $\Omega$ as
$\Psi:\Omega_{\rm p}\to\Omega\quad{\rm
s.t.}\quad\Psi(x)\big{|}_{\texttt{D}_{k,\rm p}}=\Psi_{k}\left(\Psi_{k,\rm
p}^{-1}(x)\right).$ (14d)

Since $\Psi_{k,\rm p}$ is a linear map, $\Psi\in\mathcal{X}_{\rm hf,p}$. If
the local FE mappings (14a) are invertible with positive Jacobian determinant,
we find that the restriction of $\Psi$ to $\texttt{D}_{k,\rm p}$ is a
bijection from $\texttt{D}_{k,\rm p}$ to $\texttt{D}_{k}$ with positive
Jacobian determinant (composition of bijections with positive Jacobian
determinant), for $k=1,\ldots,N_{\rm e}$: therefore, $\Psi$ is a bijection
from $\Omega_{\rm p}$ to $\Omega$. We further observe that by construction the
restriction of $\Psi$ to the vertices of $\mathcal{T}$ and to linear facets of
the mesh is the identity map: therefore, the pair $(\Omega_{\rm p},\Psi)$
satisfies Hypothesis 1.

Figure 5 provides a graphical interpretation of our construction for a curved
mesh with two elements. Figure 5(a) shows a curved mesh of the domain
$\Omega$, while Figure 5(b) shows the corresponding linearized mesh. The
mapping $\Psi$ is the unique polynomial map of degree $\kappa$ that maps each
element $\texttt{D}_{k,\rm p}$ into the corresponding curved element
$\texttt{D}_{k}$. Notice that by construction the vertices of the mesh
$\mathcal{T}$ are unchanged in $\mathcal{T}_{\rm p}$; notice also that linear
elements of $\mathcal{T}$ are preserved by $\mathcal{T}_{\rm p}$.

$\Omega$$\mathcal{T}$$\texttt{D}_{1}$$\texttt{D}_{2}$ (a)

$\Omega_{\rm p}$$\mathcal{T}_{\rm p}$$\texttt{D}_{1,\rm p}$$\texttt{D}_{2,\rm
p}$ (b)

Figure 5: construction of the polytope $\Omega_{\rm p}$ and the mapping $\Psi$
for a curved mesh with two elements.

The construction of $\Psi$ outlined above might fail for excessively coarse
meshes, as shown in the example of Figure 6. Given the domain $\Omega$
depicted in Figure 6(a), we might be tempted to consider the coarse mesh in
Figure 6(b): even if the curved mesh $\mathcal{T}$ associated with this
partition is a proper mesh of $\Omega$, we find that the resulting polytope
$\Omega_{\rm p}$ is not isomorphic to $\Omega$ up to the boundary (since
$\overline{\Omega}_{\rm p}$ is a rectangle) and $\Psi$ is not continuous up to
the boundary of ${\Omega}_{\rm p}$. To fix this issue for this geometry, we
should consider meshes with at least three points on the profile (cf. Figure
6(c)).

We here rely on polynomial discretizations of very high-order (up to
$\kappa=10$): to ensure accurate and stable computations, we rely on a nodal-
based discretization that exploits Koornwinder polynomials to represent the
local shape functions and to a tensor product of Gauss and Gauss–Radau
quadratures (see, e.g., [10, 14]). Instead, we consider regular (tensorized)
nodes $\\{\tilde{x}_{i}\\}_{i=1}^{n_{\rm
lp}}\subset\overline{\widehat{\texttt{D}}}$; alternative selections that
improve the conditioning of the interpolation matrix can be found in [10, 14],
see also [9].

In our framework, the use of high-order discretizations is motivated by two
independent considerations. First, coarse-grained meshes with relatively few
elements enable rapid searches over the elements and hence guarantee rapid
evaluations of the mapping $\Psi$ for new points in $\Omega_{\rm p}$ — this
feature is crucial in optimization-based registration as iterative algorithms
for (3) require multiple evaluations of (4) for many values of $\mathbf{a}$.
Second, since by construction $\texttt{N}(x^{\rm v};\mathbf{a})=x^{\rm v}$ for
all vertices $x^{\rm v}$ of $\Omega_{\rm p}$ and any
$\mathbf{a}\in\mathbb{R}^{M}$ (cf. section 3.2), the reduction of the number
of vertices on curved facets of $\partial\Omega$ improves the expressiveness
of the ansatz (4).

(a)

(b)

(c)

Figure 6: isolated profile in rectangular domain. (a) domain $\Omega$. (b)
inadmissible coarse-grained partition. (c) admissible coarse-grained
partition.

###  4.2 Finite element displacement space on polytopes

Given the polytope $\Omega_{\rm p}\subset\mathbb{R}^{2}$ and the linear
triangular mesh $\mathcal{T}_{\rm p}$ defined in section 4.1, we denote by
$\texttt{I}_{\rm int}\subset\\{1,\ldots,N_{\rm f}\\}$ the indices of the
internal facets and we define the union of internal facets $\mathcal{E}_{\rm
p}^{\rm int}=\bigcup_{j\in\texttt{I}_{\rm int}}\texttt{F}_{j,\rm p}$. We
further introduce the positive111The positive normal is chosen arbitrarily for
internal facets and coincides with the outward normal to $\Omega_{\rm p}$ for
boundary facets. normal $\mathbf{n}^{+}:\mathcal{E}_{\rm p}^{\rm
int}\to\mathbb{S}^{1}$. Given the field $\phi\in\mathcal{X}_{\rm hf,p}$ and an
arbitrary (possibly differential) operator $\mathfrak{L}$ (e.g.,
$\mathfrak{L}(\phi)=\phi$, $\mathfrak{L}(\phi)=\nabla\phi$, or
$\mathfrak{L}(\phi)=H(\phi)$), we define the limits

$\mathfrak{L}(\phi)^{\pm}(x)=\lim_{\epsilon\to
0^{+}}\mathfrak{L}(\phi)(x\mp\epsilon\mathbf{n}^{+}(x)),\quad\forall\,x\in\mathcal{E}_{\rm
p}^{\rm int},$ (15a) and the average and jump operators (we omit dependence on
$x$)
$\left\\{\mathfrak{L}(\phi)\right\\}=\frac{\mathfrak{L}(\phi)^{+}+\mathfrak{L}(\phi)^{-}}{2},\quad\llbracket\mathfrak{L}(\phi)\rrbracket=\left(\mathfrak{L}(\phi)^{+}-\mathfrak{L}(\phi)^{-}\right)\cdot\mathbf{n}^{+},$
(15b) for all $x\in\mathcal{E}_{\rm p}^{\rm int}$. We clarify that if
$\mathfrak{L}(\phi)$ is a matrix-valued quantity the jump operator involves a
matrix-vector product.

We introduce the space $\mathcal{U}$ for the displacement field (cf. section
2):

$\mathcal{U}=\left\\{\varphi\in[\mathcal{X}_{\rm
hf,p}]^{2}:\varphi\cdot\mathbf{n}\big{|}_{\partial\Omega_{\rm p}}=0\right\\}.$
(16)

We observe that $\mathcal{U}$ (16) is contained in the space of Lipschitz
functions in $\Omega_{\rm p}$ but it is not contained in $C^{1}(\Omega_{\rm
p})$.

In view of the definition of the penalty function in (3) in section 4.4, we
define discrete counterparts of the $H^{2}$ norm and seminorm. Since the FE
space $\mathcal{X}_{\rm hf,p}$ is not $H^{2}$-conforming, the discrete
approximation of the $H^{2}$ inner product, should involve a direct control of
the jump of the mapping derivatives across elements. We hence equip
$\mathcal{U}$ with the inner product

$\begin{array}[]{rl}\displaystyle{(w,v)=}&\displaystyle{\sum_{k=1}^{N_{\rm
e}}\int_{\texttt{D}_{k,\rm p}}\left(H(w):H(v)\,+\,w\cdot
v\right)\,dx}\\\\[8.53581pt] +&\displaystyle{\sum_{j\in\texttt{I}_{\rm
int}}\int_{\texttt{F}_{j,\rm p}}\left(\beta_{j}\llbracket\nabla
w\rrbracket\cdot\llbracket\nabla
v\rrbracket+\frac{1}{\beta_{j}}\\{H(w)\\}:\\{H(v)\\}\right)\,ds,}\\\
\end{array}$ (17)

where $H(w):H(v)=\sum_{i,j,k=1}^{2}\partial_{i,j}w_{k}\partial_{i,j}v_{k}$,
and $\beta_{j}=\sigma_{\beta}\kappa^{2}|\texttt{F}_{j,\rm p}|^{-1}$ with
$\sigma_{\beta}=10$ and $j\in\texttt{I}_{\rm int}$; we further define the
induced norm $\|\cdot\|=\sqrt{(\cdot,\cdot)}$. The choice of the inner product
is inspired by the analysis in [22, section 4] for the biharmonic equation:
the authors of [22] consider the norm

$\begin{array}[]{rl}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|w\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}\,=&\displaystyle{\Big{(}\sum_{k=1}^{N_{\rm
e}}\|\Delta w\|_{L^{2}(\texttt{D}_{k,\rm
p})}^{2}\,+\,\sum_{j\in\texttt{I}_{\rm int}}\|\sqrt{\alpha}\llbracket
w\rrbracket\|_{L^{2}(\texttt{F}_{j,\rm p})}^{2}}\\\\[8.53581pt]
+&\displaystyle{\|\sqrt{\beta}\llbracket\nabla
w\rrbracket\|_{L^{2}(\texttt{F}_{j,\rm p})}^{2}+\|1/\sqrt{\beta}\\{\Delta
w\\}\|_{L^{2}(\texttt{F}_{j,\rm p})}^{2}+\|1/\sqrt{\alpha}\\{\nabla\Delta
w\\}\|_{L^{2}(\texttt{F}_{j,\rm p})}^{2}\Big{)},}\\\ \end{array}$

for proper choices of $\alpha$ and $\beta$. Since we here consider continuous
discretizations, we omit in (17) the first and fourth terms in the facet
integral and we consider the same expression for $\beta$ as in [22];
furthermore, since we wish to consider the full $H^{2}$ norm, we replace
$\Delta w$ with $H(w)$ in the volumetric and surface integral.

We also define the seminorm

$\begin{array}[]{rl}\mathfrak{P}\left(\varphi\right)=&\displaystyle{\sum_{k=1}^{N_{\rm
e}}\int_{\texttt{D}_{k,\rm
p}}\left(H(\varphi):H(\varphi)\right)\,dx}\\\\[8.53581pt]
\,+&\displaystyle{\sum_{j\in\texttt{I}_{\rm int}}\int_{\texttt{F}_{j,\rm
p}}\left(\beta_{j}\llbracket\nabla\varphi\rrbracket\cdot\llbracket\nabla\varphi\rrbracket+\frac{1}{\beta_{j}}\\{H(\varphi)\\}:\\{H(\varphi)\\}\right)\,ds.}\\\
\end{array}$ (18)

Note that $\mathfrak{P}\left(\Phi\right)=0$ for any linear map
$\Phi=\mathbf{b}+\mathbf{A}x$; in particular,
$\mathfrak{P}\left(\texttt{id}\right)=0$. We also define the penalty

$\mathfrak{P}_{\rm brkn}\left(\varphi\right)=\sum_{k=1}^{N_{\rm
e}}\int_{\texttt{D}_{k,\rm p}}\left(H(\varphi):H(\varphi)\right)\,dx,$ (19)

which does not penalize discontinuities of the displacement derivative at
elements’ interfaces.

###  4.3 Construction of the curved mesh of $\Omega$

The polytope $\Omega_{\rm p}$ and the map $\Psi$ are explicitly linked to the
high-order mesh $\mathcal{T}$: the problem of determining $\Omega_{\rm p}$ and
the map $\Psi$ can be recast as the problem of generating a (coarse) high-
order mesh $\mathcal{T}$ of the domain of interest. We propose to first
prescribe the polytope $\Omega_{\rm p}$, then determine a linear triangulation
$\mathcal{T}_{\rm p}$ of $\Omega_{\rm p}$ and finally compute the map $\Psi$
by solving a suitable optimization problem. We introduce the parameterizations
$\\{\gamma_{k}:(0,1)\to\texttt{F}_{k}\\}_{k=1}^{N_{\rm b}}$ of the curved
edges $\\{\texttt{F}_{k}\\}_{k=1}^{N_{\rm b}}$ and the Gauss-Lobatto points
$\\{t_{i}^{\rm gl}\\}_{i=1}^{\kappa+1}\subset[0,1]$. Then, we introduce the
reference and deformed points

$x_{j_{i,k}}=\gamma_{k}(0)(1-t_{i}^{\rm gl})+\gamma_{k}(1)\,t_{i}^{\rm
gl},\quad y_{j_{i,k}}=\gamma_{k}(t_{i}^{\rm gl}),\quad{\rm
where}\;\;j_{i,k}=i+(k-1)N_{\rm b},$

for $i=1,\ldots,\kappa+1$, and $k=1,\ldots,N_{\rm b}$, and the affine space
$\widetilde{\mathcal{W}}$ with

$\widetilde{\mathcal{W}}=\left\\{\texttt{id}+\varphi\in[\mathcal{X}_{\rm
hf,p}]^{2}:\varphi\cdot\mathbf{n}\big{|}_{\partial\Omega_{\rm
p}\cap\partial\Omega}=0\right\\};$ (20)

the set $\partial\Omega_{\rm p}\cap\partial\Omega$ corresponds to the portion
of $\partial\Omega$ that is the union of linear facets and that is hence
shared with $\partial\Omega_{\rm p}$ (see Hypothesis 1). Note that the affine
space $\widetilde{\mathcal{W}}$ enables deformation of curved facets.

In conclusion, we consider the constrained optimization statement:

$\min_{\Psi\in\widetilde{\mathcal{W}}}\;\;\mathfrak{f}_{\rm
jac}(\Psi)\,+\,\mathfrak{P}_{\rm brkn}\left(\Psi\right)\quad{\rm
s.t}\;\;\|\Psi(x_{j})-y_{j}\|_{\infty}\leq\delta,\;\;j=1,\ldots,N=(\kappa+1)N_{\rm
b},$ (21)

where $\mathfrak{f}_{\rm jac}$ is introduced in (7) and is designed to enforce
the bijectivity of the mapping, while $\mathfrak{P}_{\rm brkn}$ is the
quadratic penalty term (19). The penalty (19) was found superior to (18) —
that is, it led to more regular meshes and to more rapid convergence of the
optimizer — for the model problem of section 5.

The optimization statement (21) reads as a nonlinear non-convex optimization
problem with linear inequality constraints; in our numerical implementation,
we resort to the Matlab function fmincon to solve (21) based on an interior-
point method. In the numerical experiments, we consider $\delta=10^{-6}$.

###  4.4 Penalty function

We denote by $\mathcal{T}_{\rm pb}$ the mesh used for HF calculations:
$\mathcal{T}_{\rm pb}$ is independent of the FE mesh $\mathcal{T}_{\rm p}$
used in sections 4.1 and 4.2 to define the mapping space; it might also be
associated to a different (spectral element, finite volume) discretization
method. For model reduction tasks, the registration procedure should ensure
that the deformed mesh $\Phi(\mathcal{T}_{\rm pb})$ is well-behaved (see,
e.g., [3]). The penalty $\mathfrak{f}_{\rm pen}^{\rm th}$ (8) is not tailored
to the discrete representation of the FE fields $\\{\varphi_{i}\\}_{i=1}^{M}$
that are not $H^{2}$-conforming; furthermore, the penalty does not ensure that
the deformed mesh $\Phi(\mathcal{T}_{\rm pb})$ is well-behaved. Below, we
hence derive an alternative penalty function that is well-suited for model
reduction tasks. We denote by $\\{\texttt{D}_{k}^{\rm pb}\\}_{k=1}^{N_{\rm
e}^{\rm pb}}$ and $\\{\Psi_{k}^{\rm pb}\\}_{k=1}^{N_{\rm e}^{\rm pb}}$ the
elements and the elemental mappings of $\mathcal{T}_{\rm pb}$; we use notation
$\\{\Psi_{\Phi,k}^{\rm pb}\\}_{k=1}^{N_{\rm e}^{\rm pb}}$ to refer to the
elemental mappings of the deformed mesh (see, e.g., [34, Eq. (9)]).

We hence propose the penalty:

$\mathfrak{f}_{\rm pen}(\mathbf{a})=\mathfrak{f}_{\rm jac}(\texttt{N}_{\rm
p}(\mathbf{a}))+\mathfrak{f}_{\rm
msh}(\texttt{N}(\mathbf{a}))+\mathfrak{P}\left(\texttt{N}_{\rm
p}(\mathbf{a})\right)$ (22a) where $\mathfrak{P}$ (c.f. (18)) promotes
smoothness, $\mathfrak{f}_{\rm jac}$ (cf. (7)) enforces local bijectivity, and
the function $\mathfrak{f}_{\rm msh}$ controls the quality of the deformed
mesh, $\mathfrak{f}_{\rm msh}({\Phi})=\frac{1}{|\Omega|}\sum_{k=1}^{N_{\rm
e}^{\rm pb}}\;\;|\texttt{D}_{k}^{\rm pb}|\int_{\widehat{\texttt{D}}}{\rm
exp}\left(\frac{q_{k}^{\rm msh}(\Phi)}{q_{k}^{\rm
msh}(\texttt{id})}\,-\,\kappa_{\rm msh}\right)\;dx,$ (22b) where $q_{k}^{\rm
msh}(\Phi):=\frac{1}{d^{2}}\left(\frac{\|\nabla\Psi_{\Phi,k}^{\rm pb}\|_{\rm
F}^{2}}{({\rm det}(\nabla\Psi_{\Phi,k}^{\rm pb}))_{+}^{2/d}}\right)^{2},\quad
d=2.$ (22c) We observe that the ratio $q_{k}^{\rm msh}$ measures the degree of
anisotropy of the mesh and it can thus be interpreted as a measure of the
quality of the deformed mesh — we recall that the indicator $q_{k}^{\rm msh}$
was used in [41] in the framework of DG methods. As discussed in [33], the
decision to activate the penalty terms $\mathfrak{f}_{\rm jac}$ and/or
$\mathfrak{f}_{\rm msh}$ depends on the particular way we treat parameterized
geometries — map-then-discretize or discretize-then-map: in the numerical
experiments, we activate only $\mathfrak{f}_{\rm msh}$ with $\kappa_{\rm
msh}=10$.

In (22), we employ an $L^{2}$ penalization of the Hessian rather than an
$L^{\infty}$ penalization as in (8): this choice simplifies computations and
is justified by the fact that the ansatz $\texttt{N}_{\rm p}$ consists of a
finite expansion of modes.

###  4.5 Target function

The target $\mathfrak{f}_{\mu}^{\rm tg}$ measures the degree of similarity
between the available estimate of the solution field $u_{\mu}$ to the problem
of interest and a suitable template solution or template reduced space;
$\mathfrak{f}_{\mu}^{\rm tg}$ relies on the introduction of a sensor
$s_{\mu}:=\mathfrak{s}_{\mu}(u_{\mu})$ which should highlight the coherent
structures we wish to track. We might distinguish between _point-set sensors_
and _distributed sensors_. The choice of the sensor might reflect the
connection between distributed sensors and the shock sensors that are used to
activate artificial viscosities in high-order discretization [24, 25]. We
anticipate that in the numerical experiments we combine point-set and
distributed sensors to take into account different sources of information and
to robustify the greedy procedure outlined in section 4.6 (cf. (33)).

#### Target function based on point-set sensors

Point-set sensors are based on the introduction of a scalar testing function
that selects the points of the mesh
$\\{x_{\mu,j}^{\star}\\}_{j=1}^{N_{\mu}^{\star}}$ where the solution $u_{\mu}$
satisfies a suitable user-defined criterion (cf. [12]), for any
$\mu\in\mathcal{P}$. Given the template point cloud
$\\{\bar{x}_{j}^{\star}\\}_{j=1}^{\overline{N}^{\star}}$ — which can be
prescribed _a priori_ or be chosen based on one specific value of the
parameters — and the parameter value $\mu\in\mathcal{P}$, we first rely on a
standard point-set registration (PSR) procedure that takes as inputs the point
clouds $\\{\bar{x}_{j}^{\star}\\}_{j=1}^{\overline{N}^{\star}}$ and
$\\{x_{\mu,j}^{\star}\\}_{j=1}^{N_{\mu}^{\star}}$, to determine the deformed
point cloud $\\{\widehat{x}_{j,\mu}^{\star}\\}_{j=1}^{\overline{N}^{\star}}$
that approximates (in a suitable sense) the target cloud
$\\{x_{\mu,j}^{\star}\\}_{j=1}^{N_{\mu}^{\star}}$ (see, e.g., [23] for further
details); then, we solve the optimization problem (3) with

$\mathfrak{f}_{\mu}^{\rm
tg}\left(\Phi\right)=\frac{1}{\overline{N}^{\star}}\sum_{i=1}^{\overline{N}^{\star}}\|\Phi(\bar{x}_{i}^{\star})-\widehat{x}_{i,\mu}^{\star}\|_{2}^{2}.$
(23)

As discussed in [12] and also [33], we can interpret the solution to (3) with
objective (22)-(23) as an approximate projection of the mapping returned by
the PSR procedure — which is neither guaranteed to map the boundary of the
domain $\Omega$ in itself nor to be globally bijective — onto the space of
admissible bijections in $\Omega$.

In order to speed up calculations, it is important to exploit the structure of
the mapping in (4) and precompute $\Psi^{-1}$ where needed before calling the
optimizer. In more detail, in order to speed up the evaluation of
$\mathfrak{f}_{\mu}^{\rm tg}$ in (23), we first compute and store
$\bar{z}_{j}^{\star}=\Psi^{-1}(\bar{x}_{j}^{\star})$ for
$j=1,\ldots,\overline{N}^{\star}$ and then we evaluate (23) as

$\mathfrak{f}_{\mu}^{\rm
tg}\left(\Phi\right)=\frac{1}{\overline{N}^{\star}}\sum_{i=1}^{\overline{N}^{\star}}\|\Psi\left(\Phi_{\rm
p}(\bar{z}_{i}^{\star})\right)-\widehat{x}_{i,\mu}^{\star}\|_{2}^{2}.$

Similar reasoning can be exploited for the evaluation of (22b).

#### Target function based on distributed sensors

A distributed sensor is a function of the solution field $u_{\mu}$ that
highlights the features of $u_{\mu}$ that we wish to track. Given the reduced
$n$-dimensional space $\mathcal{Z}_{n}$ embedded in a Hilbert space
$\mathcal{X}$ defined over $\Omega$, recalling (13), we find that the
parametric mapping $\Phi$ should minimize the target

$\mathfrak{f}_{\mu}^{\rm
tg,opt}\left(\Phi\right):=\min_{\zeta\in\mathcal{Z}_{n}}\int_{\Omega}\|u_{\mu}\circ\Phi-\zeta\|_{2}^{2}\,dx,$
(24)

over all values of $\mu$ in $\mathcal{P}$: $\mathfrak{f}_{\mu}^{\rm
tg,opt}\left(\Phi\right)$ is the projection error in the mapped configuration.
Note that $\mathfrak{f}_{\mu}^{\rm tg,opt}$ depends on the choice of the
reduced space $\mathcal{Z}_{n}$ whose selection is inherently coupled with the
problem of finding $\Phi$: we postpone the procedure for the construction of
the reduced space to the next section.

The explicit use of the solution $u_{\mu}$ in the optimization statement (3)
is computationally unfeasible and prone to instabilities: first, the solution
$u_{\mu}$ is typically defined over an unstructured grid for which function
evaluation at the deformed quadrature points $\\{\Phi(x_{q}^{\rm
qd})\\}_{q=1}^{N_{\rm qd}}$ is extremely expensive; second, we might exploit
prior knowledge about the problem of interest to identify a scalar function of
$u_{\mu}$ that better isolates the features we wish to track (e.g., shocks)
using registration, from the features we expect to be able to approximate
through a linear expansion of (mapped) snapshots.

Exploiting the form of $\Phi$, and the change-of-variable $x=\Psi(\xi)$, we
find

$\int_{\Omega}\|u_{\mu}\circ\Phi-\zeta\|_{2}^{2}\,dx=\int_{\Omega_{\rm
p}}\|u_{\mu}\circ\Psi\circ\Phi_{\rm
p}-\zeta\circ\Psi^{-1}\|_{2}^{2}\,J(\Psi)d\xi.$

If we replace $u_{\mu}\circ\Psi$ with a scalar sensor $s_{\mu}$ defined over
the domain $\Omega_{\rm p}$ and the reduced space $\mathcal{Z}_{n}$ for the
solution with a space (dubbed _template space_) $\mathcal{S}_{n}\subset
L^{2}(\Omega_{\rm p})$ for the sensor, we obtain

$\mathfrak{f}_{\mu}^{\rm
tg}\left(\Phi\right):=\min_{\nu\in\mathcal{S}_{n}}\int_{\Omega_{\rm
p}}\big{|}s_{\mu}\circ\Phi_{\rm p}-\nu\big{|}^{2}J(\Psi)\,dx,$ (25)

which is the target function employed in the numerical experiments.

From the definition (25), we deduce that computation of
$\mathfrak{f}_{\mu}^{\rm tg}(\Phi)$ requires to evaluate $s_{\mu}$ in the
deformed quadrature points $\\{\Phi_{\rm p}(x_{q}^{\rm qd,p})\\}_{q}$ of the
mesh $\mathcal{T}_{\rm p}$ at each iteration of the optimization algorithm for
(3). Since the solution is discontinuous, we consider a P1 discretization of
the sensor $s_{\mu}$ over a linear mesh $\mathcal{T}_{\rm p,s}$ of
$\Omega_{\rm p}$, which is generated independently of $\mathcal{T}_{\rm p}$.
Following [17], we rely on KD-trees to speed up mesh interpolation (cf. Matlab
function KDTreeSearcher): since $\mathcal{T}_{\rm p}$ is a linear mesh with a
modest number of elements, the evaluation of the sensor $s_{\mu}$ in the
deformed quadrature points is still affordable. Notice that the same fast mesh
interpolation routine should also be used to evaluate the FE mapping $\Psi$ in
(22b).

###  4.6 Construction of the parametric map

The target function (25) depends on the template space $\mathcal{S}_{n}\subset
L^{2}(\Omega_{\rm p})$. Following [34, 35], given a set of sensor snapshots
$\\{s_{\mu}:\mu\in\mathcal{P}_{\rm train}\\}$ with $\mathcal{P}_{\rm
train}=\\{\mu^{k}\\}_{k=1}^{n_{\rm train}}$, we propose an iterative procedure
that performs registration over the entire training set and then exploits the
results to update the template space $\mathcal{S}_{n}$ in a greedy fashion.
Algorithm 1 summarizes the computational procedure.

Algorithm 1 Registration algorithm

_Inputs:_ $\\{s_{\mu}:\mu\in\mathcal{P}_{\rm train}\\}$ snapshot set,
$\mathcal{S}_{n_{0}}={\rm span}\\{s_{\mu^{\star,(i)}}\\}_{i=1}^{n_{0}}$
initial template space; $\mathcal{T}_{\rm pb}$ mesh for HF computations.

_Outputs:_ ${\mathcal{S}}_{n}$ template space,
$\mathbf{W}\in\mathbf{R}^{M\times m}$ mapping space,
$\\{\mathbf{a}_{\mu^{k}}^{\star}\\}_{k}\subset\mathbb{R}^{m}$ optimal mapping
coefficients.

1:Initialization: $\mathcal{S}_{n=n_{0}}=\mathcal{S}_{n_{0}}$,
$\Xi_{\star}=\\{\mu^{\star,(i)}\\}_{i=1}^{n_{0}}$,
$\mathbf{W}=\mathbbm{1}_{M}$.

2:for $n=n_{0},\ldots,n_{\rm max}-1$ do

3:
$\left[\mathbf{a}_{\mu}^{\star},\mathfrak{f}_{\mu}^{\star}\right]\,=\,\texttt{registration}\left(s_{\mu},\,\mathcal{S}_{n},\,\mathbf{W},\,\mathcal{T}_{\rm
pb},\,\Psi,\mathbf{a}_{\mu}^{0}\right)$ for all $\mu\in\mathcal{P}_{\rm
train}$ _see Remark 4.1 for definition of $\mathbf{a}_{\mu}^{0}$_

4: $[\mathbf{W},\;\\{\mathbf{a}_{\mu}^{\rm proj}\\}_{\mu\in\mathcal{P}_{\rm
train}}]=\texttt{POD}\left(\\{\mathbf{a}_{\mu}^{\star}\\}_{\mu\in\mathcal{P}_{\rm
train}},tol_{\rm pod},(\cdot,\cdot)_{2}\right),$

5: if $\max_{\mu\in\mathcal{P}_{\rm
train}}\mathfrak{f}_{\mu}^{\star}<\texttt{tol}$ then

6: break

7: else

8: $\Xi_{\star}=\Xi_{\star}\cup\\{\mu^{\star,(n+1)}\\}$ with
$\mu^{\star,(n+1)}={\rm arg}\max_{\mu\in\mathcal{P}_{\rm
train}}\mathfrak{f}_{\mu}^{\star}$.

9: $\mathcal{S}_{n+1}={\rm span}\\{s_{\mu^{\star,(i)}}\circ\texttt{N}_{\rm
p}(\mathbf{a}_{\mu^{\star,(i)}}^{\rm proj})\\}_{i=1}^{n+1}$.

10: end if

11:end for

Given the orthonormal basis $\\{\varphi_{i}\\}_{i}$ of $\mathcal{U}$ and the
operator N, we define the orthogonal matrix $\mathbf{W}\in\mathbb{R}^{M\times
m}$ and the low-rank map
$\widehat{\texttt{N}}(\mathbf{a}):=\texttt{N}(\mathbf{W}\mathbf{a})$; by
construction, we have
$\|\mathbf{a}\|_{2}=\|\sum_{i}(\mathbf{W}\mathbf{a})\varphi_{i}\|$ for any
$\mathbf{a}\in\mathbb{R}^{m}$. Then, we introduce notation:

$\left[\mathbf{a}_{\mu}^{\star},\mathfrak{f}_{\mu}^{\star}\right]\,=\,\texttt{registration}\left(s_{\mu},\,\mathcal{S}_{n},\,\mathbf{W},\,\mathcal{T}_{\rm
pb},\,\Psi,\,\mathbf{a}_{\mu}^{0}\right)$

to refer to the function that takes as inputs (i) the target sensor
$s_{\mu}:\Omega_{\rm p}\to\mathbb{R}$, (ii) the template space
$\mathcal{S}_{n}$, (iii) the orthogonal matrix $\mathbf{W}$, (iv) the mesh
$\mathcal{T}_{\rm pb}$, (v) the geometric mapping $\Psi:\Omega_{\rm
p}\to\Omega$ and (vi) the initial guess
$\mathbf{a}_{\mu}^{0}\in\mathbb{R}^{m}$ for the optimizer, and returns (I) the
mapping coefficients $\mathbf{a}_{\mu}^{\star}$ that minimize the (reduced)
objective $\mathbf{a}\mapsto\mathfrak{f}_{\mu}^{\rm
obj}(\mathbf{W}\mathbf{a}_{\mu}^{\star})$, and (II) the value of the target
function $\mathfrak{f}_{\mu}^{\star}=\mathfrak{f}_{\mu}^{\rm
tg}(\mathbf{W}\mathbf{a})$. Note that the objective $\mathfrak{f}_{\mu}^{\rm
obj}:\mathbb{R}^{M}\to\mathbb{R}_{+}$ of (3) depends on the mesh
$\mathcal{T}_{\rm pb}$ through the term $\mathfrak{f}_{\rm msh}$ in (22b). We
also introduce the function

$[\mathbf{W}_{\rm new},\;\\{\mathbf{a}_{\mu}^{\rm
proj}\\}_{\mu\in\mathcal{P}_{\rm train}}]=\texttt{POD}\left(\\{\mathbf{W}_{\rm
old}\mathbf{a}_{\mu}^{\star}\\}_{\mu\in\mathcal{P}_{\rm train}},tol_{\rm
pod},(\cdot,\cdot)_{2}\right),$

which implements POD based on the method of snapshots with Euclidean inner
product $(\cdot,\cdot)_{2}$: the tolerance $tol_{\rm pod}>0$ drives the
selection of the number of modes $m$ based on the energy criterion

$m:=\min\left\\{m^{\prime}:\,\sum_{j=1}^{m^{\prime}}\lambda_{j}\geq\left(1-tol_{\rm
pod}\right)\sum_{i=1}^{n_{\rm train}}\lambda_{i}\right\\},$ (26)

where $\lambda_{1}\geq\ldots\geq\lambda_{n_{\rm train}}\geq 0$ are the
eigenvalues of the Gramian matrix $\mathbf{C}\in\mathbb{R}^{n_{\rm
train}\times n_{\rm train}}$ such that
$(\mathbf{C})_{k,k^{\prime}}=\mathbf{a}_{\mu^{k}}^{\star}\cdot\mathbf{a}_{\mu^{k^{\prime}}}^{\star}$.
The function POD returns also the mapping coefficients associated with the
projected displacements $\mathbf{a}_{\mu}^{\rm proj}=\mathbf{W}_{\rm
new}^{\top}\mathbf{W}_{\rm old}\mathbf{a}_{\mu}^{\star}$; the latter are used
to initialize the iterative method for the optimization problem for the
subsequent iterations.

######  Remark 4.1.

Further implementation details. Since the optimization problem is highly non-
convex, the choice of the initial condition is extremely important to avoid
convergence to unsatisfactory local minima. For
$n=n_{0}+1,n_{0}+2,\ldots,n_{\rm max}-1$, we simply use
$\mathbf{a}_{\mu}^{0}=\mathbf{a}_{\mu}^{\rm proj}$ (cf. Line 4). On the other
hand, for the first iteration, we first reorder the parameters in
$\mathcal{P}_{\rm train}$ so that $\mu^{(1)}={\rm
arg}\min_{\mu\in\mathcal{P}_{\rm train}}\|\mu-\mu^{\star,(1)}\|_{2}$ and

$\mu^{(k)}={\rm arg}\min_{\mu\in\mathcal{P}_{\rm
train}\setminus\\{\mu^{(i)}\\}_{i=1}^{k-1}}\left(\min_{\mu^{\prime}\in\\{\mu^{(i)}\\}_{i=1}^{k-1}}\|\mu-\mu^{\prime}\|_{2}\right),\quad
k=2,\ldots,n_{\rm train};$

then, we choose the initial condition as follows:

$\mathbf{a}_{\mu^{(1)}}^{0}=0,\quad\mathbf{a}_{\mu^{(k)}}^{0}={\mathbf{a}}_{\mu^{({\rm
ne}_{k})}}^{\star}\;\;{\rm with}\;\;{\rm ne}_{k}={\rm
arg}\,\min_{j=1,\ldots,k-1}\|\mu^{(j)}-\mu^{(k)}\|_{2},$

for $k=2,\ldots,n_{\rm train}$. In a previous implementation of the procedure,
we also included box constraints in the optimization statement (cf. [32,
section 3.1.2]) to control the sensitivity of the mapping coefficients to
parameter variations;

$\|\mathbf{a}_{\mu^{(k)}}^{\star}-\mathbf{a}_{\mu^{(k)}}^{0}\|_{\infty}\leq
C_{\infty}\|\mu^{(k)}-\mu^{({\rm ne}_{k})}\|_{2},\quad{\rm
with}\;\;\;C_{\infty}=10;$ (27)

however, in the numerical experiments of the present work, we found that the
solution to the unconstrained problem satisfied the constraints (27) for all
the experiments considered. We further observe that Algorithm 1 depends on
several hyper-parameters. In our tests, we set $\mathcal{S}_{n_{0}=1}={\rm
span}\\{s_{\bar{\mu}}\\}$, where $\bar{\mu}$ is the centroid of $\mathcal{P}$;
furthermore, we set $n_{\rm max}=6$, $tol_{\rm pod}=5\cdot 10^{-3}$ and
$\texttt{tol}=10^{-4}$.

######  Remark 4.2.

Generalization. Given the dataset
$\\{(\mu^{k},\mathbf{a}_{\mu^{k}}^{\star})\\}_{k=1}^{n_{\rm train}}$ as
provided by Algorithm 1, we resort to a multi-target regression algorithm to
learn a regressor $\mu\mapsto\widehat{\mathbf{a}}_{\mu}$ for the mapping
coefficients, and ultimately define the parametric mapping

${\Phi}:\Omega\times\mathcal{P}\to\Omega,\quad{\Phi}_{\mu}:=\texttt{N}\left(\mathbf{W}\widehat{\mathbf{a}}_{\mu}\right).$
(28)

We here resort to radial basis function (RBF, [38]) approximation: other
regression algorithms could also be considered. Similarly to [32, 35], to
avoid overfitting, we verify the statistical significance of the RBF
estimators. We randomly split the dataset
$\\{(\mu^{k},\mathbf{a}_{\mu^{k}}^{\star})\\}_{k=1}^{n_{\rm train}}$ into the
learning and test sets $\\{(\mu^{k},\mathbf{a}_{\mu^{k}}^{\star})\\}_{k\in
D_{\rm learn}}$ and $\\{(\mu^{j},\mathbf{a}_{\mu^{j}}^{\star})\\}_{j\in D_{\rm
test}}$ with $D_{\rm learn}\cap D_{\rm test}=\emptyset$ and $D_{\rm learn}\cup
D_{\rm test}=\\{1,\ldots,n_{\rm train}\\}$ (we here consider a $80\%$-$20\%$
learning/test split); we compute the RBF approximation
$\widehat{\mathbf{a}}:\mathcal{P}\to\mathbb{R}^{m}$ based on the learning set
and we compute the out-of-sample R-squared coefficient for each component:

$\texttt{R}_{i}^{2}=1-\frac{\sum_{j\in D_{\rm
test}}\left(\mathbf{a}_{\mu^{j}}^{\star}-\widehat{\mathbf{a}}_{\mu^{j}}\right)_{i}^{2}}{\sum_{j\in
D_{\rm
test}}\left(\mathbf{a}_{\mu^{j}}^{\star}-\bar{\mathbf{a}}\right)_{i}^{2}},\quad\bar{\mathbf{a}}=\frac{1}{|D_{\rm
learn}|}\sum_{k\in D_{\rm learn}}\,\mathbf{a}_{\mu^{k}}^{\star},\quad
i=1,\ldots,m.$ (29)

Then, we retain exclusively modes for which $\texttt{R}_{i}^{2}$ is above the
threshold $R_{\rm min}=0.70$.

######  Remark 4.3.

Parametric registration based on point-set sensors. The greedy procedure in
Algorithm 1 is motivated by the need to construct the template space
$\mathcal{S}_{n}$. If we rely on point-set sensors, we do not have to perform
multiple iterations. However, we empirically found that performing two
iterations of the for loop in Algorithm 1 does not hinder computational
efficiency — since the cost is dominated by the high-dimensional registration
problems solved during the first iteration — and has a beneficial effect on
generalization outside the training set.

##  5 Application to an inviscid flow past an array of LS89 turbine blades

We consider the problem of estimating the solution to the two-dimensional
Euler equations past an array of LS89 turbine blades.

###  5.1 Model problem

We consider the computational domain depicted in Figure 9(a); we prescribe
total temperature, total pressure and flow direction at the inflow, static
pressure at the outflow, non-penetration condition on the blade and periodic
boundary conditions on the lower and upper boundaries. We study the
sensitivity of the solution with respect to two parameters: the free-stream
Mach number ${\rm Ma}_{\infty}$ and the height of the channel $H$, i.e.
$\mu=[H,{\rm Ma}_{\infty}]$. We consider the parameter domain
$\mathcal{P}=[0.95,1.05]\times[0.9,0.95]$. We refer to [3] for a detailed
presentation of the employed nondimensionalization, HF DG formulation and
pseudo-transient continuation strategy.

Figure 7 shows the distribution of the Mach field for $\mu_{\rm
min}=[0.95,0.9]$ and $\mu_{\rm max}=[1.05,0.95]$; Figure 8(a) shows the
behavior of the Mach number on the upper side of the blade for four parameter
values, while Figure 8(b) shows the behavior of the entropy profile ${\rm
E}=\log(p)-\gamma\log(\rho)$ where $p$ is the pressure field, $\gamma=1.4$ is
the ratio of specific heats, and $\rho$ is the density field. The solution
develops a normal shock on the upper side of the blade for sufficiently large
values of ${\rm Ma}_{\infty}$ and $H$; furthermore, the entropy ${\rm E}$
exhibits several peaks that correspond to the blade wakes. The solution
develops two shocks at the trailing edge, which are highly undesirable for
turbomachinery applications: we expect that at the trailing edge viscous
effects might not be negligible; for this reason, a more thorough
investigation should rely on a model that accounts for viscous effects. We
observe that the shock location and the entropy peaks are sensitive to the
value of the parameter: this justifies the application of registration
procedures.

(a)

(b)

Figure 7: inviscid flow past an array of LS89 turbine blades. (a)-(b) Mach
field for $\mu=[0.95,0.9]$ and $\mu=[1.05,0.95]$.

(a)

(b)

Figure 8: inviscid flow past an array of LS89 turbine blades. (a) Mach profile
over the upper side of the blade for four values of the parameter in
$\mathcal{P}$. (b) entropy profile on the bottom boundary for four values of
the parameter in $\mathcal{P}$. The parameter $s$ denotes the normalized
curvilinear coordinate, $s\in[0,1]$.

###  5.2 Definition of the mapping ansatz

We deal with geometry variations through a piecewise-smooth mapping associated
with the partition in Figure 9(b). We set $H_{\rm ref}=1$ and we define the
curve $x_{1}\mapsto f_{\rm btm}(x_{1})$ that describes the lower boundary
$\Gamma_{\rm btm}$ of the domain $\Omega=\Omega(H=H_{\mathrm{ref}})$; then, we
define $\widetilde{H}>0$ such that $x_{1}\mapsto f_{\rm
btm}(x_{1})+\widetilde{H}$ and $x_{1}\mapsto f_{\rm
btm}(x_{1})+H-\widetilde{H}$ do not intersect the blade for any
$H\in[0.95,1.05]$; finally, we define the geometric mapping

$\Psi_{H}^{\rm geo}(x=[x_{1},x_{2}])=\left[\begin{array}[]{l}x_{1}\\\
\psi_{H}^{\rm geo}(x)\\\ \end{array}\right],\quad x\in\Omega,$ (30a) where
$\psi_{H}^{\rm
geo}(x)=\left\\{\begin{array}[]{ll}o_{1}(x_{1})+C(H)\left(x_{2}-o_{1}(x_{1})\right)&x_{2}<o_{1}(x_{1}),\\\
o_{2}(x_{1})+C(H)\left(x_{2}-o_{2}(x_{1})\right)&x_{2}>o_{2}(x_{1}),\\\
x_{2}&{\rm otherwise},\\\ \end{array}\right.$ (30b) with $o_{1}(x_{1})=f_{\rm
btm}(x_{1})+\widetilde{H}$, $o_{2}(x_{1})=f_{\rm btm}(x_{1})+H_{\rm
ref}-\widetilde{H}$ and $C(H)=\frac{H-H_{\rm ref}}{2\widetilde{H}}+1$.

Then, we consider computational maps (4) from the domain
$\Omega=\Omega(H=H_{\rm ref})$ to $\Omega_{\mu}=\Omega(H)$ such that

$\texttt{N}(x,\mathbf{a},\mu)=\Psi_{H}^{\rm geo}\circ\Psi\circ\texttt{N}_{\rm
p}(\mathbf{a})\circ\Psi^{-1}(x),$ (31)

for a proper choice of the mapping $\Psi$ and the polytope $\Omega_{\rm p}$.
Note that $\widetilde{\texttt{N}}(\mathbf{a})=\Psi\circ\texttt{N}_{\rm
p}(\mathbf{a})\circ\Psi^{-1}$ defines a bijection in $\Omega$.

(a)

(b)

(c)

(d)

Figure 9: inviscid flow past an array of LS89 turbine blades. (a)
computational mesh $\mathcal{T}_{\rm pb}$ for DG calculations; red dots denote
the nodes on the upper blade side; blue squares denote the nodes on the lower
boundary. (b) partition associated with the geometric map (30). (c)-(d)
coarse-grained linear mesh $\mathcal{T}_{\rm p}$ and curved mesh $\mathcal{T}$
for registration.

The polytope $\Omega_{\rm p}$ and the mapping $\Psi:\Omega_{\rm p}\to\Omega$
should be designed to ensure the periodicity constraint

$\texttt{N}(x+[0,H_{\rm
ref}],\mathbf{a},\mu)=\texttt{N}(x,\mathbf{a},\mu)+[0,H],\quad\forall\,x\in\Gamma_{\rm
btm}.$ (32)

The geometric mapping $\Psi^{\rm geo}$ satisfies the periodicity constraint;
it is thus sufficient to enforce that $\widetilde{\texttt{N}}(x+[0,H_{\rm
ref}],\mathbf{a},\mu)=\widetilde{\texttt{N}}(x,\mathbf{a},\mu)+[0,H_{\rm
ref}]$ for all $x\in\Gamma_{\rm btm}$. Towards this end, (i) we define the
polytope $\Omega_{\rm p}$ with vertices $V$ such that $V\cap\Gamma_{\rm
btm,p}=V\cap\Gamma_{\rm top,p}+[0,H_{\rm ref}]$, where $\Gamma_{\rm
btm,p},\Gamma_{\rm top,p}$ denote the lower and upper boundaries of the
polytope (i.e. for every vertex on the lower boundary there exists a
corresponding vertex on the upper boundary); (ii) we define the linear mesh
$\mathcal{T}_{\rm p}$ in Figure 9(c) with matching nodes on the periodic
boundaries222The same condition is enforced in the HF mesh.; (iii) we define
the mapping $\Psi$ by solving (21) in the affine space
$\widetilde{W}=\left\\{\texttt{id}+\varphi\in[\mathcal{X}_{\rm
hf,p}]^{2}:\varphi\cdot\mathbf{n}\big{|}_{\partial\Omega_{\rm
p}\cap\partial\Omega}=0,\;\;\varphi|_{\Gamma_{\rm
btm,p}}=\varphi|_{\Gamma_{\rm top,p}}\right\\}$ and we set ${\mathcal{U}}$
(cf. (16)) as ${\mathcal{U}}=\left\\{\varphi\in[\mathcal{X}_{\rm
hf,p}]^{2}:\varphi\cdot\mathbf{n}\big{|}_{\partial\Omega_{\rm
p}}=0,\;\;\varphi|_{\Gamma_{\rm btm,p}}=\varphi|_{\Gamma_{\rm
top,p}}\right\\}$. The resulting curved mesh is provided in Figure 9(d); in
the experiments, we consider polynomials of degree ten, which implies
$M=1359$. Note that, since $\widetilde{\mathcal{W}}$ and ${\mathcal{U}}$ are
spaces of piecewise polynomials, it suffices to enforce the periodicity
constraint at mesh nodes.

###  5.3 Definition of the sensor

We consider a target function $\mathfrak{f}^{\rm tg}$ that combines a point-
set sensor and a distributed sensor

$\mathfrak{f}_{\mu}^{\rm
tg}=\frac{1}{4}\sum_{i=1}^{4}\|\texttt{N}(\bar{x}_{i},\mathbf{a},\mu)-\widehat{x}_{i,\mu}\|_{2}^{2}\,+\,\min_{\nu\in\mathcal{S}_{n}}\int_{\Omega_{\rm
p}}\big{|}s_{\mu}\circ\texttt{N}_{\rm
p}(\mathbf{a})-\nu\big{|}^{2}J(\Psi)\,dx,$ (33)

where the distributed sensor $s_{\mu}$ is obtained from the Mach field, the
template space $\mathcal{S}_{n}$ is built adaptively using Algorithm 1, and
the salient points $\\{\widehat{x}_{i,\mu}\\}_{i}$ correspond to the first
three peaks of the entropy profile on the lower boundary (cf. Figure 8(b)) and
the shock location on the upper boundary (cf. Figure 8(a)).

We identify the entropy peaks by computing the local maxima of the entropy
profile on $\Gamma_{\rm btm,\mu}$. We detect the shock location on the upper
blade side using the following procedure: first, we compute the mean values
$\\{\overline{{\rm Ma}}_{j}\\}_{j}$ and $\\{\overline{d{\rm Ma}}_{j}\\}_{j}$
of the Mach number and its tangential derivative on each facet
$\\{\texttt{F}_{j}\\}_{j}$ of the select boundary; second, we find the index
$j^{\star}$ such that $\overline{{\rm Ma}}_{j^{\star}}>1>\overline{{\rm
Ma}}_{j^{\star}+1}$ and $\overline{d{\rm
Ma}}_{j^{\star}}<-\frac{10^{-2}}{|\texttt{F}_{j^{\star}}|}$ and we return the
estimate ${x}_{\rm shk,\mu}^{\rm raw}$ equal to the midpoint of the selected
facet.

Given the raw estimates $\\{{x}_{i,\mu}^{\rm
raw}:i=1,\ldots,4,\mu\in\mathcal{P}_{\rm train}\\}$ of the point sensors, we
apply radial basis function (RBF) regression to find smoother estimates that
facilitate the generalization step (cf. Remark 4.2); in order to ensure that
the points are on the boundary, we apply RBF to the curvilinear coordinates of
the points. As shown in Figure 8(a), the solution does not exhibit any shock
for several parameter values; we do not consider these points to train the RBF
surrogate, and we rely on the surrogate itself to find a fictitious estimate
of $\widehat{x}_{\rm shk,\mu}$ for all $\mu\in\mathcal{P}$.

###  5.4 Performance of the registration procedure

We apply Algorithm 1 based on a regular $11\times 6$ grid of parameters
$\mathcal{P}_{\rm train}$. Figure 10 compares the location of the sensors
$\\{\widetilde{x}_{i,\mu}=\left(\Psi_{H}^{\rm
geo}\right)^{-1}(\widehat{x}_{i,\mu}):i=1,\ldots,4,\mu\in\mathcal{P}_{\rm
train}\\}$ and their mapped counterparts
$\\{\widetilde{x}_{i,\mu}^{\prime}=\Phi_{\mu}^{-1}(\widehat{x}_{i,\mu}):i=1,\ldots,4,\mu\in\mathcal{P}_{\rm
train}\\}$; similarly, Figure 11 reproduces the results in Figure 8 for the
mapped solution field. We observe that the peaks of the mapped entropy field
on $\Gamma_{\rm btm}$ and the shock on the upper side of the blade are nearly
insensitive to the parameter value. For all in-sample and out-of-sample
configurations considered, the quality of the deformed mesh — which is
measured by (22c) — is comparable with the one of the original mesh.

A thorough analysis of the computational cost of the procedure is beyond the
scope of the present paper. We observe that the execution time of Algorithm 1
is $54$ minutes on a commodity Linux workstation; for comparison, the
generation of the training set of HF simulations ($66$ simulations) requires
$7$ hours and $33$ minutes on the same machine.

(a)

(b)

Figure 10: inviscid flow past an array of LS89 turbine blades. (a) sensor
points used in the registration procedure for $n_{\rm train}=66$
configurations. (b) mapped sensor points.

(a)

(b)

Figure 11: inviscid flow past an array of LS89 turbine blades. (a) mapped Mach
profile over the upper side of the blade for four values of the parameter in
$\mathcal{P}$. (b) mapped entropy profile on the bottom boundary for four
values of the parameter in $\mathcal{P}$. Unmapped profiles are provided in
Figure 8.

###  5.5 Model reduction

We rely on POD+regression to estimate the solution field. POD solution
coefficients are estimated using the same RBF regression method that is
employed for the mapping coefficients (cf. Remark 4.2). We denote by
$u_{\mu}^{\rm hf}$ the HF solution, and by $\widehat{u}_{\mu}$ the POD+RBF
estimate. To assess performance, we consider a dataset of $n_{\rm test}=20$
randomly-selected parameters $\mathcal{P}_{\rm test}$ and we compute the
$L^{2}$ maximum relative errors

$E_{\rm max}=\max_{\mu\in\mathcal{P}_{\rm
test}}\frac{\|\widehat{u}_{\mu}-{u}_{\mu}^{\rm
hf}\|_{L^{2}(\Omega_{\mu})}}{\|{u}_{\mu}^{\rm
hf}\|_{L^{2}(\Omega_{\mu})}},\quad E_{\rm max}^{\rm
bnd}=\max_{\mu\in\mathcal{P}_{\rm
test}}\frac{\|\widehat{u}_{\mu}-{u}_{\mu}^{\rm hf}\|_{L^{2}(\Gamma_{\rm
bld})}}{\|{u}_{\mu}^{\rm hf}\|_{L^{2}(\Gamma_{\rm bld})}},$ (34)

where $\Gamma_{\rm bld}$ denotes the boundary of the blade. Figure 12 shows
the results for the “linear” ROM and the “registered” ROM — both approaches
involve the application of a geometric mapping; the difference is that for the
linear ROM the mapping is chosen a priori and corresponds to $\Psi^{\rm geo}$
(30), while for the registered ROM the mapping is chosen through the
registration procedure discussed in the paper. Both linear and registered ROMs
reach a plateau for $n\approx 10$ due to the limited amount of datapoints. We
observe that registration improves performance by a factor $2.04$ for the
global error and a factor $8.04$ for the error on the blade, for the same
amount of HF data.

$5$$10$$15$$20$$10^{-4}$$10^{-3}$$10^{-2}$$n$$E_{\rm max}$ (a)

$5$$10$$15$$20$$10^{-4}$$10^{-3}$$10^{-2}$$n$$E_{\rm max}^{\rm bnd}$ (b)

Figure 12:  inviscid flow past an array of LS89 turbine blades; model
reduction. Behavior of the out-of-sample errors (34) for the linear ROM

12a

and the registered ROM

12a

for $n_{\rm train}=66$ and several values of $n$.

##  6 Summary and discussion

We addressed the problem of parametric registration for manifolds associated
to the solution to parametric PDEs: registration is designed to track coherent
structures of the solution field, and ultimately simplify the task of MOR. We
proposed and analyzed a new class of compositional maps (4) for registration
in two-dimensional domains. We provided sufficient conditions to ensure the
bijectivity of the mapping $\Phi$ (cf. Proposition 2.1 and Corollary 2.3) and
we studied the approximation properties: we found that the ansatz (4) is dense
in a meaningful subspace of diffeomorphisms only for polytopal domains
(Corollary 2.4, Lemma 3.2); on the other hand, the multi-layer generalization
(5) is dense for arbitrary domains under the assumption of small deformation
(cf. Lemma 3.3). We discussed an actionable procedure to determine the ansatz
(4) for arbitrary domains, which exploits a coarse-grained high-order FE mesh;
finally, we illustrated the performance of the method for a parametric
compressible flow past a cascade of turbine blades.

We aim to extend our approach in several directions. First, we aim to extend
Proposition 2.1 to (a class of) Lipschitz maps. Second, we wish to automate
the construction of the polytope $\Omega_{\rm p}$; we also wish to devise
effective implementations of the ansatz (5). Third, we wish to extend our
method to three-dimensional domains: this extension likely requires the
development of specialized optimization routines to cope with high-dimensional
mapping spaces (i.e., large $M$), and efficient mesh interpolation routines to
speed up the evaluation of the objective function in (3).

As extensively discussed in section 3.4, our choice of the ansatz (4) is a
compromise between approximation power and simplicity of implementation. In
the future, we aim to investigate generalizations of (4) to improve the
approximation power of our model class; in this respect, we plan to develop
registration techniques based on the ansatz (5).

Finally, we remark that registration procedures constitute a key elements of
MOR techniques based on Lagrangian approaches. The results of this paper
contribute to motivate the interest in this class of methods for problems with
coherent localized structures. Nevertheless, their inadequacy to cope with
topology changes (e.g., two shocks that merge into one shock [35, Appendix
E.1]) might justify the use of different nonlinear proposals (see, e.g., [2,
4]). It is hence of paramount importance to perform thorough comparative
studies to study the strengths and the weaknesses of Lagrangian methods and to
compare performance with other nonlinear approaches.

## Acknowledgments

The author acknowledges the support by European Union’s Horizon 2020 research
and innovation programme under the Marie Sklodowska-Curie Actions, grant
agreement 872442 (ARIA), and the support provided by Inria under the programme
Exploratory Actions (project AM2OR). The author also thanks Prof. Angelo Iollo
and Dr. Pierre Mounoud (Univ. Bordeaux), Dr. Alberto Remigi (Safran Tech), and
Prof. Masayuki Yano (University of Toronto) for fruitful discussions on
various aspects of the method. The author also thanks three anonymous
reviewers for their many comments that help improve the quality of the
manuscript.

##  Appendix A Proof of Proposition 2.1

The key challenge of the proof is that the domain $\Omega_{\rm p}$ is not in
general simply connected: it is hence not possible to directly employ the
result in [32] based on Hadamard’s global inverse function theorem. To address
this issue, (i) we first prove a variant of Proposition 2.1 in [32] (cf.
Proposition A.1) that addresses the case of simply connected domains (cf.
Proposition A.1); (ii) we prove that the restriction of $\Phi$ to
$\partial\Omega_{\rm p}$ is a diffeomorphism (cf. Lemma A.1); (iii) we prove a
technical result to glue together $C^{1}$ maps with positive Jacobian
determinants (cf. Lemma A.2); (iv) we exploit the previous results and a well-
known extension theorem in differential geometry to prove the desired result.
Below, we present a number of preliminary results that address (i), (ii) and
(iii) and then we discuss the proof of Proposition 2.1.

######  Proposition A.1.

Let $\Omega_{\rm p}\subset\mathbb{R}^{2}$ be a polytope isomorphic to the unit
ball. Consider the vector-valued function ${\Phi}:\Omega_{\rm
p}\to\mathbb{R}^{2}$ that satisfies (i) $\Phi\in
C^{1}\left(\overline{\Omega}_{\rm p};\mathbb{R}^{2}\right)$; (ii)
$|J({\Phi})(x)|=\big{|}{\rm
det}\left(\nabla\Phi(x)\right)\big{|}\geq\epsilon>0$ for all $x\in\Omega_{\rm
p}$ and a given $\epsilon>0$; (iii) ${\Phi}(\partial\Omega_{\rm
p})\subseteq\partial\Omega_{\rm p}$. Then, ${\Phi}$ is a bijection that maps
$\Omega_{\rm p}$ into $\Omega_{\rm p}$.

###### Proof.

This result is a variant of Proposition 2.1 in [32] for polytopes $\Omega_{\rm
p}$ isomorphic to the unit ball; the proof follows the exact same steps of the
one in [32]; details are omitted. In [32] we assumed that $\Phi$ is of class
$C^{1}$ in a $\delta$-neighborhood of $\Omega_{\rm p}$; since $\Omega_{\rm p}$
is a Lipschitz domain, if $\Phi$ is of class $C^{1}$ up to the boundary, then
there exists a $C^{1}$ extension of $\Phi$ to $\mathbb{R}^{2}$ (cf. [6,
section 2.5]): we can thus replace the hypothesis (i) of [32, Proposition 2.1]
with the hypothesis ${\Phi}\in C^{1}\left(\overline{\Omega}_{\rm
p};\mathbb{R}^{2}\right)$. Furthermore, we replace
$J({\Phi})(x)\geq\epsilon>0$ for all $x\in\Omega_{\rm p}$ with the more
general assumption $|J({\Phi})(x)|\geq\epsilon>0$: since $J({\Phi})$ is
continuous, this condition means that $J({\Phi})$ is either strictly positive
or strictly negative in $\Omega_{\rm p}$. ∎

######  Lemma A.1.

Let $\Phi=\texttt{id}+\varphi$ be a $C^{1}(\overline{\Omega}_{\rm
p};\mathbb{R}^{2})$ field such that $\inf_{x\in\Omega_{\rm p}}J(\Phi)>0$ and
$\varphi\cdot\mathbf{n}|_{\partial\Omega_{\rm p}}=0$. Then,
$\Phi(\partial\Omega_{\rm p})=\partial\Omega_{\rm p}$, and
$\Phi|_{\partial\Omega_{\rm p}}$ is a diffeomorphism in $\partial\Omega_{\rm
p}$.

###### Proof.

It suffices to prove that each edge $F$ is mapped in itself and that the
restriction of $\Phi$ to $F$ is a bijection. Towards this end, we define the
vertices $x_{1}^{\rm v},x_{2}^{\rm v}$ of $F$, the tangent and normal vectors
$\mathbf{t}_{\rm f},\mathbf{n}_{\rm f}$ and the parameterization $\gamma_{\rm
f}:[0,1]\to F$ such that $\gamma_{\rm f}(s)=x_{1}^{\rm v}+s\|x_{2}^{\rm
v}-x_{1}^{\rm v}\|_{2}\mathbf{t}_{\rm f}$. We also introduce the normal
vectors $\mathbf{n}_{\rm f}^{-},\mathbf{n}_{\rm f}^{+}$ of the neighboring
edges that meet $F$ at $x_{1}^{\rm v},x_{2}^{\rm v}$, respectively.

We first observe that $\varphi(x_{1}^{\rm v})=\varphi(x_{2}^{\rm v})=0$: since
$F$ meets with the neighboring edges, the pairs of normals $\mathbf{n}_{\rm
f}^{-},\mathbf{n}_{\rm f}$ and $\mathbf{n}_{\rm f}^{+},\mathbf{n}_{\rm f}$
should be linearly independent (i.e., they should span $\mathbb{R}^{2}$);
then, exploiting the continuity of $\varphi$ up to the boundary of
$\Omega_{\rm p}$ — and in particular on $\partial\Omega_{\rm p}$ — we find
$\varphi(x_{1}^{\rm v})\cdot\mathbf{n}_{\rm f}^{-}=\varphi(x_{1}^{\rm
v})\cdot\mathbf{n}_{\rm f}=0$ and $\varphi(x_{2}^{\rm v})\cdot\mathbf{n}_{\rm
f}=\varphi(x_{2}^{\rm v})\cdot\mathbf{n}_{\rm f}^{+}=0$, which imply
$\varphi(x_{1}^{\rm v})=\varphi(x_{2}^{\rm v})=0$.

By contradiction, $\Phi(F)\not\subset F$. Given $x\in F$, we find

$\varphi(x)=(\varphi(x)\cdot\mathbf{n}_{\rm f})\mathbf{n}_{\rm
f}+(\varphi(x)\cdot\mathbf{t}_{\rm f})\mathbf{t}_{\rm
f}=(\varphi(x)\cdot\mathbf{t}_{\rm f})\mathbf{t}_{\rm f}$

and thus $\Phi(\gamma_{\rm f}(s))=x_{1}^{\rm v}+\alpha(s)\mathbf{t}_{\rm f}$
with $\alpha(s)=(\Phi(\gamma_{\rm f}(s))-x_{1}^{\rm v})\cdot\mathbf{t}_{\rm
f}.$ Since $\varphi(x_{1}^{\rm v})=\varphi(x_{2}^{\rm v})=0$, we find
$\alpha(0)=0,\alpha(1)=\|x_{2}^{\rm v}-x_{1}^{\rm v}\|_{2}$: exploiting the
intermediate value theorem, we find that $[0,1]\subset\alpha([0,1])$ and thus
$F\subset\Phi(F)$. The condition $\Phi(F)\not\subset F$ implies that
$\alpha(s)\notin[0,1]$ for some $s\in[0,1]$ and thus
$\alpha^{\prime}(s^{\star})=\lVert
x_{2}^{\mathrm{v}}-x_{1}^{\mathrm{v}}\rVert_{2}\mathbf{t}_{\rm
f}^{T}\nabla\Phi(\gamma_{\rm f}(s^{\star}))\mathbf{t}_{\rm f}=0$ for some
$s^{\star}\in[0,1]$. Note that $\Phi(\gamma_{\rm f}(s))\cdot\mathbf{n}_{\rm
f}=x_{1}^{\rm v}\cdot\mathbf{n}_{\rm f}$ for any $s\in[0,1]$; if we
differentiate left- and right-hand sides of the previous identity with respect
to $s$ we obtain $\lVert
x_{2}^{\mathrm{v}}-x_{1}^{\mathrm{v}}\rVert_{2}\mathbf{n}_{\rm
f}^{\top}\nabla\Phi(\gamma_{\rm f}(s))\mathbf{t}_{\rm f}=0$ for all
$s\in[0,1]$. Therefore, if we set $x^{\star}=\gamma_{\rm f}(s^{\star})$, we
find $\nabla\Phi(x^{\star})\mathbf{t}_{\rm f}=\left(\mathbf{t}_{\rm
f}^{\top}\nabla\Phi(x^{\star})\mathbf{t}_{\rm f}\right)\mathbf{t}_{\rm
f}+\left(\mathbf{n}_{\rm f}^{T}\nabla\Phi(x^{\star})\mathbf{t}_{\rm
f}\right)\mathbf{n}_{\rm f}=0,$ which implies that $\nabla\Phi$ is singular at
$x^{\star}$ and thus $\min_{x\in\overline{\Omega}_{\rm p}}J(\Phi)\leq
J(\Phi)(x^{\star})=0$. Contradiction.

So far, we proved that the $C^{1}$ function $\alpha:[0,1]\to\mathbb{R}$ is
surjective in $[0,1]$ and is monotonic increasing. Therefore, $\alpha$ is
bijective in $[0,1]$ and thus $\Phi|_{F}$ is a bijection in $F$. ∎

######  Lemma A.2.

Let $\Omega=\Omega_{\rm ext}\setminus\bigcup_{i=1}^{N}\Omega_{{\rm int},i}$ be
a polytope that satisfies the conditions of Definition 2.1. We define
$\Omega_{\rm int}=\bigcup_{i=1}^{N}\Omega_{{\rm int},i}$. Let
$\Phi^{+}:\Omega\to\mathbb{R}^{2}$ and $\Phi^{-}:\Omega_{\rm
int}\to\mathbb{R}^{2}$ be $C^{1}$ maps up to the boundary that coincide on
$\Gamma:=\partial\Omega_{\rm int}$ and satisfy $J(\Phi^{+})(x)\geq\epsilon>0$
for all $x\in\Omega$, $J(\Phi^{-})(x)\geq\epsilon>0$ for all $x\in\Omega_{\rm
int}$. We define $\Phi:\Omega_{\rm ext}\to\mathbb{R}^{2}$ that is equal to
$\Phi^{+}$ in $\Omega$ and is equal to $\Phi^{-}$ in $\Omega_{\rm int}$. Then,
there exists a sequence $\\{\Phi_{j}\\}_{j}\subset
C^{1}(\overline{\Omega};\mathbb{R}^{2})$ such that (i) $\Phi_{j}=\Phi$ on
$\partial\Omega_{\rm ext}$, (ii)
$\lim_{j\to\infty}\|\Phi_{j}-\Phi\|_{L^{\infty}(\Omega_{\rm ext})}=0$, (iii)
$\min_{x\in\overline{\Omega}}J(\Phi_{j})\geq\frac{1}{2}\epsilon$ for
$j=1,2,\ldots$.

###### Proof.

We assume that $\Omega_{\rm int}$ is isomorphic to the unit ball (that is,
$\Omega_{{\rm int}}=\Omega_{{\rm int},1}$ and $N=1$); the proof can be
extended to $N>1$ through an inductive process, iteratively incorporating one
subdomain at each step. We define the piecewise-linear parameterization
$\boldsymbol{\gamma}:(0,1)\to\Gamma$; we denote by $\mathbf{t}$ the tangent
vector to $\Gamma$, by $\mathbf{n}$ the normal vector that points towards
$\Omega$, and by $V_{\rm int}$ the vertices of $\Omega$ on $\Gamma$; ${\rm
dist}_{\rm H}(A,B)$ is the Hausdorff distance between the sets $A$ and $B$.
Given $\delta>0$, we define the set $V_{\rm int,\delta}=\bigcup_{x\in V_{\rm
int}}\mathcal{B}_{\delta}(x)$, and the $C^{1}$ approximation $\Gamma_{\delta}$
of $\Gamma$ such that ${\rm dist}_{\rm H}(\Gamma_{\delta},\Gamma)\leq\delta$
and ${\rm dist}_{\rm H}(\Gamma_{\delta}\setminus V_{\rm
int,\delta},\Gamma\setminus V_{\rm int,\delta})=0$. Then, we introduce the
$C^{1}$ parameterization
$\boldsymbol{\gamma}_{\delta}:(0,1)\to\Gamma_{\delta}$, the corresponding unit
normal $\mathbf{n}_{\delta}$ and the set
$\Gamma_{\eta,\delta}=\left\\{\boldsymbol{\gamma}_{\delta}(u)+\eta
v\mathbf{n}_{\delta}(\boldsymbol{\gamma}_{\delta}(u))\,:\,v\in\left(-f(u),1-f(u)\right),\;\;u\in(0,1)\right\\}$,
where $f:(0,1)\to[0,1]$ is a suitable smooth function that will be introduced
below. Figure 13 provides a graphic explanation of the quantities introduced
so far.

Exploiting the previous definitions, we define
$\psi_{\eta,\delta}:\mathbb{R}^{2}\to[0,1]$ such that (i)
$\psi_{\eta,\delta}(\boldsymbol{\gamma}_{\delta}(u)+\eta
v\mathbf{n}_{\delta})=\phi(v+f(u))$ for all
$u\in(0,1),v\in\left(-f(u),1-f(u)\right)$, with $\phi(t)=-2t^{3}+3t^{2}$, (ii)
$\psi_{\eta,\delta}|_{\Omega_{\rm ext}\setminus\Gamma_{\eta,\delta}}=1$, and
(iii) $\psi_{\eta,\delta}|_{\Omega_{\rm int}\setminus\Gamma_{\eta,\delta}}=0$.
It is straightforward to verify that $\psi_{\eta,\delta}$ is of class $C^{1}$,
and that

$\nabla\psi_{\eta,\delta}(x)=\left\\{\begin{array}[]{ll}\displaystyle{\frac{1}{\eta}\left(\alpha_{\psi,\delta}(x)\mathbf{n}+o(1)\right)}&\displaystyle{\forall\,x\in\Gamma_{\eta,\delta}\setminus
V_{\rm int,\delta},}\\\\[8.53581pt]
\displaystyle{\mathcal{O}\left(\frac{\alpha_{\psi,\delta}(x)}{\eta}\mathbf{n}_{\delta}+\frac{\beta_{\psi,\delta}(x)}{\delta}\mathbf{t}_{\delta}\right)}&\displaystyle{\forall\,x\in
V_{\rm int,\delta},}\\\ \end{array}\right.$

where
$\alpha_{\psi,\delta}(x)\in(0,\max_{v\in\mathbb{R}}\phi^{\prime}(v))=(0,3/2)$
and $\beta_{\psi,\delta}(x)$ belongs to a bounded interval that is independent
of $\delta>0$. In the previous expression, we used the fact that
$\textbf{n}_{\delta}=\textbf{n}$ for $x\in\Gamma_{\eta,\delta}\setminus V_{\rm
int,\delta}$, since ${\rm dist}_{\rm H}(\Gamma_{\delta}\setminus V_{\rm
int,\delta},\Gamma\setminus V_{\rm int,\delta})=0$.

We define the positive vanishing sequences
$\\{\eta_{j}\\}_{j},\\{\delta_{j}\\}_{j}\subset\mathbb{R}_{+}$; for
$j=1,2,\ldots$, we introduce the function
$\psi_{j}=\psi_{\eta_{j},\delta_{j}}$ and the field $\Phi_{j}=\left(\Phi_{\rm
e}^{+}-\Phi_{\rm e}^{-}\right)\psi_{j}+\Phi_{\rm e}^{-}$, where $\Phi_{\rm
e}^{+},\Phi_{\rm e}^{-}$ are $C^{1}$ extensions of $\Phi^{+},\Phi^{-}$ to
$\mathbb{R}^{2}$ (cf. [6, section 2.5]). Clearly, $\\{\Phi_{j}\\}_{j}$
satisfies the conditions (i) and (ii); furthermore, $J(\Phi_{j})\geq\epsilon$
in $\Omega_{\rm ext}\setminus\Gamma_{\eta_{j},\delta_{j}}$; in the remainder,
we prove that $J(\Phi_{j})\geq\epsilon/2$ in $\Gamma_{\eta_{j},\delta_{j}}$.

Let $x\in\Gamma_{\eta_{j},\delta_{j}}\setminus V_{\rm int,\delta_{j}}$. Since
$\Phi^{+}=\Phi^{-}$ on $\Gamma$, we have $\left(\nabla(\Phi_{\rm
e}^{+}-\Phi_{\rm e}^{-})\right)\mathbf{t}=0$ for all $x\in\Gamma$. We define
$\mathbf{t}_{\Phi}=\frac{\nabla\Phi_{\rm e}^{+}\mathbf{t}}{\|\nabla\Phi_{\rm
e}^{+}\mathbf{t}\|_{2}}$ and
$\mathbf{n}_{\Phi}=[(\mathbf{t}_{\Phi})_{2},-(\mathbf{t}_{\Phi})_{1}]$. Given
$x\in\Gamma_{\eta_{j},\delta_{j}}$,
$x=\boldsymbol{\gamma}(u)+\eta_{j}v\mathbf{n}$, we define
$\bar{x}=\boldsymbol{\gamma}(u)$ and $\bar{y}=\Phi_{\rm
e}^{+}(\bar{x})=\Phi_{\rm e}^{-}(\bar{x})$. Then,

$\begin{array}[]{rl}\Phi_{\rm
e}^{\pm}(x)\approx&\displaystyle{\bar{y}+\nabla\Phi_{\rm
e}^{\pm}(\bar{x})(x-\bar{x})}\\\\[8.53581pt]
=&\displaystyle{\bar{y}+\nabla\Phi_{\rm
e}^{\pm}(\bar{x})\mathbf{t}(\bar{x})\mathbf{t}^{\top}(\bar{x})(x-\bar{x})+\nabla\Phi_{\rm
e}^{\pm}(\bar{x})\mathbf{n}(\bar{x})\mathbf{n}^{\top}(\bar{x})(x-\bar{x})}\\\\[8.53581pt]
=&\displaystyle{\bar{y}+\left(\alpha(\bar{x})\mathbf{t}_{\Phi}(\bar{x})\mathbf{t}^{\top}(\bar{x})+\beta^{\pm}(\bar{x})\mathbf{t}_{\Phi}(\bar{x})\mathbf{n}^{\top}(\bar{x})+\gamma^{\pm}\mathbf{n}_{\Phi}(\bar{x})\mathbf{n}^{\top}(\bar{x})\right)(x-\bar{x})},\\\
\end{array}$

which implies $\nabla\Phi_{\rm
e}^{\pm}=\alpha\mathbf{t}_{\Phi}\mathbf{t}^{T}+\beta^{\pm}\mathbf{t}_{\Phi}\mathbf{n}^{T}+\gamma^{\pm}\mathbf{n}_{\Phi}\mathbf{n}^{T}$,
where $\alpha,\beta^{\pm},\gamma^{\pm},\mathbf{t}_{\Phi}$ depend on $\bar{x}$.
Exploiting the latter, we find

$\begin{array}[]{l}\displaystyle{\nabla\Phi_{j}(x)=\left(\nabla\Phi_{\rm
e}^{+}(x)-\nabla\Phi_{\rm e}^{-}(x)\right)\psi_{j}(x)+\left(\Phi_{\rm
e}^{+}(x)-\Phi_{\rm e}^{-}(x)\right)(\nabla\psi_{j}(x))^{\top}+\nabla\Phi_{\rm
e}^{-}(x)}\\\\[8.53581pt]
\approx\displaystyle{\alpha(\bar{x})\mathbf{t}_{\Phi}(\bar{x})\mathbf{t}^{\top}(\bar{x})+\left(\beta^{-}(\bar{x})+\llbracket\beta(\bar{x})\rrbracket\psi_{j}\right)\mathbf{t}_{\Phi}(\bar{x})\mathbf{n}^{\top}(\bar{x})+}\\\\[8.53581pt]
\displaystyle{\left(\gamma^{-}(\bar{x})+\llbracket\gamma(\bar{x})\rrbracket\psi_{j}(\bar{x})\right)\mathbf{n}_{\Phi}(\bar{x})\mathbf{n}^{\top}(\bar{x})+v\alpha_{\psi}(\bar{x})\left(\llbracket\beta(\bar{x})\rrbracket\mathbf{t}_{\Phi}(\bar{x})+\llbracket\gamma(\bar{x})\rrbracket\mathbf{n}_{\Phi}(\bar{x})\right)\mathbf{n}^{\top}(\bar{x})},\\\
\end{array}$

for all $x\in\Gamma_{\eta}$. Since $\\{\mathbf{t},\mathbf{n}\\}$ and
$\\{\mathbf{t}_{\Phi},\mathbf{n}_{\Phi}\\}$ are orthonormal bases of
$\mathbb{R}^{2}$, we must have

$J(\Phi_{j})=\alpha\left(\gamma^{-}+\llbracket\gamma\rrbracket\left(\psi_{j}+v\alpha_{\psi}\right)\right)+o(1).$

Let $\alpha$ be positive; the case $\alpha<0$ is analogous. If
$\llbracket\gamma\rrbracket>0$, we find

$\gamma^{-}+\llbracket\gamma\rrbracket\Big{(}\underbrace{\psi_{\eta}}_{\geq
0}+\underbrace{v}_{\geq-f(u)}\underbrace{\alpha_{\psi}}_{\leq
3/2}\big{)}\geq\gamma-\frac{3}{2}\llbracket\gamma\rrbracket f(u),$

which exceeds $\gamma^{-}/2$ if
$f(u)\leq\frac{\gamma^{-}}{3\llbracket\gamma\rrbracket}$. If
$\llbracket\gamma\rrbracket<0$, we find that
$J(\Phi_{\eta})\geq\alpha\frac{\gamma_{+}}{2}$ provided that $f(u)\geq
1-\frac{\gamma^{+}}{3|\llbracket\gamma\rrbracket|}$. Note that for
$\llbracket\gamma\rrbracket=0$ the bounds reduce to $-\infty\leq
f(u)\leq\infty$: we can hence find a smooth function $f$ such that
$J(\Phi_{j})\geq\frac{\alpha}{2}\min\\{\gamma^{+},\gamma^{-}\\}\geq\frac{\epsilon}{2}$.

Let $x\in\Gamma_{\eta_{j},\delta_{j}}\cap V_{\rm int,\delta_{j}}$; we assume
that $\eta_{j}=\mathcal{O}(\delta_{j})$ for $j=1,2,\ldots$. Exploiting the
same argument as before, we find $\nabla\Phi_{\rm e}^{+}=\nabla\Phi_{\rm
e}^{-}$ for $\bar{x}\in V_{\rm int}$, which implies that $\Phi_{\rm
e}^{+}(x)=\Phi_{\rm e}^{-}(x)+o(\delta)$ for any $x\in V_{\rm int}$.
Therefore, by considering a Taylor expansion centered in $x\in V_{\rm int}$,
we find

$\begin{array}[]{l}\nabla\Phi_{j}(x)=\displaystyle{\left(\nabla\Phi_{\rm
e}^{+}(x)-\nabla\Phi_{\rm e}^{-}(x)\right)\psi_{j}(x)+\left(\Phi_{\rm
e}^{+}(x)-\Phi_{\rm e}^{-}(x)\right)(\nabla\psi_{j}(x))^{\top}+\nabla\Phi_{\rm
e}^{-}(x)}\\\\[8.53581pt]
=\displaystyle{o(1)\mathcal{O}(1)+o(\delta)\mathcal{O}\left(\frac{\alpha_{\psi,\delta}(x)}{\eta}\mathbf{n}_{\delta}(x)+\frac{\beta_{\psi,\delta}(x)}{\delta}\mathbf{t}_{\delta}(x)\right)+\nabla\Phi_{\rm
e}^{-}(\bar{x})\overset{\rm(i)}{=}\nabla\Phi_{\rm e}^{-}(\bar{x})+o(1)}\\\
\end{array}$

which implies $J(\Phi_{j})(x)\gtrsim\epsilon$. In (i), we exploited the
hypothesis $\eta_{j}=\mathcal{O}(\delta_{j})$ for $j=1,2,\ldots$. ∎

(a)

(b)

Figure 13: proof of Lemma A.2. (a) example of regular polytope. (b)
visualization of $\Gamma$, $\Gamma_{\delta}$, $\Gamma_{\delta,\eta}$, $V_{\rm
int,\eta}$ in the proximity of one vertex.

###### Proof.

(_Proposition 2.1_). Exploiting Lemma A.1, we find that
$\Phi|_{\partial\Omega_{\rm int}}$ is a diffeomorphism in $\partial\Omega_{\rm
int}$. Therefore, exploiting Cerf’s theorem (see, e.g., [16]), we find that,
for $i=1,\ldots,N$, there exists a diffeomorphism $\Psi_{i}:\Omega_{{\rm
int},i}\to\Omega_{{\rm int},i}$ such that $\Psi_{i}=\Phi$ on
$\partial\Omega_{{\rm int},i}$ and $|J(\Psi_{i})|$ does not vanish in
$\overline{\Omega}_{{\rm int},i}$. If we define $\Omega_{\rm
int}=\bigcup_{i=1}^{N}\Omega_{{\rm int},i}$ and $\Psi:\Omega_{\rm
int}\to\Omega_{\rm int}$ such that $\Psi|_{\partial\Omega_{{\rm
int},i}}=\Psi_{i}$, we hence find that $\Psi$ is a diffeomorphism in
$\Omega_{\rm int}$. Exploiting the same argument as in Lemma A.2, we find that
$\nabla\Psi=\nabla\Phi$ on the vertices of $\partial\Omega_{\rm int}$;
therefore, $J(\Psi)$ is strictly positive in $\overline{\Omega}_{\rm int}$. We
set $\epsilon=\\{\inf_{x\in\Omega_{\rm p}}J(\Phi),\inf_{x\in\Omega_{\rm
int}}J(\Psi)\\}$. We now define the extended map $\Phi_{\rm ext}:\Omega_{\rm
ext}\to\mathbb{R}^{2}$ that is equal to $\Phi$ in $\Omega_{\rm p}$ and is
equal to $\Psi$ in $\Omega_{\rm int}$. Below, we prove that $\Phi_{\rm ext}$
is a bijection (injective and surjective) in $\Omega_{\rm ext}$: since $\Psi$
is bijective in $\overline{\Omega}_{\rm int}$, this implies that $\Phi$ is
bijective in $\Omega_{\rm p}$. Towards this end, we introduce a sequence
$\\{\Phi_{j}\\}_{j}\subset C^{1}(\overline{\Omega}_{\rm ext};\mathbb{R}^{2})$
such that (i) $\Phi_{j}=\Phi$ on $\partial\Omega_{\rm ext}$, (ii)
$\lim_{j\to\infty}\delta_{j}=0$ with $\delta_{j}=\|\Phi_{j}-\Phi_{\rm
ext}\|_{L^{\infty}(\Omega_{\rm ext})}=0$, (iii)
$\min_{x\in\overline{\Omega}_{\rm ext}}J(\Phi_{j})\geq\frac{1}{2}\epsilon$ for
$j=1,2,\ldots$ (cf. Lemma A.2). For any $j>0$, $\Phi_{j}$ satisfies the
hypotheses of Proposition A.1: it is hence bijective in ${\Omega}_{\rm ext}$.

Injectivity of $\Phi_{\rm ext}$. By contradiction, assume that $\Phi_{\rm
ext}$ is not injective, that is, there exist $x,y\in{\Omega}_{\rm ext}$ such
that $\Phi_{\rm ext}(x)=\Phi_{\rm ext}(y)$. We distinguish between the case in
which $x,y\in\partial\Omega_{\rm int}$ and the case in which
$x\notin\partial\Omega_{\rm int}$ (the case $y\notin\partial\Omega_{\rm int}$
is analogous to the latter): since $\Phi|_{\partial\Omega_{\rm int}}$ is a
diffeomorphism in $\partial\Omega_{\rm int}$ (cf. Lemma A.1), the first
scenario is impossible; we hence focus on the second case.

We denote by $B_{x}$ the ball centered in $x$ of radius $\eta>0$. Since
$x\notin\partial\Omega_{\rm int}$, there exists $\eta>0$ such that $y\notin
B_{x}$ and $\partial\Omega_{\rm int}\cap B_{x}=\emptyset$. Therefore,
recalling notation introduced in Lemma A.2, there exists $j_{0}\in\mathbb{N}$
such that $\Gamma_{\eta_{j},\delta_{j}}\cap B_{x}=\emptyset$ for all
$j>j_{0}$: this implies that, for any $j>j_{0}$,
$\Phi_{j}(x^{\prime})=\Phi(x^{\prime})$ for any $x^{\prime}\in B_{x}$;
furthermore, since $\Phi$ is a local diffeomorphism, there exist $c>0$ such
that
$\|\Phi_{j}(x^{\prime})-\Phi_{j}(x)\|_{2}=\|\Phi(x^{\prime})-\Phi(x)\|_{2}\geq
c\eta$ for any $x^{\prime}\in\partial B_{x}$; finally, since $\Phi_{j}$ is a
global diffeomorphism, we also have
$\|\Phi_{j}(x^{\prime})-\Phi_{j}(x)\|_{2}\geq c\eta$ for any $x^{\prime}\notin
B_{x}$.

Exploiting the previous identities and the reverse triangle inequality, we
find

$\begin{array}[]{rl}0=&\|\Phi(y)-\Phi(x)\|_{2}=\|\Phi(y)-\Phi_{j}(x)\|_{2}=\|\Phi(y)-\Phi_{j}(y)+\Phi_{j}(y)-\Phi_{j}(x)\|_{2}\\\\[8.53581pt]
\geq&\|\Phi(y)-\Phi_{j}(y)\|_{2}-\|\Phi_{j}(y)-\Phi_{j}(x)\|_{2}\\\\[8.53581pt]
\Rightarrow&\|\Phi_{j}(y)-\Phi_{j}(x)\|_{2}\leq\|\Phi(y)-\Phi_{j}(y)\|_{2}.\\\
\end{array}$

Since $\|\Phi_{j}(y)-\Phi_{j}(x)\|_{2}\geq c\eta$ and
$\|\Phi(y)-\Phi_{j}(y)\|_{2}\leq\delta_{j}$, we find $c\eta\leq\delta_{j}$ for
any $j>j_{0}$, and thus $c\eta=0$. Contradiction.

Surjectivity of $\Phi_{\rm ext}$. In order to prove that $\Phi_{\rm ext}$ is
surjective in ${\Omega}_{\rm ext}$, that is, $\Phi_{\rm ext}({\Omega}_{\rm
ext})={\Omega}_{\rm ext}$, we first observe that $\Phi_{\rm
ext}(\overline{\Omega}_{\rm ext})\subset\overline{\Omega}_{\rm ext}$. We have
indeed that, for any $y\in\overline{\Omega}_{\rm ext}$,

${\rm dist}\left(y,\Phi_{\rm ext}({\Omega}_{\rm ext})\right)={\rm
dist}\left(\Phi_{j}(x_{j}),\Phi_{\rm ext}({\Omega}_{\rm ext})\right)\leq{\rm
dist}\left(\Phi_{j}(x_{j}),\Phi_{\rm ext}(x_{j})\right)\leq\delta_{j},$

where $x_{j}\in\overline{\Omega}_{\rm ext}$ is the preimage of $y$ through
$\Phi_{j}$; since the latter holds for any $j$, we must have ${\rm
dist}\left(y,\Phi_{\rm ext}(\overline{\Omega}_{\rm ext})\right)=0$. Second, we
observe that $\overline{\Omega}_{\rm ext}\subset\Phi_{\rm
ext}(\overline{\Omega}_{\rm ext})$. We find indeed that for any
$x\in\overline{\Omega}_{\rm ext}$

${\rm dist}\left(\Phi_{\rm ext}(x),\overline{\Omega}_{\rm ext}\right)={\rm
dist}\left(\Phi_{\rm ext}(x),\Phi_{j}(\overline{\Omega}_{\rm
ext})\right)\leq\delta_{j},$

for $j=1,2,\ldots$, which implies ${\rm dist}\left(\Phi_{\rm
ext}(x),\overline{\Omega}_{\rm ext}\right)=0$. In conclusion, we find
$\Phi_{\rm ext}(\overline{\Omega}_{\rm ext})=\overline{\Omega}_{\rm ext}$:
since $\Phi_{\rm ext}(\partial\Omega_{\rm ext})=\partial\Omega_{\rm ext}$, we
must have $\Phi_{\rm ext}(\Omega_{\rm ext})=\Omega_{\rm ext}$. ∎

##  Appendix B Notation

To ease the presentation, we distinguish between basic definitions,
definitions associated with the parametric problem of interest, definitions
associated with the optimization statement, and definitions associated with
the finite element approximation of the compositional map. Below, $\omega$
denotes a generic domain in $\mathbb{R}^{2}$, which is either $\Omega$ or
$\Omega_{\rm p}$ in the main body of the manuscript. To shorten notation, the
explicit dependence on the domain $\omega$ is omitted if not strictly
necessary.

Basic definitions.

  * •

$\Omega\subset\mathbb{R}^{2}$ Lipschitz domain.

  * •

$\texttt{id}:\mathbb{R}^{2}\to\mathbb{R}^{2}$ identity map.

  * •

$\mathfrak{B}(\omega)$ set of Lipschitz bijections from the domain $\omega$ in
itself.

  * •

$\mathfrak{D}(\omega)$ set of diffeomorphisms from the domain $\omega$ in
itself.

  * •

$J(\Phi)={\rm det}(\nabla\Phi)$ Jacobian determinant of the vector-valued
field $\Phi:\omega\to\mathbb{R}^{2}$; $H(\Phi)$ Hessian of $\Phi$.

  * •

$\mathbf{n}_{\omega}:\partial\omega\to\mathbb{S}^{1}=\\{x\in\mathbb{R}^{2}:\|x\|_{2}=1\\}$
outward normal to the domain $\Omega$.

  * •

${\mathfrak{U}}_{0}(\omega)=\\{\varphi\in
C^{1}(\overline{\omega};\mathbb{R}^{2})\,:\,\varphi\cdot\mathbf{n}|_{\partial\omega}=0\\}$.

  * •

$\Omega_{\rm p}\subset\mathbb{R}^{2}$ polytope isomorphic to the domain
$\Omega$.

  * •

$V=\\{x_{i}^{\rm v}\\}_{i=1}^{N_{\rm v}}$ minimal set of vertices of
$\Omega_{\rm p}$.

  * •

$\Psi:\Omega_{\rm p}\to\Omega$ bijection from the polytope $\Omega_{\rm p}$ to
$\Omega$.

  * •

$\\{F_{j}^{\rm p}\\}_{j=1}^{N_{\rm f,p}}$ facets of $\Omega_{\rm p}$;
$\\{F_{j}=\Psi(F_{j}^{\rm p})\\}_{j=1}^{N_{\rm f,p}}$ mapped facets.

  * •

$V_{\rm ang}(\omega)$ angular points of $\partial\omega$ (discontinuities of
the normal $\mathbf{n}_{\omega}$).

  * •

$\overline{\omega}$ closure of $\omega$.

  * •

${\rm dist}_{\rm H}(\omega,\omega^{\prime})$ Hausdorff distance between
$\omega$ and $\omega^{\prime}$.

  * •

${\rm col}(\mathbf{W})$ subspace of $\mathbb{R}^{M}$ spanned by the columns of
the matrix $\mathbf{W}\in\mathbb{R}^{M\times m}$.

Parametric problem.

  * •

$\mu$ vector of parameters in the parameter domain
$\mathcal{P}\subset\mathbb{R}^{P}$.

  * •

$u:\Omega\times\mathcal{P}\to\mathbb{R}^{D_{\rm u}}$ parametric field of
interest, $u_{\mu}(\cdot):=u(\cdot,\mu)$.

  * •

$\mathcal{M}=\\{u_{\mu}:\mu\in\mathcal{P}\\}$ solution manifold.

  * •

$\widetilde{\mathcal{M}}=\\{{u}_{\mu}\circ\Phi_{\mu}^{-1}:\mu\in\mathcal{P}\\}$
mapped solution manifold.

  * •

$\mathcal{T}_{\rm pb}$ high-fidelity mesh with elements
$\\{\texttt{D}_{k}^{\rm pb}\\}_{k=1}^{N_{\rm e}^{\rm pb}}$ and facets
$\\{\Psi_{k}^{\rm pb}\\}_{k=1}^{N_{\rm e}^{\rm pb}}$ employed to approximate
the elements of $\widetilde{\mathcal{M}}$.

Optimization statement (I): general definitions and analysis.

  * •

Minimization statement: (3).

  * •

$\mathfrak{f}^{\rm tg}:\mathfrak{B}(\Omega)\times\mathcal{P}\to\mathbb{R}$
target function.

  * •

$\mathfrak{f}_{\rm pen}:\mathbb{R}^{M}\to\mathbb{R}$ penalty function.

  * •

$\texttt{N}:\mathbb{R}^{M}\to{\rm Lip}(\Omega;\mathbb{R}^{2})$ mapping
operator.

  * •

$A$ (cf. (6)) admissible set for the operator N.

  * •

$A_{\rm jac}$ (cf. (6)) maps with strictly positive Jacobian determinant.

  * •

(4) compositional maps. The pair $(\Omega_{\rm p},\Psi)$ satisfies Hypothesis
1.

  * •

$\mathfrak{f}_{\rm pen}^{\rm th}$ (cf. (8)) penalty function used for the
analysis of section 2

  * •

(5) multi-layer compositional maps. The pairs $(\Omega_{{\rm
p},1},\Psi_{1}),\ldots,(\Omega_{{\rm p},\ell},\Psi_{\ell})$ satisfy Hypothesis
1.

Finite-element discretization of the mapping operator (cf. section 4).

  * •

$\mathbb{P}_{\kappa}(\mathbb{R}^{2})$ space of two-variable polynomials of
total degree less or equal to $\kappa$.

  * •

$\widehat{\texttt{D}}=\\{x\in(0,1)^{2}:\sum_{i=1}^{2}(x)_{i}<1\\}$ reference
(or master) element. $\\{\ell_{i}^{\rm fe}\\}_{i=1}^{n_{\rm lp}}$ Lagrangian
basis of $\mathbb{P}_{\kappa}(\mathbb{R}^{2})$ associated with the nodes
$\\{\tilde{x}_{i}\\}_{i=1}^{n_{\rm
lp}}\subset\overline{\widehat{\texttt{D}}}$; the nodes include the vertices
$\\{\tilde{x}_{i}^{\rm p}\\}_{i=1}^{3}$ of the triangle
$\widehat{\texttt{D}}$.

  * •

$\\{\ell_{i}^{\rm fe,p}\\}_{i=1}^{3}$ Lagrangian basis of $\mathbb{P}_{1}$
associated with the vertices $\\{\tilde{x}_{i}^{\rm p}\\}_{i=1}^{3}$ of
$\widehat{\texttt{D}}$.

  * •

$\mathcal{T}$ curved mesh employed to define the polytope $\Omega_{\rm p}$
(cf. (14b)) and the mapping $\Psi$ (cf. (14d));
$\\{\texttt{D}_{k}\\}_{k=1}^{N_{\rm e}}$ elements of the mesh,
$\\{\texttt{F}_{j}\\}_{j=1}^{N_{\rm f}}$ facets of the mesh, $\\{x_{i,k}^{\rm
hf}\\}_{i,k}$ nodes of the mesh ($x_{i,k}^{\rm hf}$ is the $i$-th node of the
$k$-th element of the mesh).

  * •

$\Psi_{k}$ FE mapping from $\widehat{\texttt{D}}$ to the $k$-th element of the
mesh, $k=1,\ldots,N_{\rm e}$ (cf. (14a)).

  * •

$\Psi_{k,\rm p}$ FE mapping from $\widehat{\texttt{D}}$ to the $k$-th
linearized element of the mesh, $k=1,\ldots,N_{\rm e}$ (cf. (14a)).

  * •

$\mathcal{T}_{\rm p}$ linear mesh with elements $\\{\texttt{D}_{k,\rm
p}\\}_{k=1}^{N_{\rm e}}$ and facets $\\{\texttt{F}_{j,\rm p}\\}_{k=1}^{N_{\rm
f}}$.

  * •

$\mathcal{E}_{\rm p}^{\rm int}=\bigcup_{j\in\texttt{I}_{\rm
int}}\texttt{F}_{j,\rm p}$ internal facets of $\mathcal{T}_{\rm p}$.

  * •

$\mathbf{n}^{+}:\mathcal{E}_{\rm p}^{\rm int}\to\mathbb{S}^{1}$ positive
normal.

  * •

$\mathcal{U}$ (16) FE displacement space equipped with the inner product (17).

Optimization statement (II): practical definitions.

  * •

$\mathfrak{f}_{\rm pen}$ penalty function (22a): $\mathfrak{P}$ (18) smoothing
term; $\mathfrak{f}_{\rm jac}$ (7) bijectivity-enforcing penalty term;
$\mathfrak{f}_{\rm msh}$ (22b) discrete bijectivity-enforcing penalty term;

  * •

(21) optimization problem for the construction of $\Psi$:
$\widetilde{\mathcal{W}}$ (20) search space; $\mathfrak{P}_{\rm brkn}$
smoothing term (19).

  * •

Target function based on point-set sensor (23); target function based on
distributed sensor (25), $\mathcal{S}_{n}\subset L^{2}(\Omega_{\rm p})$
template space for the distributed sensor $s:\mathcal{P}\to L^{2}(\Omega_{\rm
p})$.

## References

  * [1] J. H. Argyris, I. Fried, and D. W. Scharpf.  The TUBA family of plate elements for the matrix displacement method.  The Aeronautical Journal, 72(692):701–709, 1968.
  * [2] J. Barnett, C. Farhat, and Y. Maday.  Neural-network-augmented projection-based model order reduction for mitigating the Kolmogorov barrier to reducibility.  Journal of Computational Physics, 492:112420, 2023.
  * [3] N. Barral, T. Taddei, and I. Tifouti.  Registration-based model reduction of parameterized PDEs with spatio-parameter adaptivity.  Journal of Computational Physics, page 112727, 2023.
  * [4] F. Black, P. Schulze, and B. Unger.  Efficient wildland fire simulation via nonlinear model order reduction.  Fluids, 6(8):280, 2021.
  * [5] T. Blickhan.  A registration method for reduced basis problems using linear optimal transport.  arXiv preprint arXiv:2304.14884, 2023.
  * [6] A. Brudnyi and Y. Brudnyi.  Methods of geometric analysis in extension and trace problems: volume 1, volume 102.  Springer Science & Business Media, 2011.
  * [7] V. Camion and L. Younes.  Geodesic interpolating splines.  In International workshop on energy minimization methods in computer vision and pattern recognition, pages 513–527. Springer, 2001.
  * [8] Y. Cao, M. I. Miller, R. L. Winslow, and L. Younes.  Large deformation diffeomorphic metric mapping of vector fields.  IEEE transactions on medical imaging, 24(9):1216–1230, 2005.
  * [9] J. S. Hesthaven.  From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex.  SIAM Journal on Numerical Analysis, 35(2):655–676, 1998.
  * [10] J. S. Hesthaven and T. Warburton.  Nodal discontinuous Galerkin methods: algorithms, analysis, and applications.  Springer Science & Business Media, 2007.
  * [11] A. Iollo and D. Lombardi.  Advection modes by optimal mass transfer.  Physical Review E, 89(2):022923, 2014.
  * [12] A. Iollo and T. Taddei.  Mapping of coherent structures in parameterized flows by learning optimal transportation with Gaussian models.  Journal of Computational Physics, 471:111671, 2022.
  * [13] M. Jacobs and F. Léger.  A fast approach to optimal transport: The back-and-forth method.  Numerische Mathematik, 146(3):513–544, 2020.
  * [14] G. E. Karniadakis and S. Sherwin.  Spectral/hp element methods for computational fluid dynamics.  Oxford University Press on Demand, 2005.
  * [15] S. G. Krantz and H. R. Parks.  The implicit function theorem: history, theory, and applications.  Springer Science & Business Media, 2002.
  * [16] F. Laudenbach.  About the diffeomorphisms of the 3-sphere and a famous theorem of Cerf ($\Gamma_{4}=0$).  arXiv preprint arXiv:2306.12837, 2023.
  * [17] E. Luke, E. Collins, and E. Blades.  A fast mesh deformation method using explicit interpolation.  Journal of Computational Physics, 231(2):586–601, 2012.
  * [18] J. Ma, J. Wu, J. Zhao, J. Jiang, H. Zhou, and Q. Z. Sheng.  Nonrigid point set registration with robust transformation learning under manifold regularization.  IEEE transactions on neural networks and learning systems, 30(12):3584–3597, 2018.
  * [19] J. E. Marsden and T. J. Hughes.  Mathematical foundations of elasticity.  Courier Corporation, 1994.
  * [20] M. A. Mirhoseini and M. J. Zahr.  Model reduction of convection-dominated partial differential equations via optimization-based implicit feature tracking.  Journal of Computational Physics, page 111739, 2022.
  * [21] R. Mojgani and M. Balajewicz.  Arbitrary Lagrangian Eulerian framework for efficient projection-based reduction of convection dominated nonlinear flows.  In APS Division of Fluid Dynamics Meeting Abstracts, pages M1–008, 2017.
  * [22] I. Mozolevski, E. Süli, and P. R. Bösing.  hp-version a priori error analysis of interior penalty discontinuous galerkin finite element approximations to the biharmonic equation.  Journal of Scientific Computing, 30(3):465–491, 2007.
  * [23] A. Myronenko and X. Song.  Point set registration: Coherent point drift.  IEEE transactions on pattern analysis and machine intelligence, 32(12):2262–2275, 2010.
  * [24] F. Nicoud and F. Ducros.  Subgrid-scale stress modelling based on the square of the velocity gradient tensor.  Flow, turbulence and Combustion, 62(3):183–200, 1999.
  * [25] P.-O. Persson and J. Peraire.  Sub-cell shock capturing for discontinuous Galerkin methods.  In 44th AIAA aerospace sciences meeting and exhibit, page 112, 2006\.
  * [26] G. Peyré and M. Cuturi.  Computational optimal transport: With applications to data science.  Foundations and Trends® in Machine Learning, 11(5-6):355–607, 2019.
  * [27] M. Ruzhansky and M. Sugimoto.  On global inversion of homogeneous maps.  Bulletin of Mathematical Sciences, 5(1):13–18, 2015.
  * [28] N. Sarna and P. Benner.  Data-driven model order reduction for problems with parameter-dependent jump-discontinuities.  Computer Methods in Applied Mechanics and Engineering, 387:114168, 2021.
  * [29] G. R. Shubin, A. Stephens, and H. M. Glaz.  Steady shock tracking and Newton’s method applied to one-dimensional duct flow.  Journal of Computational Physics, 39(2):364–374, 1981.
  * [30] L. Sirovich.  Turbulence and the dynamics of coherent structures. I. Coherent structures.  Quarterly of applied mathematics, 45(3):561–571, 1987.
  * [31] M. L. Staten, S. J. Owen, S. M. Shontz, A. G. Salinger, and T. S. Coffey.  A comparison of mesh morphing methods for 3D shape optimization.  In Proceedings of the 20th international meshing roundtable, pages 293–311. Springer, 2012.
  * [32] T. Taddei.  A registration method for model order reduction: data compression and geometry reduction.  SIAM Journal on Scientific Computing, 42(2):A997–A1027, 2020.
  * [33] T. Taddei.  An optimization-based registration approach to geometry reduction.  arXiv preprint arXiv:2211.10275, 2022.
  * [34] T. Taddei and L. Zhang.  Registration-based model reduction in complex two-dimensional geometries.  Journal of Scientific Computing, 88(3):79, 2021.
  * [35] T. Taddei and L. Zhang.  Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs.  ESAIM: Mathematical Modelling and Numerical Analysis, 55(1):99–130, 2021.
  * [36] P. Tonon, R. A. K. Sanches, K. Takizawa, and T. E. Tezduyar.  A linear-elasticity-based mesh moving method with no cycle-to-cycle accumulated distortion.  Computational Mechanics, 67(2):413–434, 2021.
  * [37] S. Volkwein.  Model reduction using proper orthogonal decomposition.  Lecture Notes, Institute of Mathematics and Scientific Computing, University of Graz. see http://www. uni-graz. at/imawww/volkwein/POD. pdf, 1025, 2011.
  * [38] H. Wendland.  Scattered data approximation, volume 17.  Cambridge university press, 2004.
  * [39] M. J. Zahr and P.-O. Persson.  An optimization-based approach for high-order accurate discretization of conservation laws with discontinuous solutions.  Journal of Computational Physics, 365:105–134, 2018.
  * [40] M. J. Zahr and P.-O. Persson.  An r-adaptive, high-order discontinuous Galerkin method for flows with attached shocks.  In AIAA Scitech 2020 Forum, page 0537, 2020.
  * [41] M. J. Zahr, A. Shi, and P.-O. Persson.  Implicit shock tracking using an optimization-based high-order discontinuous Galerkin method.  Journal of Computational Physics, 410:109385, 2020.