# High-Order Numerical Integration on Domains Bounded by Intersecting Level
Sets

Lauritz Beck1,2, Florian Kummer1,2

( 1 Chair of Fluid Dynamics, Technical University Darmstadt, Germany
2 Graduate School Computational Engineering, Technical University Darmstadt,
Germany
September 3, 2024)

###### Abstract

We present a high-order method that provides numerical integration on volumes,
surfaces, and lines defined implicitly by two smooth intersecting level sets.
To approximate the integrals, the method maps quadrature rules defined on
hypercubes to the curved domains of the integrals. This enables the numerical
integration of a wide range of integrands since integration on hypercubes is a
well known problem. The mappings are constructed by treating the isocontours
of the level sets as graphs of height functions. Numerical experiments with
smooth integrands indicate a high-order of convergence for transformed Gauss
quadrature rules on domains defined by polynomial, rational, and trigonometric
level sets. We show that the approach we have used can be combined readily
with adaptive quadrature methods. Moreover, we apply the approach to
numerically integrate on difficult geometries without requiring a low-order
fallback method.

##  1 Introduction

In engineering, the natural sciences, or economics, partial differential
equations (PDEs) are prevalently used to model complex processes. Because they
facilitate numerical solutions with high accuracy, high-order methods for
solving PDEs have long been a focus of research. Many high-order numerical
methods solving PDEs are based on discretizing the weak formulations of the
PDEs. This requires accurate numerical integration since the weak formulation
involves an integral on the domain of the PDE. The integrand of the weak
formulation belongs to a function space characteristic to the numerical
method. Common choices for high-order methods are polynomial, rational or
trigonometric function spaces [1, 2, 3].

In multi-phase problems like fluid structure interactions or fluid flows with
free surfaces, separate domains are governed by different PDEs. A flexible
approach to model the geometry of the domains is to define the boundaries of
the domains implicitly by the zero-isocontours of smooth level sets. The level
sets often are approximated by functions of the function space of the
numerical method. For example, this approach is applied in level set
methods[4, 5] or Nitsche type methods[6, 7, 8]. Using level sets offers two
advantages. First, the boundaries can be moved by simply advecting the level
sets. Second, a change of topology of the domains is unproblematic for the
model. However, to create a high-order numerical method, this approach has to
be combined with a method that provides high-order numerical integration on
the domains bounded by level sets.

There exist many multi-phase problems where the boundaries of the domains
intersect, forming multi-phase contact lines. A droplet sitting on a soft
substrate features a three-phase contact line situated at the intersection of
the droplet’s surface and the substrate’s surface, for example. To model these
multi-phase problems numerically, it is sensible to represent each surface by
a level set, so that the geometry of the contact line is formed by the
intersection of their zero-isocontours. This is particularly suitable, when
the multi-phase problem contains boundary conditions that require accurate
geometries, for example surface tension[9]. Consequently, discretizing the
problem requires numerical integration on domains that are bounded by
intersecting level sets.

Blending multiple smooth level sets and then applying an integration method
for domains bounded by single level sets is a straightforward approach to
numerically integrate on domains bounded by multiple level sets. In order to
blend multiple level sets, they can be combined into a single level set by
either multiplying them or by choosing their minimum or maximum. Although this
is a viable approach in combination with robust low-order methods, the
resulting level set generally is not smooth, which renders it unfit for most
methods of high-order. Therefore, either a method that can resolve nonsmooth
level sets or a method for intersecting level sets is required.

A variety of high-order methods to integrate on domains bounded by level sets
have been developed. One approach is to extend the domain of integration to a
domain for which a quadrature rule is known. This is achieved by injecting a
function into the integral that models the geometry of the domain. The
injected function can be a dirac delta or a Heaviside function[10, 11, 12], or
it can be defined through a coarea formula[13]. Respectively relying on error
cancelation or a bounded level set, it can be challenging to develop
convergent schemes and importantly, they do not cover nonsmooth geometries or
intersecting level sets.

Another approach relies on reconstructing the surfaces to map a known
quadrature rule to the domain of the integral. While some methods refine an
initial approximation through perturbation and correction [14], other methods
of this approach rely on parametrized curves[15, 16, 17, 18, 19]. Interest in
this approach is enhanced by the fact that it overlaps approaches from
numerical methods with explicit surface representations, see[20, 21] for
example. Some of these methods can resolve kinks, even so reconstructing the
implicit surfaces is cumbersome and bounds the accuracy of the numerical
method by the accuracy of the approximation of the geometry.

When the geometry of the level set surface can be reconstructed, methods
allowing unsmooth surfaces can also be derived via integration by parts. This
has been demonstrated for piecewise parameterized boundaries[22] or trimmed
parametric surfaces[23]. In spite of a favorable order of convergence, they
also suffer from the disadvantages of the reconstruction approach.

Following a different approach, moment fitting methods construct a quadrature
rule by solving a system of equations. The linear moment fitting methods
define quadrature weights through a linear system of equations derived for
preset quadrature nodes. For example, the system of equations can be assembled
by applying the divergence theorem [24, 25], or by utilizing piecewise
parameterized boundaries[26]. Nonlinear methods follow the same approach but
do not require a preset position of the quadrature nodes[27, 28]. Although
these approaches often result in a comparatively low number of quadrature
nodes, they involve the assembly and solution of a system of equations which
is often computationally expensive. In addition, the approaches are limited to
smooth level sets, except for the method relying on reconstruction[26].
Nevertheless, the reconstruction of the parameterized boundary presents a
significant obstacle.

A unique approach is to recast the zero-isocontour as the graph of an implicit
height function [29, 30]. The graph is defined recursively, decomposing the
integral into nested, one dimensional integrals. Although this offers an
efficient high-order integration method for single smooth level sets, some
geometries require a low-order fallback method which can reduce efficiency and
may deteriorate the convergence order of the method. If the level set is
restricted to polynomials, a method without the need for a low-order fallback
was presented recently that can resolve geometries of intersecting level
sets[31]. However, a quadrature method for domains bounded by general, smooth
intersecting level sets is missing.

In this paper we propose a high-order method that provides numerical
integration on volumes, surfaces, and lines defined implicitly by two
intersecting level sets. Augmenting the approach of Saye[29], the method is
based on implicit height functions and restricts the level sets to smooth
function spaces only. The implicit height functions are used to map hypercubes
to the curved domains of the integrals, offering a general approach to
construct quadrature rules. Since we conceptually integrate on hypercubes,
adaptive quadrature methods can be applied readily. Moreover, the approach can
be applied to integrals defined by single level sets, removing the need for a
low-order fallback method. The method is also available as source code which
is publicly available on github[32].

This paper is organized as follows: After explaining the principle of the
method, section 2 introduces its main component, nested mappings. Section 3
details the two algorithms involved in creating nested mappings before
utilizing them to derive numerical quadrature rules in section 4. In section
5, numerical experiments are conducted and discussed. Finally, section 6
briefly summarizes our results and concludes with further discussions.

##  2 The Gist of it

Figure 1:  Implicit definition of the domains of integration. The level sets
$\alpha$ and $\beta$ partition hyperrectangle $K$ (a, b), defining the
surfaces $\Gamma_{\alpha},\Gamma_{\beta}$, and the boundary lines
$\partial\Gamma_{\alpha},\partial\Gamma_{\beta}$ as subsets of their zero-
isocontours (c). The region with negative sign of $\alpha$ and $\beta$
specifies the domain of integration $\Omega$ (d).

We are interested in numerically integrating a function
$g(\boldsymbol{x}),\mathbb{R}^{3}\rightarrow\mathbb{R},$ on an implicitly
defined three-dimensional domain $\Omega\subset\mathbb{R}^{3}$,

$\int_{\Omega}g(\boldsymbol{x})\,\text{d}V\approx\sum_{i}w_{i}g(\boldsymbol{x}_{i}),$
(1)

with a quadrature rule, consisting of a set of nodes
$\boldsymbol{x}_{i}\in\mathbb{R}^{d}$ and corresponding weights
$w_{i}\in\mathbb{R}$. The domain $\Omega$ lies inside a three-dimensional
hyperrectangle $K$ and is the intersection of the domains $\Omega_{\alpha}$
and $\Omega_{\beta}$,

$\Omega=\Omega_{\alpha}\cap\Omega_{\beta},\qquad\Omega_{\alpha}=\\{\boldsymbol{x}\in
K|\alpha(\boldsymbol{x})\leq 0\\},\quad\Omega_{\beta}=\\{\boldsymbol{x}\in
K|\beta(\boldsymbol{x})\leq 0\\},$ (2)

defined by the sign of two level set functions
$\alpha(\boldsymbol{x}):K\rightarrow\mathbb{R}$ and
$\beta(\boldsymbol{x}):K\rightarrow\mathbb{R}$ (see Figure 1). Further, we are
interested in numerically integrating on the two-dimensional smooth surfaces
$\gamma_{\alpha}$ and $\gamma_{\beta}$,

$\oint_{\gamma\in\\{\gamma_{\alpha},\gamma_{\beta}\\}}g(\boldsymbol{x})\,\text{d}S,\quad\gamma_{\alpha}=\Gamma_{\alpha}\cap\Omega_{\beta},\quad\gamma_{\beta}=\Gamma_{\beta}\cap\Omega_{\alpha},$
(3)

defined by the zero-isocontours of the level set functions,

$\Gamma_{\alpha}=\\{\boldsymbol{x}\in
K|\alpha(\boldsymbol{x})=0\\},\quad\Gamma_{\beta}=\\{\boldsymbol{x}\in
K|\beta(\boldsymbol{x})=0\\},$ (4)

and on their one-dimensional smooth intersection $\partial\gamma$,

$\fint_{\partial\gamma}g(\boldsymbol{x})\,\text{d}L,\quad\partial\gamma=\Gamma_{\alpha}\cap\Gamma_{\beta}.$
(5)

###  2.1 Composition

The general idea is to simplify the domain of the integral by transforming it
into a hypercube. In case of the volume integral (1), we map a hypercube
$C=[-1,1]^{3}$ to $\Omega$ with a mapping $T$,

$T:C\rightarrow\Omega,\quad T(\tilde{x},\tilde{y},\tilde{z})=(x,y,z),$ (6)

encoding $\Omega$’s geometry in the Jacobian determinant
$|\text{det}(DT(\tilde{x}))|=J(T(\tilde{x}))$,

$\int_{\Omega}g(\boldsymbol{x})\,\text{d}V=\int_{C}g(T(\boldsymbol{\tilde{x}}))\,|\text{det}(DT(\boldsymbol{\tilde{x}}))|\,\text{d}V=\int_{C}g(T(\boldsymbol{\tilde{x}}))\,J(T(\boldsymbol{\boldsymbol{\tilde{x}}}))\,\text{d}V.$
(7)

Since it is a cartesian product of intervals, a range of quadrature rules
suits $C$, e.g. tensorized Gaussian quadrature
$\boldsymbol{\tilde{x}}_{i},\tilde{w}_{i}$,

$\int_{K}g(T(\boldsymbol{\tilde{x}}))\,J(T(\boldsymbol{\tilde{x}}))\,\text{d}V\approx\sum_{i}\tilde{w}_{i}g(T(\boldsymbol{\tilde{x}}_{i}))J(T(\boldsymbol{\tilde{x}}_{i})).$
(8)

After transforming back to $\Omega$,

$\int_{\Omega}g(\boldsymbol{x})\,\text{d}V\approx\sum
w_{i}g(\boldsymbol{x}_{i}),\quad\boldsymbol{x}_{i}=T(\boldsymbol{\tilde{x}}_{i}),\quad
w_{i}=J(\boldsymbol{x}_{i})\tilde{w}_{i}.$ (9)

we receive a quadrature rule depending only on $T$.

Figure 2: Different steps involved in computing a quadrature rule on domain
$\Omega$, colored in dark green (1). After mapping $\Omega$ to a hypercube $C$
in two steps (1-2), the integral may be evaluated via standard tensorized
quadrature, e.g. tensorized Gaussian quadrature $\hat{Q}$ (3). By embedding
$\hat{Q}$ in $\Omega$ with mappings $A$ and $B$ through transforming back (4,
5), the final quadrature rule $Q$ is generated.

We define $T$ as a composition of mappings for single surfaces, $T=A\circ B$,
utilizing that $\Omega$ is bounded by $\Gamma_{\alpha}$ and $\Gamma_{\beta}$.
As illustrated in Figure 2, we proceed in five steps:

  1. 1.

Find a mapping
$A:C\rightarrow\Omega_{\alpha},A(\tilde{x},\tilde{y},\tilde{z})=(x,y,z)$ and
map the domain $\Omega_{\alpha}$ defined by level set $\alpha$ to the
hypercube $C=[-1,1]^{3}$, s.t. $\Omega_{\alpha}=A(C)$. Map level set $\beta$
to the hypercube to receive
$\tilde{\beta}(\tilde{x},\tilde{y},\tilde{z})=\beta(A(\tilde{x},\tilde{y},\tilde{z}))$.

  2. 2.

Find a mapping
$B:C\rightarrow\Omega_{\tilde{\beta}},B(\hat{x},\hat{y},\hat{z})=(\tilde{x},\tilde{y},\tilde{z}),$
and map the domain $\Omega_{\tilde{\beta}}=\\{\boldsymbol{\tilde{x}}\in
C|\tilde{\beta}(\boldsymbol{\tilde{x}})<0\\}$ defined by the transformed level
set $\tilde{\beta}$ to the hypercube $C=[-1,1]^{3}$, s.t.
$\Omega_{\tilde{\beta}}=B(C)$.

  3. 3.

Solve the transformed integral numerically on the hypercube $C$

$\int_{\Omega}g(\boldsymbol{x})\,\text{d}V=\int_{C}g(A(B(\boldsymbol{\hat{x}})))\,|\text{det}(DA(B(\boldsymbol{\hat{x}})))|\,\text{d}V\approx\sum_{i}\hat{w}_{i}g(T(\boldsymbol{\hat{x}}_{i}))J_{A}(T(\boldsymbol{\hat{x}}_{i}))J_{B}(B(\boldsymbol{\hat{x}}_{i}))$
(10)

using a quadrature rule with nodes
$\hat{Q}=\\{\boldsymbol{\hat{x}}_{1},\boldsymbol{\hat{x}}_{2},...\\}$ and
weights $\hat{w}_{i}$. The Jacobian determinant of the composition $T$ is
separated into its components

$|\text{det}(DA(B(\boldsymbol{\hat{x}})))|=|\text{det}(DA\,DB)|=J_{A}(T(\boldsymbol{\hat{x}}))J_{B}(B(\boldsymbol{\hat{x}})).$
(11)

  4. 4.

Transform back to receive quadrature nodes
$\tilde{Q}=\\{\boldsymbol{\tilde{x}}_{1},\boldsymbol{\tilde{x}}_{2},...\\}=B(\hat{Q})$
and weights $\tilde{w}_{i}=\hat{w}_{i}J_{B}(\boldsymbol{\tilde{x}}_{i})$,

$\sum_{i}\hat{w}_{i}g(T(\boldsymbol{\hat{x}}_{i}))J_{B}(B(\boldsymbol{\hat{x}}_{i}))J_{A}(T(\boldsymbol{\hat{x}}_{i}))=\sum_{i}\hat{w}_{i}g(A(\boldsymbol{\tilde{x}}_{i}))J_{B}(\boldsymbol{\tilde{x}}_{i})J_{A}(A(\boldsymbol{\tilde{x}}_{i}))=\sum_{i}\tilde{w}_{i}g(A(\boldsymbol{\tilde{x}}_{i}))J_{A}(A(\boldsymbol{\tilde{x}}_{i})).$
(12)

  5. 5.

Transform back quadrature nodes
$Q=\\{\boldsymbol{x}_{1},\boldsymbol{x}_{2},...\\}=A(\tilde{Q})$, yielding the
final quadrature rule,

$\int_{\Omega}g(\boldsymbol{x})\,\text{d}V\approx\sum_{i}\tilde{w}_{i}g(A(\boldsymbol{\tilde{x}}_{i}))J_{A}(A(\boldsymbol{\tilde{x}}_{i}))=\sum_{i}\tilde{w}_{i}g(\boldsymbol{x}_{i})J_{A}(\boldsymbol{x}_{i})=\sum_{i}w_{i}g(\boldsymbol{x}_{i}),$
(13)

with nodes $\boldsymbol{x}_{i}$ and scaled weights
$w_{i}=J_{A}(\boldsymbol{x}_{i})\tilde{w}_{i}$.

By simply composing mappings, constructing a quadrature rule for domains
bounded by two level set isocontours is reduced to repeatedly mapping a
hypercube to domains bounded by a single level set isocontour. The question
that remains is: how to find such mappings?

###  2.2 Nested Mapping

Representing the backbone of the approach, let the nested mapping
$M:C\rightarrow\Omega_{M}$,

$M(\tilde{x},\tilde{y},\tilde{z})=\frac{1}{2}\begin{pmatrix}(m_{x}^{\uparrow}-m_{x}^{\downarrow})\tilde{x}+m_{x}^{\uparrow}+m_{x}^{\downarrow}\\\
(m_{y}^{\uparrow}(x)-m_{y}^{\downarrow}(x))\tilde{y}+m_{y}^{\uparrow}(x)+m_{y}^{\downarrow}(x)\\\
(m_{z}^{\uparrow}(x,y)-m_{z}^{\downarrow}(x,y))\tilde{z}+m_{z}^{\uparrow}(x,y)+m_{z}^{\downarrow}(x,y)\end{pmatrix}=\begin{pmatrix}x\\\
y\\\ z\end{pmatrix},$ (14)

be a composition of height functions,

$m_{x}^{\uparrow},\quad m_{x}^{\downarrow},\quad
m_{y}^{\uparrow}=m_{y}^{\uparrow}(x),\quad
m_{y}^{\downarrow}=m_{y}^{\downarrow}(x),\quad
m_{z}^{\uparrow}=m_{z}^{\uparrow}(x,y),\quad
m_{z}^{\downarrow}=m_{z}^{\downarrow}(x,y),$ (15)

that maps a three-dimensional hypercube $C$ to the codomain $\Omega_{M}$.
Nesting its components, $M$ is defined recursively: Evaluating the second
entry of $M$ requires its first entry and evaluating the third entry of $M$
requires its second entry.

Figure 3: The nested mapping $M$ maps the hypercube $C$ to codomain
$\Omega_{M}$ (a). The codomain is defined by a set of height functions
$m_{x}^{\uparrow}$, $m_{x}^{\downarrow}$, $m_{y}^{\uparrow}(x)$,
$m_{y}^{\downarrow}(x)$, $m_{z}^{\uparrow}(x,y)$, $m_{z}^{\downarrow}(x,y)$,
confining $\Omega_{M}$ along each coordinate axis (b).

Its height functions do not intersect

$\displaystyle m_{x}^{\downarrow}<m_{x}^{\uparrow},\quad
m_{y}^{\downarrow}<m_{y}^{\uparrow},\quad
m_{z}^{\downarrow}<m_{z}^{\uparrow},$ (16)

and confine $\Omega_{M}$ along the three coordinate axes as shown in Figure 3.
First, the constant height functions $m_{x}^{\uparrow}$ and,
$m_{x}^{\downarrow}$ confine $\Omega_{M}$ along the $x$-axis. Second, the
univariate height functions $m_{y}^{\downarrow}(x)$ and $m_{y}^{\uparrow}(x)$
confine $\Omega_{M}$ along the $y$-axis and, third, the bivariate height
functions $m_{z}^{\uparrow}(x,y)$ and $m_{z}^{\downarrow}(x,y)$ confine
$\Omega_{M}$ along the $z$-axis.

Inheriting the nested mappings’s properties, the Jacobian matrix $DM$,

$DM=\frac{1}{2}\begin{pmatrix}m^{\uparrow}_{x}-m^{\downarrow}_{x}&0&0\\\
M_{xy}\partial_{\tilde{x}}x&m_{y}^{\uparrow}-m_{y}^{\downarrow}&0\\\
M_{xz}\partial_{\tilde{x}}x+M_{yz}\partial_{\tilde{x}}y&M_{yz}\partial_{\tilde{y}}y&m_{z}^{\uparrow}-m_{z}^{\downarrow}\end{pmatrix}=\begin{pmatrix}\partial_{\tilde{x}}x&0&0\\\
\partial_{\tilde{x}}y&\partial_{\tilde{y}}y&0\\\
\partial_{\tilde{x}}z&\partial_{\tilde{y}}z&\partial_{\tilde{z}}z\end{pmatrix},$
(17)

is structured recursively. Its entries, $M_{xy},M_{xz},M_{yz}$,

$\displaystyle
M_{xy}=\partial_{x}y(x,\tilde{y})=(\partial_{x}m_{y}^{\uparrow}-\partial_{x}m_{y}^{\downarrow})\tilde{y}+\partial_{x}m_{y}^{\uparrow}+\partial_{x}m_{y}^{\downarrow},$
(18) $\displaystyle
M_{xz}=\partial_{x}z(x,y,\tilde{z})=(\partial_{x}m_{z}^{\uparrow}-\partial_{x}m_{z}^{\downarrow})\tilde{z}+\partial_{x}m_{z}^{\uparrow}+\partial_{x}m_{z}^{\downarrow},$
(19) $\displaystyle
M_{yz}=\partial_{y}z(x,y,\tilde{z})=(\partial_{y}m_{z}^{\uparrow}-\partial_{y}m_{z}^{\downarrow})\tilde{z}+\partial_{y}m_{z}^{\uparrow}+\partial_{y}m_{z}^{\downarrow},$
(20)

require the first partial derivatives of the height functions. However, to
evaluate the Jacobian determinant $J_{M}$,

$J_{M}(\tilde{x},\tilde{y},\tilde{z})=\det
DM=\frac{1}{8}(m^{\uparrow}_{x}-m^{\downarrow}_{x})(m_{y}^{\uparrow}-m_{y}^{\downarrow})(m_{z}^{\uparrow}-m_{z}^{\downarrow}),$
(21)

partial derivatives are not required.

To increase the possible shapes of $\Omega_{M}$, the mapping can be oriented
differently by rearranging $M$,

$M^{\prime}=Q^{T}M(Q(\tilde{x},\tilde{y},\tilde{z})),\quad
Q\in\\{0,1\\}^{3\times 3},\quad\det(Q)=1,$ (22)

with a permutation matrix $Q$. By rearranging, $\Omega_{M}$ can be confined
along the $x$-axis by a bivariate surface for example. Naturally, the mapping
can be restricted to two-dimensional sets by dropping component $z$.

###  2.3 Construction

Our approach to construct a mapping, $A:C\rightarrow\Omega_{\alpha}$, is to
match the codomain of a nested mapping to $\Omega_{\alpha}$. To match
$\Omega_{\alpha}$, we disassemble its surface $\Gamma_{\alpha}$ into a set of
height functions. This is possible if $\Gamma_{\alpha}$ and the intersections
of $\Gamma_{\alpha}$ with suitable faces, edges and corners of the enclosing
hyperrectangle $K$, are graphs on their associated coordinate planes.

To illustrate, let $\Omega_{\alpha}$ be a subset of the cube
$K^{3}=[-1,1]^{3}$, such that its surface is a graph,

$\Gamma_{\alpha}=\left\\{\gamma(\boldsymbol{x})\,|\,\boldsymbol{x}\in\Pi_{z}(\Gamma_{\alpha})\right\\},\quad\gamma(\boldsymbol{x})=\begin{pmatrix}x\\\
y\\\ a_{z}^{\uparrow}(x,y)\end{pmatrix},$ (23)

on its orthogonal projection $\Pi_{z}(\Gamma_{\alpha})$ onto the $xy$-plane
(see figure 4). Hence, by embedding $a^{\uparrow}_{z}(x,y)$ in $K^{3}$, the
graph’s function determines the first height function of the nested mapping.
Forming the graph’s domain, the orthogonal projection flattens direction
$h\in\\{x,y,z\\}$,

$\Pi_{h}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d},\quad\Pi_{h}(\boldsymbol{x})=(\mathbb{1}-\boldsymbol{e}_{h}\otimes\boldsymbol{e}_{h})\boldsymbol{x}.$
(24)

The projection utilizes height vector $\boldsymbol{e}_{h}$ which is $1$ at
component $h$ and $0$ otherwise.

Further, let the intersection of the surface and the cube’s lower face $K^{2}$
be a graph on the orthogonal projection of the intersection onto the
$xz$-plane,

$\Gamma_{\alpha}\cap
K^{2}=\left\\{\gamma(\boldsymbol{x})\,|\,\boldsymbol{x}\in\Pi_{y}(\Gamma_{\alpha}\cap
K^{2})\right\\},\quad\gamma(\boldsymbol{x})=\begin{pmatrix}x\\\
a_{y}^{\downarrow}(x)\\\ z\end{pmatrix}.$ (25)

Repeating the approach, the graph’s function $\gamma$ subsequently specifies
the second height function of the nested mapping, $a^{\downarrow}_{y}(x)$.

Finally, the graph of the intersection of the cube’s edge $K^{1}$ and
$\Gamma_{\alpha}$ on its projection onto the $yz$-plane,

$\Gamma_{\alpha}\cap
K^{1}=\left\\{\gamma(\boldsymbol{x})\,|\,\boldsymbol{x}=\Pi_{x}(\Gamma_{\alpha}\cap
K^{1})\right\\},\quad\gamma(\boldsymbol{x})=\begin{pmatrix}a_{x}^{\uparrow}\\\
y\\\ z\end{pmatrix},$ (26)

determines the height function $a^{\uparrow}_{x}$.

When combined, the height functions define the nested mapping,

$A(\tilde{x},\tilde{y},\tilde{z})=\frac{1}{2}\begin{pmatrix}(a^{\uparrow}_{x}+1)\tilde{x}+a^{\uparrow}_{x}-1\\\
(1-a^{\downarrow}_{y}(x))\tilde{y}+1+a^{\downarrow}_{y}(x)\\\
(a^{\uparrow}_{z}(x,y)+1)\tilde{z}+a^{\uparrow}_{z}(x,y)-1\end{pmatrix}=\begin{pmatrix}x\\\
y\\\ z\end{pmatrix}.$ (27)

where each axis is respectively confined by a height function and a face, edge
or point of $K^{3}$.

Figure 4: Disassembly of the implicit surface $\Gamma_{\alpha}$ of domain
$\Omega_{\alpha}$ into graphs of the height functions
$a_{x}^{\uparrow},a_{y}^{\downarrow}$ and $a_{z}^{\uparrow}$, from which a
nested mapping can be constructed (a): Height function $a_{z}^{\uparrow}(x,y)$
graphs surface $\Gamma_{\alpha}\cap K^{3}$ over the projection of the surface
onto the $xy$-plane (b). Height function $a_{y}^{\downarrow}$ graphs surface
line $\Gamma_{\alpha}\cap K^{2}$ over the projection of the surface line onto
of the $xz$-plane (c). Height function $a_{x}^{\uparrow}$ graphs
$\Gamma_{\alpha}\cap K^{1}$ over a point of the $yz$-plane (d).

##  3 Algorithmic Assembly

In view of the fact that only mappings for domains bounded by a single level
set are required to construct quadrature rules, this section concentrates on
mapping $C$ to $\Omega_{\alpha}$. Algorithmically, two separate components are
involved in finding suitable mappings. The first component of the algorithm
partitions $K$ into suitable hyperrectangles, ensuring that $\Gamma_{\alpha}$
is representable by a graph in each. The second component defines a set of
nested mappings by tessellating each hyperrectangle with the nested mappings’
codomains.

Both components rely on nested bodies; A data structure that allows the clear
geometric decomposition of a domain bounded by a graph. A nested body
$G=(K_{G},H)$ is a strict rooted binary tree combining a set of bodies $K_{G}$
with a set of height directions $H$. It hierarchically structures the bodies
which are compact sets embedded in $\mathbb{R}^{d}$, e.g. a square embedded in
$\mathbb{R}^{3}$. To help distinguish between them, all bodies of equal
dimension $i$ are collected in the sets $\mathbb{K}^{i}$. Precisely, a nested
body $G$ has the following properties:

  1. 1.

Each node of the binary tree is a body $k\in K_{G}$.

  2. 2.

A parent node $k\in\mathbb{K}^{i}$ with dimension $i$ has two child nodes
$k^{\downarrow},k^{\uparrow}\in\mathbb{K}^{i-1}$ with dimension $i-1$.

  3. 3.

Each node $k$ of level $l$ of the binary tree is a $d-l$-dimensional body,
$k\in\mathbb{K}^{d-l}$.

  4. 4.

The nodes of the deepest feasible level $l=d$ consist of zero-dimensional
points.

  5. 5.

Each level $l<d$ corresponds to a height direction $h_{d-l}\in\\{x,y,...\\}$.

  6. 6.

The set $H$ gathers the height directions of every level,
$H=\\{h_{d},h_{d-1},...\\}$.

###  3.1 Graph of the Surface

Algorithm 1 recursively creates a nested body $G=(K_{G},H)$ consisting of
hyperrectangles. Beginning with a nested body consisting of a single level
containing hyperrectangle $K$, it gradually deepens $G$ by selecting suitable
height directions. Corresponding to the selected height direction, the
algorithm adds a new level of hyperrectangles to $G$ while upholding the
invariant: For each $k\in K_{G}$ with associated height direction $h\in H$,
$\Gamma_{\alpha}\cap k$ is a graph on $\Pi_{h}(\Gamma_{\alpha}\cap k)$.

In each recursion, Algorithm 1 constructs a set of hyperrectangles $P$. After
gathering the lowest level of hyperrectangles of $G$ in $P$, it removes each
hyperrectangle $k$ from $P$ satisfying one of two conditions: First, if the
hyperrectangle does not intersect $\Gamma_{\alpha}$, $\Gamma_{\alpha}\cap
k=\emptyset$. Second, if the entire hyperrectangle is a subset of
$\Gamma_{\alpha}$, $\Gamma_{\alpha}\cap k=k$. It removes them because the
invariant holds unrestrictedly for these hyperrectangles and all contained
faces, regardless of height direction. For all remaining hyperrectangles $k\in
P$, Algorithm 1 chooses a common height direction $h$ and checks if the
invariant holds. If it cannot be guaranteed, Algorithm 1 subdivides $G$ and
restarts or designates $G$ for a brute force method. If the invariant can be
guaranteed, it attaches the faces of $P$ as a new level of hyperrectangles to
$G$.

Each new level of $G$ consists of the lower faces and the upper faces $F$ in
direction $h$ of the hyperrectangles $P$,

$F=\text{face}^{\downarrow}_{h}(P)\cup\text{face}^{\uparrow}_{h}(P).$ (28)

The lower face, $\text{face}_{h}^{\downarrow}(k)$, and upper face,
$\text{face}_{h}^{\uparrow}(k)$, of a hyperrectangle $k\in\mathbb{K}^{i},i>0,$

$\displaystyle\text{face}_{h}^{\downarrow}:\quad\mathbb{K}^{i}\rightarrow\mathbb{K}^{i-1},\quad\text{face}_{h}^{\downarrow}(k)$
$\displaystyle=\left\\{\boldsymbol{x}\in
k\,|\,\boldsymbol{x}\cdot\boldsymbol{e}_{h}\leq\boldsymbol{y}\cdot\boldsymbol{e}_{h},\forall\boldsymbol{y}\in
k)\right\\},$ (29)
$\displaystyle\text{face}_{h}^{\uparrow}:\quad\mathbb{K}^{i}\rightarrow\mathbb{K}^{i-1},\quad\text{face}_{h}^{\uparrow}(k)$
$\displaystyle=\left\\{\boldsymbol{x}\in
k\,|\,\boldsymbol{x}\cdot\boldsymbol{e}_{h}\geq\boldsymbol{y}\cdot\boldsymbol{e}_{h},\forall\boldsymbol{y}\in
k)\right\\},$ (30)

are subsets of the hyperrectangle’s boundary.

As long as $P$ contains hyperrectangles new levels are attached to $G$ until
its deepest level consists of points. When Algorithm 1 has terminated,
$\Gamma_{\alpha}$ is a graph whose domain is either a hyperrectangle or a
subset of a hyperrectangle bounded by graphs. If the domain is bounded by
graphs, their domains are in turn either hyperrectangles or subsets of
hyperrectangles bounded by graphs, and so on. Illustrating the involved
hyperrectangles, Figure 5 provides an example of a nested body created by
Algorithm 1 for a planar $\Gamma_{\alpha}$.

Algorithm 1 Creates a nested body by gradually adding new levels of faces to
it. The surface $\Gamma_{\alpha}$ is a graph on the hyperrectangles of the
nested body.

1:Nested body $G=(K_{G},H)$ and level set $\alpha$. $G$ consists of at least
one level of hyperrectangles.

2:Nested body $G$ consisting of hyperrectangles. For each hyperrectangle $k\in
K_{G}$ with height direction $h\in H$ which is not a point,
$\Gamma_{\alpha}\cap k$ is a graph on $\Pi_{h}(\Gamma_{\alpha}\cap k)$.

3:procedure Graph($G$, $\alpha$)

4: $P\leftarrow$ the set of hyperrectangles of the deepest level of $G$

5: Remove hyperrectangles $\\{k\in P|k\cap\Gamma_{\alpha}=k\lor
k\cap\Gamma_{\alpha}=\emptyset\\}$ from $P$

6: if $|P|>0$ then

7: $h\leftarrow$ height direction of $P$

8: for all $k\in P$ do

9: if $\Gamma_{\alpha}\cap k$ is not a graph over $\Pi_{h}(\Gamma_{\alpha}\cap
k)$ then

10: if Subdivision count below threshold then

11: Split the nested body, $G\rightarrow G^{\prime},G^{\prime\prime}$

12: Restart and call Graph($G^{\prime},\alpha$) and
Graph($G^{\prime\prime},\alpha$)

13: else

14: Desginate $G$ for a brute force method and return

15: end if

16: end if

17: end for

18: Set the height direction of the deepest level of $G$ to $h$

19: $F\leftarrow$
$\text{face}_{h}^{\downarrow}(P)\cup\text{face}_{h}^{\uparrow}(P)$

20: Add $F$ as new deepest level to $G$

21: Graph($G,\alpha$)

22: end if

23:end procedure

Figure 5: A nested body $G$ created by Algorithm 1 for a planar
$\Gamma_{\alpha}$, hierarchically linking a hyperrectangle to its faces.
Creating a strict rooted binary tree, the faces of each body are attached to
each body unless a body has a sign $\in\\{-,0,+\\}$ imposed by level set
$\alpha$. The bodies of each level are successively defined by the level’s
height direction $h_{i}$.

So far, we have ommitted the inner workings of the algorithm’s components
which we will discuss in the following. To decide if a hyperrectangle should
be removed from $P$, we assess if $\Gamma_{\alpha}\cap k=\emptyset$ or
$\Gamma_{\alpha}\cap k=k$ on the hyperrectangle $k$. This is equivalent to
evaluating the sign $\in\\{-,0,+,\pm\\}$ of the level set on $k$,

$\Gamma_{\alpha}\cap
k=\emptyset\Leftrightarrow\text{sign}_{\alpha}(k)\in\\{-,+\\},\quad\Gamma_{\alpha}\cap
k=k\Leftrightarrow\text{sign}_{\alpha}(k)=0.$ (31)

The sign,

$\text{sign}_{\alpha}(k)=\begin{cases}-,\quad&\text{if
}\min\limits_{\boldsymbol{x}\in
k}\alpha(\boldsymbol{x})<0\land\max\limits_{\boldsymbol{x}\in
k}\alpha(\boldsymbol{x})\leq 0\\\ 0,\quad&\text{if
}\min\limits_{\boldsymbol{x}\in
k}\alpha(\boldsymbol{x})=0\land\max\limits_{\boldsymbol{x}\in
k}\alpha(\boldsymbol{x})=0\\\ +,\quad&\text{if }\min\limits_{\boldsymbol{x}\in
k}\alpha(\boldsymbol{x})\geq 0\land\max\limits_{\boldsymbol{x}\in
k}\alpha(\boldsymbol{x})>0\\\ \pm,\quad&\text{else}\end{cases},$ (32)

indicates if $\alpha$ can be uniformly bounded on a set $k$, not necessarily a
hyperrectangle.

To evaluate the sign numerically, we approximate $\alpha(k)$ with a Bézier
patch of degree $n$,

$\alpha(k)\approx p(k),\quad
p(\boldsymbol{x})=\sum_{i}^{m}q_{i}b_{i}(\boldsymbol{u}(\boldsymbol{x})),\quad\boldsymbol{u}:k\rightarrow[-1,1]^{\dim(k)}$
(33)

with $m=(n+1)^{\dim(k)}$ control points $q_{i}$ and basis functions $b_{i}$.
The basis functions are Bernstein polynomials, tensorized to match the
dimension of $k$. The control points are chosen so that $p$ interpolates
$\alpha$ at equidistant data points, spanning a mesh including the boundary of
$k$. Since the Bézier patch lies in their convex hull, the control points
bound $p(k)$,

$\min_{i}q_{i}\leq\min_{\boldsymbol{x}\in
k}p(\boldsymbol{x})\leq\max_{\boldsymbol{x}\in
k}p(\boldsymbol{x})\leq\max_{i}q_{i}.$ (34)

Utilizing the bound, we estimate the sign of $k$,

$\text{sign}_{\alpha}(k)\approx\text{sign}_{p}(k)\approx\begin{cases}-,\quad&\text{if
}\min\limits_{i}q_{i}<-\epsilon\land\max\limits_{i}q_{i}\leq\epsilon\\\
0,\quad&\text{if
}\min\limits_{i}q_{i}\geq-\epsilon\land\max\limits_{i}q_{i}\leq\epsilon\\\
+,\quad&\text{if
}\min\limits_{i}q_{i}\geq-\epsilon\land\max\limits_{i}q_{i}>\epsilon\\\
\pm,\quad&\text{else}\end{cases}.$ (35)

with a precision threshold $\epsilon$, necessitated by the limited precision
of machine numbers. Relying on the convex hull of the approximation, the
estimate may wrongly assess the body’s sign to be undetermined, $\pm$, in some
cases. Fortunately, this only causes the algorithm to add unnecessary faces to
$G$ without compromising the invariant. In rare cases, the estimate may
wrongly determine a sign $\in\\{-,0,+\\}$ if $\alpha$ lies outside the convex
hull of $p$ which could compromise correctness. However, practice has shown
that $n=3$ suffices to reliably assign a sign to a body in a sufficiently
resolved numerical grid.

As a way to determine a suitable height direction $h$ for the set of
hyperrectangles $P$,

$h=\underset{j}{\text{argmax}\,}\boldsymbol{y}_{j},\quad\boldsymbol{y}=\frac{1}{|P|}\sum_{k\in
P}\nabla_{k}\alpha(\boldsymbol{c}_{k}),$ (36)

the index of the greatest entry of $\boldsymbol{y}$ is selected. The vector
$\boldsymbol{y}$ denotes the average of the projections of each gradient
$\nabla_{k}\alpha$ onto each hyperrectangle $k\in P$. Each projected gradient
is evaluated at the center $\boldsymbol{c}_{k}$ of its associated
hyperrectangle.

To decide if $\Gamma_{\alpha}\cap k$ is a graph over
$\Pi_{h}(\Gamma_{\alpha}\cap k)$, we bound the gradient $\partial_{h}\alpha$
of the level set. It is a graph, if the component $h$ of $\partial_{h}\alpha$
does not change its sign on a hyperrectangle $k$,

$\text{sign}_{\partial_{h}\alpha}(k)\neq\pm.$ (37)

Once again, we interpolate $\partial_{h}\alpha$ with a Bézier patch to
evaluate its sign. Offering a conservative approximation, this estimate may
deem surfaces, which are representable as graphs, unrepresentable and lead to
unneeded subdivisions. However, it does not incorrectly evaluate the surface
as representable as a graph which would contradict the invariant.

Finally, if $\Gamma_{\alpha}\cap k$ is not a graph over
$\Pi_{h}(\Gamma_{\alpha}\cap k)$, the algorithm either subdivides $K$ and
restarts or employs a brute force method, depending on the number of previous
subdivisions. Since a linear polynomial is always a graph over the faces of a
hyperrectangle, approximating $\alpha$ by a linear polynomial represents a
convenient brute force method. If combined with a sufficiently large number of
subdivisions, high accuracy can be maintained.

###  3.2 Tessellation of the Domain

Algorithm 2 tesselates the root body $R$ of a nested body $G$ created by
Algorithm 1,

$R=\dot{\bigcup}V,\quad
V=\left\\{v|\text{sign}_{\alpha}(v)\in\\{-,+\\}\right\\},$ (38)

so that each tile $v\in V$ is the codomain of a nested mapping and is either a
subset of $R\cap\Omega_{\alpha}$ or $R\setminus\Omega_{\alpha}$. In order to
create $V$, Algorithm 2 recursively splits $R$ into suitable subsets which are
structured in nested bodies. Importantly, the subset $V^{-}$ of $V$,

$V^{-}=\left\\{v\in V|\text{sign}_{\alpha}(v)=-\right\\},$ (39)

is directly required for numerical integration since it tesselates the domain
of the integral $R\cap\Omega_{\alpha}$.

Algorithm 2 Tessellates a hyperrectangle with the codomains of nested mappings
and structures each codomain in a nested body $W$. The algorithm creates the
nested body by successively assigning signs to the nodes of $W$, deriving the
sign of a node from its children or by splitting it along $\Gamma_{\alpha}$.

1:Nested body $W=(K_{W},H)$ and level set $\alpha$. Every node $w\in K_{W}$
with height direction $h$ and children $w^{\uparrow}$ and $w^{\downarrow}$
satisfies the following:

  1. 1.

If $\text{sign}_{\alpha}(w)=\pm$ , then
$\text{sign}_{\alpha}(w^{\downarrow})\neq\pm$ and
$\text{sign}_{\alpha}(w^{\uparrow})\neq\pm.$

  2. 2.

$\Gamma_{\alpha}\cap w$ is a graph on $\Pi_{h}(\Gamma_{\alpha}\cap w)$.

  3. 3.

$\Pi_{h}(w)=\Pi_{h}(w^{\downarrow})=\Pi_{h}(w^{\uparrow})$.

2:Nested body $W$. Every node $w\in K_{W}$ has $\text{sign}(w)\in\\{-,0,+\\}$.

3:procedure Tessellate($W$, $\alpha$)

4: for all sets of domains $K^{i}_{W}$ of $W$ from second deepest level to
level 0 do

5: for all $w\in K^{i}_{W}$ with sign $\pm$ do

6: $(w,w^{\downarrow},w^{\uparrow})\leftarrow\left(w,\,\text{first child of
}w,\,\text{second child of }w\right)$

7: if the sign of $w^{\downarrow}$ is $+$ and the sign of $w^{\uparrow}$ is
$-$, or vice versa then

8: Split the nested domain $W\rightarrow W^{\prime},W^{\prime\prime}$ along
$\Gamma_{\alpha}\cap w$ so that $w=w^{\prime}\dot{\cup}w^{\prime\prime}$

9: Update height functions of $W^{\prime}$ and $W^{\prime\prime}$

10: Sign of $w^{\prime}\leftarrow$ sign of $w^{\downarrow}$

11: Sign of $w^{\prime\prime}\leftarrow$ sign of $w^{\uparrow}$

12: Restart and call Tessellate($W^{\prime},\alpha$) and
Tessellate($W^{\prime\prime},\alpha$)

13: else

14: Sign of $w\leftarrow$ sign of child of $w$ that is not $0$

15: end if

16: end for

17: end for

18:end procedure

When Algorithm 2 is called for the first time, it receives a nested body $G$
created by Algorithm 1 as input, $W=G$. In each subsequent iteration,
Algorithm 2 processes the nested body $W$. It traverses the bodies of $W$ with
undetermined sign, $\pm$, in reverse level order starting from the second
deepest level, assigning a sign $\in\\{-,0,+\\}$ to each. Visiting each
undetermined body $w$ of a level before moving upwards, it derives the sign of
$w$ from the signs of its children, splitting $W$ when required. When finally
the root has been assigned a sign, it is the codomain of a nested mapping.

The sign of an undetermined body $w\in K_{W}$ with height direction $h$ and
children $w^{\downarrow}$ and $w^{\uparrow}$,

$(w,w^{\downarrow},w^{\uparrow})_{W},$ (40)

can be determined from the signs of its children if it fullfills three
conditions. First, both children have a determined sign,

$\text{sign}_{\alpha}(w^{\downarrow})\neq\pm,\quad\text{sign}_{\alpha}(w^{\uparrow})\neq\pm.$
(41)

Second, $\Gamma_{\alpha}\cap w$ is a graph over $\Pi_{h}(\Gamma_{\alpha}\cap
w)$. Third, the projections in height direction $h$ of $w$ and its two
children $w^{\downarrow}$ and $w^{\uparrow}$ coincide,

$\Pi_{h}(w)=\Pi_{h}(w^{\downarrow})=\Pi_{h}(w^{\uparrow}).$ (42)

If the three conditions hold, there are three different cases specifying $w$’s
sign: First, if the body’s children $w^{\downarrow}$ and $w^{\uparrow}$ have
the same sign $s\in\\{-,+\\}$, then $\Gamma_{\alpha}\cap w=\emptyset$ and
$\text{sign}_{\alpha}(w)=s$. Second, if $w^{\downarrow}$ has a sign
$s\in\\{-,+\\}$ and $w^{\uparrow}$ has sign $0$, then $\Gamma_{\alpha}\cap
w=w^{\uparrow}$ and $\text{sign}_{\alpha}(w)=s$, or vice versa. Third, if
$w^{\downarrow}$ and $w^{\uparrow}$ have a different sign $\in\\{-,+\\}$, then
$\Gamma_{\alpha}\cap w$ lies between the children and the sign of $w$ stays
undertermined, $\text{sign}_{\alpha}(w)=\pm$.

To assign a sign in the third case, the algorithm splits $W$ into two nested
bodies, $W\rightarrow W^{\prime},\,W^{\prime\prime}$. It divides $w$ along
$\Gamma_{\alpha}$ into two bodies,
$w={w}^{\prime}\dot{\cup}{w}^{\prime\prime}$, and inserts a new child
$w^{0}=\Gamma_{\alpha}\cap w$ with $\text{sign}_{\alpha}(w^{0})=0$,

$(w,w^{\downarrow},w^{\uparrow})_{W}\rightarrow(w^{\prime},w^{\downarrow},w^{0})_{W^{\prime}},(w^{\prime\prime},w^{0},w^{\uparrow})_{W^{\prime\prime}}.$
(43)

The body-children triple of $w$ is split in two,

$\Pi_{h}(w^{\prime})=\Pi_{h}(w^{\downarrow})=\Pi_{h}(w^{0}),\quad\Pi_{h}(w^{\prime\prime})=\Pi_{h}(w^{0})=\Pi_{h}(w^{\uparrow}).$
(44)

so that the projections of the new triples align (see Figure 6). After the
split, the second case applies to the new bodies ${w}^{\prime}$ and
${w}^{\prime\prime}$. Hence,
$\text{sign}_{\alpha}({w}^{\prime})=\text{sign}_{\alpha}(w^{\downarrow})$ and
$\text{sign}_{\alpha}({w}^{\prime\prime})=\text{sign}_{\alpha}(w^{\uparrow})$.

In order to split $W$, not only $w$ has to be divided. To create the other
bodies of $W^{\prime}$ and $W^{\prime\prime}$, the nodes that are part of a
level higher than $w$’s are copied and added to $W^{\prime}$ and
$W^{\prime\prime}$ in a first step. In a second step, every remaining node $v$
of $K_{W}$ is split, $v=v^{\prime}\dot{\cup}v^{\prime\prime}$, and inserted
into the respective nested body. Each $v$ is split so that the projections in
height direction $h$ of $v$,

$\Pi_{h}(v^{\prime})=\Pi_{h}(v^{\prime\downarrow})=\Pi_{h}(v^{\prime\uparrow}),\quad\Pi_{h}(v^{\prime\prime})=\Pi_{h}(v^{\prime\prime\downarrow})=\Pi_{h}(v^{\prime\prime\uparrow}),$
(45)

of $v^{\prime}$ and its children $v^{\prime\uparrow}$ and
$v^{\prime\downarrow}$, and $v^{\prime\prime}$ and its children
$v^{\prime\prime\uparrow}$ and $v^{\prime\prime\downarrow}$ align. In
particular, this means that each node $v$ of the level of $w$ receives a new
child $v^{0}$ mirroring $w^{0}$,

$(v,v^{\downarrow},v^{\uparrow})_{W}\rightarrow(v^{\prime},v^{\downarrow},v^{0})_{W^{\prime}},(v^{\prime\prime},v^{0},v^{\uparrow})_{W^{\prime\prime}},$
(46)

as illustrated by Figure 6.

Figure 6: In the first step of splitting the nested body $W$ into $W^{\prime}$
and $W^{\prime\prime}$, the body $w^{0}\subseteq\Gamma_{\alpha}$ is inserted
into the body $w$. The resulting bodies $w^{\prime}$ and $w^{\prime\prime}$
have a sign $\in\\{-,+\\}$. In the second step, the siblings of $w$ are split
by inserting the body $v^{0}$, which mirrors $w_{0}$.

As long as an undetermined body is assigned a sign in each recursion and the
three conditions are met by every body-children triple of $W$, the algorithm
will terminate correctly. The three conditions are true when Algorithm 2 is
called for the first time because they hold for each nested body created by
Algorithm 1. Consequently, all three conditions remain true while the
algorithm is iterating and $W$ is not split. When $W$ is split, it is easy to
see that the second and the third condition are true for every body-children
triple of $W^{\prime}$ and $W^{\prime\prime}$; The second condition holds
because the bodies of $W^{\prime}$ and $W^{\prime\prime}$ are subsets of the
bodies of $W$, and the third condition is enforced by construction. However,
the first condition must be actively maintained by the algorithm.

To assert the first condition, each inserted body $v^{0}$ is assigned a sign
after the split. In order to determine the sign, it suffices to evaluate
$\alpha$ at any point $\boldsymbol{x}\in v^{0}$. This is the only time when
explicitly computing a point on the body is required. To compute such a point,
height functions are used.

Each nested body $W$ has $d$ pairs of height functions,

$a_{d-i}^{\downarrow},a_{d-i}^{\uparrow}:\mathbb{R}^{d}\rightarrow\mathbb{R},\quad
i=0,1,...,d-1.$ (47)

Initially, each height function is constant,

$a_{j}^{\downarrow}(\boldsymbol{x})=\min_{\boldsymbol{y}\in
R}(e_{h}\cdot\boldsymbol{y}),\quad
a_{j}^{\uparrow}(\boldsymbol{x})=\max_{\boldsymbol{y}\in
R}(e_{h}\cdot\boldsymbol{y}),\quad h=h_{j},$ (48)

and correlates to a face of the hyperrectangle $R$ in the root of $W$. When a
body $w$ of $W$ is split, the height functions
$a_{j}^{\downarrow}=a^{\downarrow},a_{j}^{\uparrow}=a^{\uparrow}$ of $w$’s
level are adjusted to match the shape of $W^{\prime}$ and $W^{\prime\prime}$
by inserting $a^{0}:\mathbb{R}^{d}\rightarrow\mathbb{R}$,

$(a^{\downarrow},a^{\uparrow})_{W}\rightarrow(a^{\downarrow},a^{0})_{W^{\prime}}\,,(a^{0},a^{\uparrow})_{W^{\prime\prime}}\,.$
(49)

The inserted height function $a^{0}(\boldsymbol{x})$ describes the surface
$\Gamma_{\alpha}\cap w$,

$\Gamma_{\alpha}\cap
w=a^{0}(w),\quad\alpha(\Pi_{h}(\boldsymbol{x})+a^{0}(\boldsymbol{x})e_{h})=0,\quad
a^{\downarrow}(\boldsymbol{x})\leq a^{0}(\boldsymbol{x})\leq
a^{\uparrow}(\boldsymbol{x}),$ (50)

defined by the zero-isocontour of $\alpha$. To evaluate the implicit function,
the Newton method,

$a^{0}_{n+1}=a^{0}_{n}-\frac{\alpha(\boldsymbol{x}_{n})}{\partial_{a^{0}}\alpha(\boldsymbol{x}_{n})}=a^{0}_{n}-\frac{\alpha(\boldsymbol{x}_{n})}{\nabla\alpha(\boldsymbol{x}_{n})\cdot
e_{h}},\quad\boldsymbol{x}_{n}=\Pi_{h}(\boldsymbol{x})+a^{0}_{n}(\boldsymbol{x})e_{h},$
(51)

is used, where the subscript $n$ denotes the iteration. Since $\alpha$ is
monotone in height direction $h$ on $w$, the method generally convergences.
Algorithm 1 provides monotonicity of $\alpha$ since $w$ is a subset of a face
$k$ of $G$,

$\nabla\alpha(x)\cdot e_{h}\neq 0,\quad\forall x\in w\subseteq k,$ (52)

where the gradient’s component $h$ can only assume zero on the boundary of
$k$. However, safeguarding is still required, e.g. restricting step-size.

##  4 Numerical Integration

As described in section 2, we integrate numerically by repeatedly finding
nested mappings that map a unit cube to domains bounded by single level sets.
To generate a set of nested mappings $A_{\Omega}$ for $\Omega_{\alpha}$
bounded by a single level set, Algorithm 1 first divides $K$ into subsets on
which $\Gamma_{\alpha}$ is a graph. Then Algorithm 2 tessellates
$\Omega_{\alpha}$ on each subset, creating the set of nested bodies $V^{-}$.
Finally, the mappings $A_{\Omega}=\\{A_{0},A_{1},...\\}$ are constructed from
the set $V^{-}$ utilizing that each nested body $V_{i}\in V^{-}$ is the
codomain of a nested mapping, $A_{i}:C\rightarrow V_{i}$.

The nested mapping $A\in A_{\Omega}$ can be extracted from the height
functions of its nested body constructed by Algorithm 2. Without loss of
generality, we assume that the height directions are ordered, so that
$h_{1}=x,h_{2}=y,h_{3}=z$. Accordingly, the height functions
$a_{i}^{\downarrow},a_{i}^{\uparrow}$ returned by Algorithm 2 define the
constant, univariate and bivariate height functions,

$\displaystyle a_{x}{\downarrow}=a_{1}^{\downarrow}(0,0,0),\quad
a_{y}^{\downarrow}(x)=a_{2}^{\downarrow}(x,0,0),\quad
a_{z}^{\downarrow}(x,y)=a_{3}^{\downarrow}(x,y,0),$ (53) $\displaystyle
a_{x}^{\uparrow}=a_{1}^{\uparrow}(0,0,0),\quad
a_{y}^{\uparrow}(x)=a_{2}^{\uparrow}(x,0,0),\quad
a_{z}^{\uparrow}(x,y)=a_{3}^{\uparrow}(x,y,0),$ (54)

of $A$,

$A(\tilde{x},\tilde{y},\tilde{z})=\frac{1}{2}\begin{pmatrix}(a_{x}^{\uparrow}-a_{x}^{\downarrow})\tilde{x}+a_{x}^{\uparrow}+a_{x}^{\downarrow}\\\
(a_{y}^{\uparrow}(x)-a_{y}^{\downarrow}(x))\tilde{y}+a_{y}^{\uparrow}(x)+a_{y}^{\downarrow}(x)\\\
(a_{z}^{\uparrow}(x,y)-a_{z}^{\downarrow}(x,y))\tilde{z}+a_{z}^{\uparrow}(x,y)+a_{z}^{\downarrow}(x,y)\end{pmatrix}=\begin{pmatrix}x\\\
y\\\ z\end{pmatrix}.$ (55)

Required in the entries of $A$’s Jacobian matrix, we determine the partial
derivatives of the height functions describing surface $\Gamma_{\alpha}$ by
differentiating its implicit definition,

$\frac{\text{d}}{\text{d}x}\left(\alpha(x,y,a_{z}(x,y))=0\right)\Leftrightarrow\partial_{x}\alpha+\partial_{z}\alpha\,\partial_{x}a_{z}=0,\quad\frac{\text{d}}{\text{d}y}\left(\alpha(x,y,a_{z}(x,y))=0\right)\Leftrightarrow\partial_{y}\alpha+\partial_{z}\alpha\,\partial_{y}a_{z}=0.$
(56)

For brevity, the height function $a_{z}$ denotes $a_{z}^{\uparrow}$ or
$a_{z}^{\downarrow}$. Rearranging yields the partial derivatives of $a_{z}$,

$\partial_{x}a_{z}=-\frac{\partial_{x}\alpha}{\partial_{z}\alpha},\quad\partial_{y}a_{z}=-\frac{\partial_{y}\alpha}{\partial_{z}\alpha},$
(57)

which are fractions of the partial derivatives of $\alpha$. Similarly,
differentiating the implicit surface line,

$\frac{\text{d}}{\text{d}x}\left(\alpha(x,a_{y}(x),0)=0\right)\Leftrightarrow\partial_{x}\alpha+\partial_{y}\alpha\partial_{x}a_{y}=0\Leftrightarrow\partial_{x}a_{y}=-\frac{\partial_{x}\alpha}{\partial_{y}\alpha},$
(58)

determines the partial derivative $\partial_{x}a_{y}$, where again $a_{y}$
denotes $a_{y}^{\uparrow}$ or $a_{y}^{\downarrow}$ for brevity.

###  4.1 Volume Integrals

To integrate over domains bounded by a single level set $\alpha$, the volume
integral,

$\int_{\Omega_{\alpha}}g(\boldsymbol{x})\,\text{d}V=\sum_{A\in
A_{\Omega}}\int_{C}g(A(\tilde{\boldsymbol{x}}))J_{A}(\tilde{\boldsymbol{x}})\,\text{d}V\approx\sum_{A\in
A_{\Omega}}\sum_{i}w_{A,i}g(\boldsymbol{x}_{A,i}),$ (59)

is transformed into a sum of volume integrals over hypercube $C$ by the nested
mappings of $A_{\Omega}$. The resulting volume integrals are evaluated with
transformed quadrature nodes $\boldsymbol{x}_{T,i}$ and weights $w_{T,i}$,

$w_{A,i}=J_{A}(\tilde{x}_{i})\tilde{w}_{i},\quad\boldsymbol{x}_{A,i}=A(\tilde{\boldsymbol{x}}_{i}),$
(60)

derived from the quadrature nodes $\tilde{\boldsymbol{x}}_{i}$ and weights
$\tilde{w}_{i}$ of a quadrature rule on $C$.

Applying the same ansatz twice, the volume integral confined by two level sets
$\alpha$ and $\beta$,

$\int_{\Omega}g(\boldsymbol{x})\,\text{d}V=\sum_{A\in A_{\Omega}}\sum_{B\in
B_{\Omega}}\int_{C}g(T(\hat{\boldsymbol{x}}))J_{T}(\hat{\boldsymbol{x}})\,\text{d}V\approx\sum_{A\in
A_{\Omega}}\sum_{B\in B_{\Omega}}\sum_{i}w_{T,i}g(\boldsymbol{x}_{T,i}),$ (61)

is transformed by the composition of mappings $T$,

$T:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3},\quad T=A\circ B.$ (62)

The codomains of the nested mappings $B_{\Omega}=\\{B_{1},B_{2},...\\}$
tesselate $\Omega_{\tilde{\beta}}$ defined by
$\tilde{\beta}(\tilde{\boldsymbol{x}})=\beta(A(\tilde{\boldsymbol{x}}))$.
Again, the resulting quadrature nodes and weights,

$w_{T,i}=J_{T}(\hat{\boldsymbol{x}}_{i})\hat{w}_{i}=J_{A}(B(\hat{\boldsymbol{x}}_{i}))J_{B}(\hat{\boldsymbol{x}}_{i})\hat{w}_{i},\quad\boldsymbol{x}_{T,i}=T(\hat{\boldsymbol{x}}_{i}).$
(63)

are derived from a quadrature rule on $C$, composed of nodes
$\hat{\boldsymbol{x}}_{i}$ and weights $\hat{w}_{i}$.

Since the mappings $B_{\Omega}$ are bounded by the transformed level set
$\tilde{\beta}$, evaluating the height function $b^{0}$ of $B\in B_{\Omega}$,

$b^{0}:\tilde{\beta}(\Pi_{h}(\tilde{\boldsymbol{x}})+b^{0}(\tilde{\boldsymbol{x}})\boldsymbol{e}_{h})=0$
(64)

defined by the level set’s zero-isocontour requires the Jacobian matrix $DA$.
The Jacobian matrix is required in the newton method

$\quad
b^{0}_{n+1}=b^{0}_{n}-\frac{\tilde{\beta}(\tilde{\boldsymbol{x}}_{n})}{\partial_{b^{0}}\tilde{\beta}(\tilde{\boldsymbol{x}}_{n})}=b^{0}_{n}-\frac{\beta(A(\tilde{\boldsymbol{x}}_{n}))}{\nabla\beta(A(\tilde{\boldsymbol{x}}_{n}))\,DA\cdot\boldsymbol{e}_{h}},\quad\boldsymbol{x}_{n}=\Pi_{h}(\tilde{\boldsymbol{x}})+b^{0}_{n}(\tilde{\boldsymbol{x}})\boldsymbol{e}_{h},$
(65)

where the composition $\tilde{\beta}$ is differentiated.

###  4.2 Surface Integrals

Let the set $A_{\Gamma}\subseteq A_{\Omega}$ group all mappings with a
codomain $V_{i}$ adjacent to the surface $\Gamma_{\alpha}$, so that
$\Gamma_{\alpha}\cap V_{i}\neq\emptyset$. Moreover, let the mapping $E_{A}$
embed a unit square in the face of the domain of the mapping $A\in A_{\Gamma}$
that maps to $\Gamma_{\alpha}$. More precisely, if we assume without loss of
generality that $\Gamma_{\alpha}$ forms the upper face of the codomain of $A$
and $h_{3}=z$,

$A(C^{\uparrow})\subseteq\Gamma_{\alpha},\quad
C^{\uparrow}=\text{face}_{z}^{\uparrow}(C),$ (66)

then $E_{A}$ is the function,

$E_{A}:C^{2}\rightarrow C^{\uparrow},\quad E_{A}(x,y)=(x,y,1)^{T},$ (67)

that embedds a unit square $C^{2}$ on the upper face of $C$.

Facilitated by $E_{A}$, the composition of mappings $T$,

$T:\mathbb{R}^{2}\rightarrow\mathbb{R}^{3},\quad T=A\circ E_{A},$ (68)

maps a two-dimensional unit square onto $\Gamma_{\alpha}$. For each mapping
$A\in A_{\Gamma}$, the composition transforms the unit square to the surface
integral,

$\oint_{\Gamma_{\alpha}}g(\boldsymbol{x})\,\text{d}S=\sum_{A\in
A_{\Gamma}}\int_{C^{2}}g(T(\tilde{\boldsymbol{x}}))\,G_{T}(\tilde{\boldsymbol{x}})\,\text{d}S\approx\sum_{A\in
A_{\Gamma}}\sum_{i}w_{T,i}g(\boldsymbol{x}_{T,i}).$ (69)

Subsequently, the surface quadrature nodes $\boldsymbol{x}_{T,i}$ and weights
$w_{T,i}$ are derived,

$w_{T,i}=G_{T}(\tilde{\boldsymbol{x}}_{i})\tilde{w}_{i},\quad\boldsymbol{x}_{T,i}=T(\tilde{\boldsymbol{x}}_{i}),$
(70)

from a set of quadrature nodes $\tilde{\boldsymbol{x}}_{i}$ and weights
$\tilde{w}_{i}$ on $C^{2}$. Since $T$ is an embedding, the Gram determinant
which generalizes the Jacobian determinant,

$G_{T}(\tilde{\boldsymbol{x}})=\sqrt{\det\left(DT(\tilde{\boldsymbol{x}})^{T}\,DT(\tilde{\boldsymbol{x}})\right)},$
(71)

encodes the geometry of the surface.

The integral on a surface of a domain confined by two level sets can be found
by injecting the embedding $E_{A}$ into the composition of mappings $A\circ B$
of the volume integral. We compose $T$,

$T:\mathbb{R}^{2}\rightarrow\mathbb{R}^{3},T=A\circ E_{A}\circ B,$ (72)

and transform $\gamma_{\alpha}$ to a set of unit squares,

$\oint_{\gamma_{\alpha}}g(\boldsymbol{x})\,\text{d}S=\sum_{A\in
A_{\Gamma}}\sum_{B\in
B_{\Omega}}\int_{C^{2}}g(T(\hat{\boldsymbol{x}}))\,G_{T}(\hat{\boldsymbol{x}})\,\text{d}S\approx\sum_{A\in
A_{\Gamma}}\sum_{B\in B_{\Omega}}\sum_{i}w_{T,i}g(\boldsymbol{x}_{T,i}),$ (73)

to construct a quadrature rule over $\gamma_{\alpha}$. The composition $T$
transforms the integral in two steps. First, the mapping $A\circ E_{A}$ maps a
unit square onto a subset of the surface $\Gamma_{\alpha}$. Then, after
transforming $\tilde{\beta}=\beta(A\circ E_{A})$ to this unit square, the
transformed integral is treated as a two-dimensional volume integral.
Consequently, the codomains of the mappings $B_{\Omega}$ tesselate the two-
dimensional domain $\Omega_{\tilde{\beta}}$.

Similiarly, a quadrature rule over $\gamma_{\beta}$ is provided by the
composition $T$,

$T:\mathbb{R}^{2}\rightarrow\mathbb{R}^{3},T=A\circ B\circ E_{B}.$ (74)

The composition restricts the domain of integration to $\Omega_{\alpha}$ and
then maps the surface defined by $\tilde{\beta}$ to $C^{2}$,

$\oint_{\gamma_{\beta}}g(\boldsymbol{x})\,\text{d}S=\sum_{A\in
A_{\Omega}}\sum_{B\in
B_{\Gamma}}\int_{C^{2}}g(T(\tilde{\boldsymbol{x}}))\,G_{T}(\tilde{\boldsymbol{x}})\,\text{d}S\approx\sum_{A\in
A_{\Omega}}\sum_{B\in B_{\Gamma}}\sum_{i}w_{T,i}g(\boldsymbol{x}_{T,i}),$ (75)

using the mappings $E_{B}$ of each mapping $B\in B_{\Gamma}$. Corresponding to
the definition of $A_{\Gamma}$, the set $B_{\Gamma}\subseteq B_{\Omega}$
groups the mappings with a codomain in $\Omega_{\tilde{\beta}}$ that is
adjacent to $\tilde{\beta}$.

###  4.3 Line Integrals

Applying the composition of transformations $T$,

$T:\mathbb{R}\rightarrow\mathbb{R}^{3},T=A\circ E_{A}\circ B\circ E_{B},$ (76)

yields the line quadrature rule,

$\fint_{\partial\gamma}g(\boldsymbol{x})\,\text{d}L=\sum_{A\in
A_{\Gamma}}\sum_{B\in
B_{\Gamma}}\int_{C^{1}}g(T(\hat{\boldsymbol{x}}))\,G_{T}(\hat{\boldsymbol{x}})\,\text{d}L\approx\sum_{A\in
A_{\Gamma}}\sum_{B\in B_{\Gamma}}\sum_{i}w_{T,i}g(\boldsymbol{x}_{T,i}).$ (77)

The composition transforms the integral in two steps. First, the mapping
$A\circ E_{A}$ maps a unit square $C^{2}$ onto the surface $\Gamma_{\alpha}$.
Then, $B\circ E_{B}$ maps a unit line $C^{1}$ to the surface defined by
$\tilde{\beta}=\beta(A\circ E_{A})$ in $C^{2}$.

The line quadrature nodes $\boldsymbol{x}_{T,i}$ and weights $w_{T,i}$,

$w_{T,i}=G_{T}(\tilde{\boldsymbol{x}}_{i})\tilde{w}_{i},\quad\boldsymbol{x}_{T,i}=T(\tilde{\boldsymbol{x}}_{i}),$
(78)

are derived from a one-dimensional set of quadrature nodes
$\tilde{\boldsymbol{x}}_{i}$ and weights $\tilde{w}_{i}$ on $C^{1}$.

##  5 Experiments

We conduct five numerical experiments to investigate the method[35]. In the
first experiment, we evaluate the edge length, surface area and volume of a
domain with a sharp, oscillating edge. Offering a curved topology with a
Gaussian curvature of zero, the domain’s shape is defined by trigonometric
level sets. In the second experiment we investigate a shape with non-zero
Gaussian curvature and evaluate the edge length, surface area and volume of a
lens defined by two intersecting spheres. Examining both an intricate shape
and integrand in the third experiment, we evaluate the line and surface
integral of a trigonometric function on a curved toric section. In the last
two experiments, we investigate shapes described by a single level set that
require a great amount of quadrature nodes due to subdivisions. We reduce the
amount of quadrature nodes by inserting a second level set.

In each experiment we measure the error $e$,

$e=\left|G-\sum_{i}^{c^{3}}\sum_{T\in
T_{i}}\sum_{j}^{n^{d}}w_{T,j}g\left(\boldsymbol{x}_{T,j}\right)\right|\quad
G=\fint_{\partial\gamma}g(\boldsymbol{x})\,\text{d}L,\oint_{\gamma_{\beta}}g(\boldsymbol{x})\,\text{d}S,\int_{\Omega}g(\boldsymbol{x})\,\text{d}V,$
(79)

which is the difference between the exact integral and its approximation. For
the sake of brevity, $e$ denotes the error for every domain of integration,
regardless if its a volume ($d=3$) a surface ($d=2$) or a line ($d=1$).
Sharing a basic setting, each experiment is conducted in a cube $K=[-1,1]^{3}$
which is subdivided into $c^{3}$ uniform cells. To compute the quadrature
rules for each cell, we create a set of mappings $T_{i}$ for each $i$-th cell.
Subsequently each mapping $T$ of $T_{i}$ is applied to transform a tensorized
Gauss quadrature rule constisting of $n^{d}$ quadrature nodes, constructing
the transformed quadrature nodes, $\boldsymbol{x}_{T,j}$, and weights,
$w_{T,j}$.

###  5.1 Oscillating Edge

Let the level sets $\alpha$ and $\beta$,

$\displaystyle\alpha(x,y,z)=z-\frac{1}{5}\sin\left(\frac{20\pi x}{11}\right),$
(80) $\displaystyle\beta(x,y,z)=y-\frac{1}{5}\sin\left(\frac{20\pi
x}{11}\right),$ (81)

define a trigonometric domain in the cube $K$. The intersection of the level
sets creates a sharp edge, oscillating around the $x$-axis. Focusing on the
geometry of the domain of integration, we integrate a constant function,

$g(x,y,z)=1,$ (83)

and evaluate the length of $\partial\gamma$, the surface area of
$\gamma_{\beta}$ and the volume of $\Omega$. Choosing a constant $g$ reduces
the integrand of each transformed integral to the Jacobian determinant and
removes any influence of $g$.

To compute the exact solution, we precisely evaluate the integrals with an
established quadrature method [33],

$\displaystyle\fint_{\partial\gamma}1\,\text{d}L=2.9018098242473137628716230441128...\,,$
(84)
$\displaystyle\oint_{\gamma_{\beta}}1\,\text{d}S=2.5048230500093248969863804012397...\,,$
(85)
$\displaystyle\int_{\Omega}1\,\text{d}V=2.0431849934260147426243415665995...\,,$
(86)

using the analytic expressions of $\gamma_{\alpha}$ and $\gamma_{\beta}$.

Figure 7: Two trigonometric surfaces $\gamma_{\alpha}$ and $\gamma_{\beta}$
that form an oscillating edge at their intersection. Each filled circle
denotes a node of a transformed Gauss quadrature rule, $n=3$. The radius of
each node is proportional to its weight. The left picture shows the quadrature
rule obtained for the oscillating edge on a grid consisting of $c^{3}=8^{3}$
uniform cells (a). The transformed quadrature rules on the surfaces (b) and in
the volume (c) are depicted for a single cell $k$ on the right.

Figure 7 depicts a section of the domains of integration
$\partial\gamma,\gamma_{\alpha},\gamma_{\beta}$ and $\Omega$ in proximity of
$\partial\gamma$ and the corresponding quadrature rules for $n=3$ and $c=8$.
The nodes of the line quadrature rule of $\partial\gamma$ are located on the
edge, showing the characteristic spacing of the underlying Gauss rule. Groups
of $n$ nodes are positioned with varying distance, depending on the length of
the intersection of the edge and the cell. By construction, each intersection
contains at least $n$ nodes. Similarly, the nodes of the surface quadrature
rules of $\gamma_{\alpha},\gamma_{\beta}$ are located on the surface and the
nodes of the volume quadrature rule of $\Omega$ are located in the volume.

Figure 8: The error $e$ obtained during the evaluation of the integrals of the
oscillating edge. Each curve is associated to $n$ underlying Gauss nodes. When
the number of cells $c$ per coordinate axis is increased, the error converges
with the order depicted in the bottom left of each graph.

If $c$ is successively increased, the observed convergence order is in good
accordance with the expected convergence order of $2n$, Figure 8 shows. It
shows the error $e$ and the order of convergence obtained while evaluating the
line, surface and volume integrals.

###  5.2 Spherical Lens

Creating a quarter lens with a sharp edge, let the zero-isocontours of
$\alpha$ and $\beta$ describe two intersecting spheres,

$\displaystyle\alpha(x,y,z)=(x+1)^{2}+(y+1)^{2}+\left(z+\frac{49}{100}\right)^{2}-\left(\frac{9}{10}\right)^{2},$
(87)
$\displaystyle\beta(x,y,z)=(x+1)^{2}+(y+1)^{2}+\left(z-\frac{51}{100}\right)^{2}-\left(\frac{9}{10}\right)^{2},$
(88) $\displaystyle K=[-1,1]^{3},$ (89)

in cube $K$.

Focusing on the geometry of the domain of integration again, we evaluate the
length of $\partial\gamma$, the surface area of $\gamma_{\beta}$ and the
volume of $\Omega$ and choose $g$ accordingly,

$g(x,y,z)=1.$ (90)

The exact edge length, surface area and volume,

$\displaystyle\fint_{\partial\gamma}1\,\text{d}L=\frac{\sqrt{14}}{10}\pi,$
(91) $\displaystyle\oint_{\gamma_{\beta}}1\,\text{d}S=\frac{9}{50}\pi,$ (92)
$\displaystyle\int_{\Omega}1\,\text{d}V=\frac{23}{375}\pi,$ (93)

can be found analytically.

Figure 9: A lens created by intersecting spheres, defined by level sets
$\alpha$ and $\beta$ (a). The red lines on the surface $\gamma_{\beta}$
outline the codomains of each surface mapping created on a grid consisting of
$c^{3}=10^{3}$ uniform cells. The right side (b) shows a detailed view of the
Jacobian determinant, $\det T$, of a surface mapping belonging to a surface
patch bounded by $\partial\gamma$ located in cell $k$.

Showing the surfaces $\gamma_{\alpha}$ and $\gamma_{\beta}$ of the lens,
Figure 9 illustrates two properties of the surface mappings constructed on a
grid with $c^{3}=10^{3}$ cells. First, it outlines the codomains of the
surface mappings of $\gamma_{\alpha}$ of each cell with red lines. The
codomains tesselate $\gamma_{\alpha}$ and vary in shape, depending on the
relative position of the cells and the surface. Second, the figure illustrates
how the geometry of a surface patch in a cell containing a part of the edge is
encoded in the Jacobian determinant. Encoding the transformation to the unit
square, the Jacobian determinant is non-polynomial.

Figure 10: Surface area error of the spherical lens under grid refinement and
grid refinement in combination with adaptive quadrature. A linear trend is
indicated by a dotted line for each transformed Gauss rule $n$. The right
panel shows the total amount of quadrature nodes of each rule. The dotted
lines depict the amount of nodes created while refining the grid and the solid
lines depict the amount of nodes created by additionally applying adaptive
quadrature.

When the grid is refined, the convergence order of the surface error is as
expected for $c$ below approximately 32, the left panel of Figure 10 shows.
After $c$ has reached approximately 32 however, the convergence order
deteriorates and for some values of $c$ the error is significantly higher than
the trend line. This effect becomes more pronounced when $n$ is increased. The
fact that the error decreases when $n$ is increased and the value of $c$ is
held constant, suggests that the irregular error is caused by the Jacobian
determinant $\det T$ in the integrand. Depending on its surrounding cell, the
Jacobian determinant $\det T$ of a mapping may be a ill-behaved non-polynomial
function. When this integrand is badly approximated by the polynomial Gauss
quadrature rule, the error in these cells converges with increasing $n$, but
lies over the average error level nonetheless.

To recover accuracy, we apply adaptive quadrature in the domains of the
mappings. Since the domain of each mapping is a $i$-dimensional hypercube $C$,
many methods can be applied, such as Gauss-Kronrod quadrature or methods
specially tailored to rational functions. For simplicity, we apply a
straightforward adaptive quadrature method. The quadrature rule $Q$ belonging
to the mapping is constructed by placing tensorized Gauss nodes and weights in
$C$, so that $|Q|=n^{i}$. If the estimated error $e_{Q}$ of the quadrature
rule is greater than a threshold $\tau$,

$e_{Q}>\tau,$ (94)

$C$ is subdivided into $2^{i}$ hypercubes of equal size. Repeating the
process, the quadrature rule $\tilde{Q}$ belonging to the mapping is
constructed by placing tensorized Gauss nodes and weights in every new
hypercube, so that $|\tilde{Q}|=2^{i}n^{i}$. If the estimated error in one of
the new hypercubes is too great still, it is subdivided again. And so on until
a predetermined depth is reached.

Error estimation is a critical component of an adaptive quadrature method and
many error erstimation schemes have been presented [34]. We choose a local
relative error,

$e_{Q}=\left|\frac{\sum_{i}g\left(\boldsymbol{x}_{Q,i}\right)w_{Q,i}-\sum_{i}g\left(\tilde{\boldsymbol{x}}_{\tilde{Q},i}\right)\tilde{w}_{\tilde{Q},i}}{\sum_{i}g\left(\boldsymbol{x}_{Q,i}\right)w_{Q,i}}\right|$
(95)

relating the result of quadrature rule $Q$ to the result of the subdivided
quadrature rule $\tilde{Q}$.

Figure 11:  Line length error and volume error of the spherical lens under
grid refinement and grid refinement in combination with adaptive quadrature. A
linear trend is indicated by a dotted line for each transformed Gauss rule
$n$. The right panel shows the total amount of quadrature nodes of each rule.
The dotted lines depict the amount of nodes created while refining the grid
and the solid lines depict the the amount of nodes created by additionally
applying adaptive quadrature.

The middle panel of Figure 10 depicts the effect of adaptive quadrature on the
surface integrals of the lens with the thresholds $\tau_{n}$ of each $n$. It
lists the error for the values $\tau_{1}=0.1,\tau_{2}=0.001,\tau_{3}=1\times
10^{-5},\tau_{4}=1\times 10^{-7}$, showing that accuracy can be recovered with
the adapative quadrature method. However, it decreases the error of lower $c$
to greater extent, flattening its trend and leading to a convergence order
which is lower than $2n$. This is also reflected by the total number of
quadrature nodes created which is shown for each $n$ in the right panel of
Figure 10. For lower values of $c$, the number of quadrature nodes of the
adaptive quadrature method lies above the number of quadrature nodes of the
method without adaptation. For higher values of $c$, the number of quadrature
nodes created adaptively slightly lies above the number of nodes created
without adaptation in the cases where the error of the method without
adaptation produces a disproportionate error. Otherwise they are almost equal.

Deviation from the expected error at certain resolutions $c$ can also be
observed in the computation of the edge length and the volume, illustrated by
Figure 11. In each case, adaptive quadrature is applied as a remedy,
effectively removing unwanted deviation of the error. The thresholds
$\tau_{n}$ of each $n$ used for adaptive quadrature of the line are,
$\tau_{1}=0.002,\tau_{2}=1\times 10^{-5},\tau_{3}=1\times
10^{-8},\tau_{4}=1\times 10^{-11}$. The respective thresholds used for the
volume are, $\tau_{1}=1,\tau_{2}=0.05,\tau_{3}=0.001,\tau_{4}=1\times
10^{-4}$. When choosing the threshold $\tau$ carefully, the amount of
quadrature nodes increases only slightly, the right panel of Figure 10 and
Figure 11 shows.

Figure 12: The error $e$ obtained during the evaluation of the integrals of
the spherical lens for different numbers of subdivisions $s$ on a grid of
$c=10$ cells per coordinate axis. To increase accuracy, $s^{i}$ quadrature
rules $n$ are placed in each $i$-dimensional domain of each mapping. The
number in the bottom left of each graph denotes the order of convergence.

The previous experiments show that the error decreases when increasing $n$ and
holding $c$ constant. Motivated by this, we choose a grid, $c=10$, and
directly refine the quadrature rules in each domain $C$ of the grid’s
mappings. Each $i$-dimensional domain $C$ is subdivided evenly into $s^{i}$
hypercubes and a tensorized Gauss rule $n$ placed in each. Providing the trend
of the error when successively increasing $s$, Figure 12 suggests an order of
convergence of about $2n$ for each curve $n$.

###  5.3 Toric Section

Respectively describing a torus and a tilted, wavy surface, the level sets,

$\displaystyle\alpha(x,y,z)=\left(\sqrt{x^{2}+y^{2}}-\frac{3}{5}\right)^{2}+z^{2}-\left(\frac{3}{10}\right)^{2},$
(96)
$\displaystyle\beta(x,y,z)=-\frac{1}{2}y+\frac{7}{8}x-\frac{1}{5}\sin(\frac{20\pi}{11}x)\sin\left(\frac{35\pi}{22}y+\frac{10\pi}{11}z\right),$
(97)

generate a curved surface $\gamma_{\beta}$ that is skewed to the grid (see
Figure 13). The function $g$ is a combination of two trigonometric functions,

$g(x,y,z)=\frac{1}{5}\sin(\frac{20\pi}{11}x)\sin\left(\frac{35\pi}{22}y+\frac{10\pi}{11}z\right).$
(98)

Since $\gamma_{\beta}$ and $g$ are line-symmetric with respect to the
$x$-axis, the integrals of $g$ along $\partial\gamma$ and over
$\gamma_{\beta}$,

$\displaystyle\fint_{\partial\gamma}g(\boldsymbol{x})\,\text{d}L=0,$ (99)
$\displaystyle\oint_{\gamma_{\beta}}g(\boldsymbol{x})\,\text{d}S=0,$ (100)

vanish.

Figure 13:  A tilted wavy surface intersects a torus, generating a curved
toric section (a). The red lines on the surface $\gamma_{\beta}$ outline the
codomains of each surface mapping created on a grid consisting of
$c^{3}=30^{3}$ uniform cells. The right (b) shows a line quadrature rule on
the the boundary line $\partial\gamma$ with $n=1$. The projections of the line
onto the $xy$-, $xz$\- and $yz$-plane are depicted in grey.

The curved intersection $\partial\gamma$ of the torus and the wavy sheet is
depicted in the right picture of Figure 13. Its intricate shape is outlined by
the quadrature nodes of a line quadrature rule. The line quadrature rule is
immersed in a grid of $c^{3}=30^{3}$ cells and consists of $n=1$ Gauss nodes
per nested mapping. Although the shape of $\partial\gamma$ is complicated, the
nodes are placed evenly along the boundary. From the projections of
$\partial\gamma$ onto the $xy$-, $xz$\- and $yz$-plane Figure 13 it can be
seen, that the boundary is line-symmetric.

Figure 14: The error $e$ obtained during the evaluation of the line and
surface integrals of the toric section with $n$ underlying Gauss nodes. The
set of mappings created on a grid of $c^{3}=30^{3}$ cells is kept constant. To
increase accuracy, $s^{i}$ quadrature rules $n$ are placed in each
$i$-dimensional domain of each mapping. The number in the bottom left of each
graph denotes the order of convergence.

The function $g$ is a composition of trigonometric functions that creates a
wavy, smooth field throughout $K$. We evaluate the line and surface integral
of $g$ with a quadrature rule with $n$ nodes on a grid of $c^{3}=30^{3}$
cells. When subdividing the quadrature rules direclty in each domain of the
grid’s mappings, we measure a convergence rate of the error of about $2n$, as
illustrated by Figure 14. This shows that high order convergence can be
achieved for intricate geometries and functions.

###  5.4 Thin Cylindrical Sheet

Depending on the relative position of the zero-isocontour and $K$, the
construction of an exact nested mapping may require a high number of
subdivisions, even for a quadratic level set. For example, let the surface of
a quadratic level set $\alpha$,

$\displaystyle\alpha(x,y,z)=\tilde{\alpha}(x,y,z)^{2}-l^{2}\quad\tilde{\alpha}(x,y,z)=\left(\frac{4}{5}x-y\right)$
(101)

be two parallel planar surfaces separated by distance $l$, creating a slanted
sheet. Both components of the linear gradient of $\alpha$,

$\nabla\alpha=2\tilde{\alpha}\nabla\tilde{\alpha}=2\tilde{\alpha}\begin{pmatrix}\frac{4}{5}\\\
-1\\\ 0\end{pmatrix}$ (102)

contain the zero-isocontour $\Gamma_{\tilde{\alpha}}$,

$\Gamma_{\tilde{\alpha}}=\left\\{\boldsymbol{x}\in
K\,|\,\tilde{\alpha}(\boldsymbol{x})=0\right\\}$ (103)

which is pinched between the two planar surfaces in $\Gamma_{\alpha}$.

$l$ | Subdivision | Contour Splitting
---|---|---
0.1 | 2426 | 2
0.03 | 29010 | 2
0.01 | 248162 | 2
Amount of volume quadrature nodes

Figure 15: Slanted sheet created by two parallel planar surfaces separated by
distance $l$ of a single level set $\alpha$. Pinched between the surfaces, the
zero-isocontour of both gradients $\Gamma_{\tilde{\alpha}}$ is indicated by a
dashed red line. The subsets of $K$ created by Algorithm 1 with unrestricted
subdivision are outlined in the $xy$-Plane. Comparing the subdivision and
contour splitting approach, the table on the right lists the number of nodes
contained in volume quadrature rules, $n=1$, for different values of $l$.

When Algorithm 1 tries to represent the surfaces by a graph, it will select a
height direction $h$ and subdivide $K$ until the zero-isocontour of
$\partial_{h}\alpha$ exclusively lies in subsets of $K$ not containing
$\gamma_{\alpha}$ (see Figure 15). Since the the zero-isocontour of
$\partial_{h}\alpha$, $\Gamma_{\tilde{\alpha}}$, is pinched between the two
surfaces of $\Gamma_{\alpha}$, decreasing distance $l$ causes a substantial
increase of subdivisions during the construction of quadrature rules. If
however $\Omega_{\alpha}$ is split along $\Gamma_{\tilde{\alpha}}$ before
creating the quadrature rule, subdivisions are not required and the number of
nodes constructed is reduced greatly. This is demonstrated by the table on the
right side of Figure 15 which shows the number of nodes created for a volume
integral for different values of $l$. In the case of $l=0.01$, splitting along
the zero-isocontour leads to approximately $62^{3}$ fewer nodes.

Reducing the amount of quadrature nodes is realized in two steps. First, we
split $\Omega_{\alpha}$ along the zero-isocontour of $\partial_{h}\alpha$ into
$\Omega^{\prime}_{\alpha}$ and $\Omega^{\prime\prime}_{\alpha}$,

$\Omega_{\alpha}=\Omega^{\prime}_{\alpha}\dot{\cup}\Omega^{\prime\prime}_{\alpha},$
(104)

after the algorithm has chosen a height direction $h$. The two new domains
result from intersecting $\Omega_{\alpha}$ with the domains $\Omega_{\beta}$
and $\Omega_{\beta^{-}}$,

$\Omega^{\prime}_{\alpha}=\Omega_{\alpha}\cap\Omega_{\beta},\quad\Omega^{\prime\prime}_{\alpha}=\Omega_{\alpha}\cap\Omega_{\beta^{-}},\quad\beta=\partial_{h}\alpha,\quad\beta^{-}=-\partial_{h}\alpha.$
(105)

defined by the level sets $\beta$ and $\beta^{-}$. Second, we construct
quadrature rules on the two new domains, applying the approach for two level
sets. Because the zero-isocontour of $\partial_{h}\alpha$ lies on their
boundary, subdivisions are not required on the new domains.

Figure 16: A thin cylindrical sheet with two concentric surfaces
$\Gamma_{\tilde{\alpha}-l}$ and $\Gamma_{\tilde{\alpha}+l}$, encasing a zero-
isocontour $\Gamma_{\tilde{\alpha}}$ of the level set’s gradient drawn as a
red, transparent surface. A section of the surface quadrature nodes, $n=1$, of
both surfaces is displayed in the middle. The quadrature nodes are created on
a grid containing $c^{3}=20^{3}$ cells, indicated by grey lines. A convergence
plot of the error of the surface area is shown on the right, where $s^{i}$
quadrature rules $n$ are placed in each $i$-dimensional domain of each
mapping.

This ansatz can be applied generally to domains of smooth levels sets with
derivatives unfavorable for subdivision. For example, let $\alpha$ be a thin
cylindrical sheet,

$\alpha(x,y,z)=\tilde{\alpha}(x,y,z)^{2}-l^{2}\quad\tilde{\alpha}(x,y,z)=x^{2}+y^{2}+\left(z+\frac{6}{5}\right)^{-2}\quad
l=\frac{1}{50}$ (106)

that has two concentric surfaces $\Gamma_{\tilde{\alpha}-l}$ and
$\Gamma_{\tilde{\alpha}+l}$,

$\Gamma_{\tilde{\alpha}-l}=\\{\boldsymbol{x}\in
K\,|\,\tilde{\alpha}(\boldsymbol{x})-l=0\\},\quad\Gamma_{\tilde{\alpha}+l}=\\{\boldsymbol{x}\in
K\,|\,\tilde{\alpha}(\boldsymbol{x})+l=0\\},$ (107)

that are curved and cylindrical. The gradient of $\alpha$ contains the zero-
isocontour $\Gamma_{\tilde{\alpha}}$ of the curved cylinder $\tilde{\alpha}$,

$\nabla\alpha(x,y,z)=2\tilde{\alpha}\,\nabla\tilde{\alpha},$ (108)

which lies between $\Gamma_{\tilde{\alpha}-l}$ and $\Gamma_{\tilde{\alpha}+l}$
(see Figure 16).

By first splitting the domain along $\Gamma_{\tilde{\alpha}}$, quadrature
rules for the thin sheet defined by a single level set can be found without
excessive subdivision. The middle panel of Figure 16 shows the quadrature
nodes of both surfaces of a quadrature rule with $n=1$. They are obtained in a
grid consisting of $c^{3}=20^{3}$ cells. The amount of quadrature nodes
depicted in the figure is approximately proportional to the number of cells
intersecting the surfaces. This reflects a low number of subdivisions,
although both surfaces are close and intersect the same cells.

In order to measure the order of convergence, we evaluate the surface areas of
$\Gamma_{\tilde{\alpha}-l}$ and $\Gamma_{\tilde{\alpha}+l}$,

$g(x,y,z)=1,$ (109)

by applying the splitting ansatz. The exact surface areas can be found
analytically,

$\displaystyle\oint_{\Gamma_{\tilde{\alpha}}-l}1\,\text{d}S=6.775163182554902237379363684639...,$
(110)
$\displaystyle\oint_{\Gamma_{\tilde{\alpha}}+l}1\,\text{d}S=6.2377886792891965132267781181652...,$
(111)

using that the surfaces are a revolution arround the $z$-axis. We create the
nested mappings required to create the quadrature rule in a grid consisting of
$c^{3}=20^{3}$ cells and successively subdivide the domains of each nested
mapping. The resulting error can be seen in the right panel of Figure 16.
Despite the complex geometry of the integral’s domain, the error obtained for
$n=1,2,3,4$ converges exponentially with an order of approximately $1.5n$.

###  5.5 Wavy Cylinder

If evaluation of $g$ is computationally expensive, constructing a minimal
amount of quadrature nodes is advantageous. Unfortunately, not only thin
sheets immersed in a coarse grid can cause a high number of subdivisions
leading to many quadrature nodes, like it is the case in subsection 5.4. Even
well resolved bodies can cause a high number of subdivisions when exact
mappings are required.

Let $\alpha$ describe a wavy cylinder,

$\alpha(x,y,z)=x^{2}+\left(y+\frac{1}{5}\sin(\pi
z)\right)^{2}-\left(\frac{2}{3}+l\right)^{2},$ (112)

that is immersed in a grid of $c^{3}=6^{3}$ cells. A small perturbation $l$ of
the radius leads to a high number of quadrature nodes due to subdivision in
the cells located at the intersection of the zero-isocontour
$\Gamma_{\tilde{\alpha}}$ of $\partial_{y}\alpha$,

$\partial_{y}\alpha(x,y,z)=2\left(y+\frac{1}{5}\sin(\pi z)\right),$ (113)

and $\Gamma_{\alpha}$, Figure 17 shows for $l=0.001$.

As described in subsection 5.4, this is caused when $\Gamma_{\tilde{\alpha}}$
is pinched between two subsets of $\Gamma_{\alpha}$. For example, let us look
at the three cells $k$ that are located at the intersection of
$\Gamma_{\tilde{\alpha}}$ and $\Gamma_{\alpha}$ depicted in Figure 17. Due to
the slight shift of the wavy cylinder in $x$-direction, $k$ contains a thin
region of $\Omega_{\alpha}$. On the intersection of $\Omega_{\alpha}$ and the
face $p^{\downarrow}=\text{face}^{\downarrow}_{x}(p)$ of the cell $p\in k$,
the line $\Gamma_{\tilde{\alpha}}\cap p^{\downarrow}$ is pinched between two
close surface sections $\Gamma_{\alpha}\cap p^{\downarrow}$. When Algorithm 1
constructs the graph on $p$, it successively subdivides $p$ until the face
$p^{\prime\downarrow}$ of each subcell $p^{\prime}$ intersects only one of the
three lines,

$p^{\prime\downarrow}=\text{face}^{\downarrow}_{x}(p^{\prime}),\quad
p^{\prime\downarrow}\cap\Gamma_{\tilde{\alpha}}=\emptyset\,\lor\,p^{\prime\downarrow}\cap\Gamma_{\alpha}=\emptyset,$
(114)

resulting in a high number of subdivisions.

Figure 17: A wavy cylinder in a grid consisting of $c^{3}=6^{3}$ cells that is
slightly shifted. Each cell of the grid is indicated by grey lines. The two
panels display the surface quadrature nodes of a subset of cells $k$ of the
grid. The surface quadrature nodes are constructed by (a) using subdivision
and (b) splitting the cells along the zero-isocontour of $\partial_{y}\alpha$
which is indicated by a red, transparent surface.

To decrease the amount of subdivisions required to find exact mappings, we
split the cells along the zero-isocontour of $\partial_{y}\alpha$ like in
subsection 5.4. Figure 17 shows the surface quadrature nodes in the three
cells $k$ for both the approach with subdividing and the approach with contour
splitting. The two sets of surface quadrature nodes displayed in the figure
are constructed for the same distance $l=0.001$ with $n=1$ and $n=2$. By
comparing the surface quadrature nodes constructed by both approaches, it can
be seen that contour splitting results in substantially fewer surface
quadrature nodes for both values of $n$. It is clear from the results of
subsection 5.4, that the difference in the amount of quadrature nodes of the
two approaches further increases, when $l$ is decreased.

##  6 Summary and Outlook

In this paper we present a high-order method that provides numerical
integration of smooth functions on volumes, surfaces, and lines defined
implicitly by two intersecting level sets. Our approach allows for general
smooth level sets and can map arbitrary quadrature rules defined on hypercubes
to the curved domains of the integrals. Conceptually, the mappings required to
create the quadrature rules are constructed by composing two nested mappings.
The first nested mapping maps a hypercube to the domain defined by the first
level set. The second nested mapping maps a hypercube to the domain defined by
the second level set function composed with the first nested mapping which
mirrors the domain of the integral. Algorithmically, each nested mapping is
generated by recasting the zero-isocontours of the level sets as the graphs of
height functions. This is achieved by augmenting a method developed for the
integration on domains bounded by single level sets[29].

The practicality and usefullness of this approach is demonstrated in various
numerical experiments. We apply the approach to compute the volume, the
surface area, and the intersection line length of an oscillating edge formed
by two trigonometric level sets and show high-order convergence of the error.
In the investigation of the same properties of a spherical lens defined by two
quadratic level sets, we measure high-order convergence when the quadrature
rules are refined on a constant grid. Fluctuating error levels introduced by
the relative position of the grid and the domain of integration are
successfully removed by applying an adaptive quadrature scheme, recovering
high accuracy and showing that adapative quadrature conveniently suits the
method. Finally, to demonstrate high-order convergence of the error of an
integral on a complex non-polynomial geometry with an intricate integrand, a
trigonometric function is integrated on a toric section.

In addition, two experiments that target the internal subdivision scheme of
the method are conducted. In the first experiment, a level set is inserted
between the surfaces of a thin cylindrical sheet to construct a high-order
surface quadrature rule without requiring internal subdivisions. The thin
cylindrical sheet is defined by a rational level set and is immersed in a
coarse grid which otherwise would require the method to subdivide excessively.
In practical applications, internal subdivisions should be avoided because
they cause a high number of quadrature nodes which increases computational
cost. Moreover, they reduce accuracy when a low-order fallback method is
applied to limit the subdivision depth. In the second experiment, a level set
is added to reduce the internal subdivisions required for integration on a
wavy cylinder. This shows that integrating on domains defined by two level
sets can be advantageous, even for simple geometries.

Beyond the close investigation of the method’s general computational cost,
there are two future directions important for practical applications. First,
since the method draws from quadrature methods defined on hypercubes, a wide
range of methods for adaptive subdivision could be applied to create
quadrature rules taylored to specific integrands. Adaptive quadrature enables
quadrature rules with a preset accuracy which often is more pratical than
providing a certain convergence order. Second, the dependency on low-order
fallback methods could be removed by inserting level sets algorithmically.
Especially when integration is performed many times on a static domain, the
reduced number of quadrature nodes could lead to a reduction of computation
time.

Another promising goal is to construct mappings which add minimal difficulty
to the transformed integral. So far, each mapping is created by linearly
interpolating between the boundaries of the integral’s domain. When the
quadrature rule is subsequently mapped to the domain of the integral, the
additional term appearing in the integral can be unnecessarily complex. This
could be mitigated by choosing a different interpolation approach to increase
accuracy while keeping the same number of quadrature nodes.

In an upcoming publication, the method of this paper will be combined with an
extended discontinuous Galerkin method to investigate multi-phase problems
with contact lines by the authors. Accurately modeling the complex behaviour
of dynamic contact lines is a problem of ongoing research which could benefit
from precisely resolved geometries.

###  6.1 Acknowledgements

This research was funded by Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) - 422800359, and supported by the Graduate School CE
within the Centre for Computational Engineering at TU Darmstadt.

## References

  * [1] Hong, Q., Wang, F., Wu, S., Xu, J., 2019. A unified study of continuous and discontinuous Galerkin methods. Sci. China Math. 62, 1-32. https://doi.org/10.1007/s11425-017-9341-1
  * [2] Nguyen, V.P., Anitescu, C., Bordas, S.P.A., Rabczuk, T., 2015. Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation 117, 89-116. https://doi.org/10.1016/j.matcom.2015.05.008
  * [3] Gottlieb, D., Turkel, E., 1983. Spectral methods for time dependent partial differential equations. Numerical Methods in Fluid Dynamics 1127.
  * [4] Olsson, E., Kreiss, G., 2005. A conservative level set method for two phase flow. Journal of Computational Physics 210, 225-246. https://doi.org/10.1016/j.jcp.2005.04.007
  * [5] Osher, S., Fedkiw, R.P., 2001. Level Set Methods: An Overview and Some Recent Results. Journal of Computational Physics 169, 463-502. https://doi.org/10.1006/jcph.2000.6636
  * [6] Annavarapu, C., Hautefeuille, M., Dolbow, J.E., 2012. A robust Nitsche’s formulation for interface problems. Computer Methods in Applied Mechanics and Engineering 225-228, 44-54. https://doi.org/10.1016/j.cma.2012.03.008
  * [7] Kummer F. Extended discontinuous Galerkin methods for two-phase flows: the spatial discretization. International Journal for Numerical Methods in Engineering. 2017;109(2):259-289. doi:https://doi.org/10.1002/nme.5288
  * [8] Groß, S., Reusken, A., 2007. An extended pressure finite element space for two-phase incompressible flows with surface tension. Journal of Computational Physics, Special Issue Dedicated to Professor Piet Wesseling on the occasion of his retirement from Delft University of Technology 224, 40-58. https://doi.org/10.1016/j.jcp.2006.12.021
  * [9] Wang, Y., Oberlack, M., 2011. A thermodynamic model of multiphase flows with moving interfaces and contact line. Continuum Mech. Thermodyn. 23, 409. https://doi.org/10.1007/s00161-011-0186-9
  * [10] Wen, X., 2010. High Order Numerical Methods to Three Dimensional Delta Function Integrals in Level Set Methods. SIAM J. Sci. Comput. 32, 1288-1309. https://doi.org/10.1137/090758295
  * [11] Wen, X., 2009. High order numerical methods to two dimensional delta function integrals in level set methods. Journal of Computational Physics 228, 4273-4290. https://doi.org/10.1016/j.jcp.2009.03.004
  * [12] Wen, X., 2008. High Order Numerical Quadratures to One Dimensional Delta Function Integrals. SIAM J. Sci. Comput. 30, 1825-1846. https://doi.org/10.1137/070682939
  * [13] Drescher, L., Heumann, H., Schmidt, K., 2017. A High Order Method for the Approximation of Integrals Over Implicitly Defined Hypersurfaces. SIAM J. Numer. Anal. 55, 2592-2615. https://doi.org/10.1137/16M1102227
  * [14] Lehrenfeld, C., 2016. High order unfitted finite element methods on level set domains using isoparametric mappings. Computer Methods in Applied Mechanics and Engineering 300, 716-733. https://doi.org/10.1016/j.cma.2015.12.005
  * [15] Bochkov, D., Gibou, F., 2019. Solving Poisson-type equations with Robin boundary conditions on piecewise smooth interfaces. Journal of Computational Physics 376, 1156-1198. https://doi.org/10.1016/j.jcp.2018.10.020
  * [16] Cheng, K.W., Fries, T.-P., 2010. Higher-order XFEM for curved strong and weak discontinuities. International Journal for Numerical Methods in Engineering 82, 564-590. https://doi.org/10.1002/nme.2768
  * [17] Fries, T.-P., Omerović, S., 2016. Higher-order accurate integration of implicit geometries. International Journal for Numerical Methods in Engineering 106, 323-371. https://doi.org/10.1002/nme.5121
  * [18] Fries, T.P., Omerović, S., Schöllhammer, D., Steidl, J., 2017. Higher-order meshing of implicit geometries-Part I: Integration and interpolation in cut elements. Computer Methods in Applied Mechanics and Engineering 313, 759-784. https://doi.org/10.1016/j.cma.2016.10.019
  * [19] Pan, S., Hu, X., Adams, N.A., 2018. A Consistent Analytical Formulation for Volume Estimation of Geometries Enclosed by Implicitly Defined Surfaces. SIAM J. Sci. Comput. 40, A1523-A1543. https://doi.org/10.1137/17M1126370
  * [20] Engvall, L., Evans, J.A., 2017. Isogeometric unstructured tetrahedral and mixed-element Bernstein-Bézier discretizations. Computer Methods in Applied Mechanics and Engineering 319, 83-123. https://doi.org/10.1016/j.cma.2017.02.017
  * [21] Antolin, P., Wei, X., Buffa, A., 2022. Robust numerical integration on curved polyhedra based on folded decompositions. Computer Methods in Applied Mechanics and Engineering 395, 114948. https://doi.org/10.1016/j.cma.2022.114948
  * [22] Gunderman, D., Weiss, K., Evans, J.A., 2021b. Spectral Mesh-Free Quadrature for Planar Regions Bounded by Rational Parametric Curves. Computer-Aided Design 130, 102944. https://doi.org/10.1016/j.cad.2020.102944
  * [23] Gunderman, D., Weiss, K., Evans, J.A., 2021a. High-accuracy mesh-free quadrature for trimmed parametric surfaces and volumes. Computer-Aided Design 141, 103093. https://doi.org/10.1016/j.cad.2021.103093
  * [24] Müller, B., Kummer, F., Oberlack, M., 2013. Highly accurate surface and volume integration on implicit domains by means of moment-fitting. International Journal for Numerical Methods in Engineering 96, 512-528. https://doi.org/10.1002/nme.4569
  * [25] Schwartz, P., Percelay, J., Ligocki, T., Johansen, H., Graves, D., Devendran, D., Colella, P., Ateljevich, E., 2015. High-accuracy embedded boundary grid generation using the divergence theorem. Communications in Applied Mathematics and Computational Science 10, 83-96. https://doi.org/10.2140/camcos.2015.10.83
  * [26] Sommariva, A., Vianello, M., 2009. Gauss-Green cubature and moment computation over arbitrary geometries. Journal of Computational and Applied Mathematics 231, 886-896. https://doi.org/10.1016/j.cam.2009.05.014
  * [27] Bui, H.-G., Schillinger, D., Meschke, G., 2020. Efficient cut-cell quadrature based on moment fitting for materially nonlinear analysis. Computer Methods in Applied Mechanics and Engineering 366, 113050. https://doi.org/10.1016/j.cma.2020.113050
  * [28] Hubrich, S., Düster, A., 2019. Numerical integration for nonlinear problems of the finite cell method using an adaptive scheme based on moment fitting. Computers & Mathematics with Applications 77, 1983-1997. https://doi.org/10.1016/j.camwa.2018.11.030
  * [29] Saye RI. High-Order Quadrature Methods for Implicitly Defined Surfaces and Volumes in Hyperrectangles. SIAM J Sci Comput. 2015;37(2):A993-A1019. doi:10.1137/140966290
  * [30] Cui, T., Leng, W., Liu, H., Zhang, L., Zheng, W., 2020. High-order Numerical Quadratures in a Tetrahedron with an Implicitly Defined Curved Interface. ACM Trans. Math. Softw. 46, 3:1-3:18. https://doi.org/10.1145/3372144
  * [31] Saye, R.I., 2022. High-order quadrature on multi-component domains implicitly defined by multivariate polynomials. Journal of Computational Physics 448, 110720. https://doi.org/10.1016/j.jcp.2021.110720
  * [32] Beck L., F. Kummer, 2023. IntersectingQuadrature, GitHub repository. https://github.com/RumBecken/IntersectingQuadrature
  * [33] Shampine LF. Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics. 2008;211(2):131-140. doi:10.1016/j.cam.2006.11.021
  * [34] Gonnet P. A review of error estimation in adaptive quadrature. ACM Comput Surv. 2012;44(4):22:1-22:36. doi:10.1145/2333112.2333117
  * [35] Beck L., F. Kummer, 2023. Data of Numerical Experiments, TU datalib item. https://doi.org/10.48328/tudatalib-1175

##  Appendix A Hessian

The hessian $H_{\tilde{\beta}}(\tilde{x},\tilde{y},\tilde{z})$ of the
composition
$\tilde{\beta}(\tilde{\boldsymbol{x}})=\beta(A(\tilde{\boldsymbol{x}}))$,

$H_{\tilde{\beta}}(\tilde{\boldsymbol{x}})=H_{\beta}(A(\tilde{\boldsymbol{x}}))=DA^{T}\,H_{\beta}\,DA+\sum_{i\in\\{x,y,z\\}}\partial_{i}\beta
H_{i},$ (115)

requires the hessians of $H_{i}$ of $A$.

Exhibiting increasing complexity, the Hessians $H_{x},H_{y}$ and $H_{z}$
reflect the recursive structure of the nested mapping $M$. While the Hessian
$H_{x}$ is composed of zeros,

$H_{x}=\begin{pmatrix}0&0&0\\\ 0&0&0\\\ 0&0&0\end{pmatrix},$ (116)

the hessian $H_{y}$,

$H_{y}=\frac{1}{2}\begin{pmatrix}M_{\tilde{x}xy}\,\partial_{\tilde{x}}x&(\partial_{x}m_{y}^{\uparrow}-\partial_{x}m_{y}^{\downarrow})\partial_{\tilde{x}}x&0\\\
(\partial_{x}m_{y}^{\uparrow}-\partial_{x}m_{y}^{\downarrow})\partial_{\tilde{x}}x&0&0\\\
0&0&0\\\
\end{pmatrix}=\begin{pmatrix}\partial_{\tilde{x}\tilde{x}}y&\partial_{\tilde{x}\tilde{y}}y&0\\\
\partial_{\tilde{y}\tilde{x}}y&0&0\\\ 0&0&0\\\ \end{pmatrix}.$ (117)

depends on the entry $\partial_{\tilde{x}}x$ of the nested mapping’s Jacobian
matrix and on first partial derivatives of its height functions. In addition,
its entry $M_{\tilde{x}xy}$ requires the second partial derivatives
$\partial_{xx}m_{y}^{\uparrow}$ and $\partial_{xx}m_{y}^{\downarrow}$,

$M_{\tilde{x}xy}=\partial_{\tilde{x}}M_{xy}=\left((\partial_{xx}m_{y}^{\uparrow}-\partial_{xx}m_{y}^{\downarrow})\tilde{y}+\partial_{xx}m_{y}^{\uparrow}+\partial_{xx}m_{y}^{\downarrow}\right)\partial_{\tilde{x}}x.$
(118)

Finally, the hessian $H_{z}$ consists of entries of $H_{y}$ and $DM$,

$H_{z}=\frac{1}{2}\begin{pmatrix}M_{\tilde{x}xz}\partial_{\tilde{x}}x+M_{\tilde{x}yz}\partial_{\tilde{x}}y+M_{yz}\partial_{\tilde{x}\tilde{x}}y&M_{\tilde{x}yz}\partial_{\tilde{y}}y+M_{yz}\partial_{\tilde{x}\tilde{y}}y&M_{\tilde{x}\tilde{z}z}\\\
M_{\tilde{x}yz}\partial_{\tilde{y}}y+M_{yz}\partial_{\tilde{x}\tilde{y}}y&M_{\tilde{y}yz}\partial_{\tilde{y}}y&(\partial_{y}m_{z}^{\uparrow}-\partial_{y}m_{z}^{\downarrow})\partial_{\tilde{y}}y\\\
M_{\tilde{x}\tilde{z}z}&(\partial_{y}m_{z}^{\uparrow}-\partial_{y}m_{z}^{\downarrow})\partial_{\tilde{y}}y&0\\\
\end{pmatrix}.$ (119)

Its entries $M_{\tilde{x}xz},M_{\tilde{x}yz},M_{\tilde{y}yz}$ and
$M_{\tilde{x}\tilde{z}z}$

$\displaystyle M_{\tilde{x}xz}$
$\displaystyle=\partial_{\tilde{x}}M_{xz}(x,y,\tilde{z})$
$\displaystyle=\left((\partial_{xx}m_{z}^{\uparrow}-\partial_{xx}m_{z}^{\downarrow})\tilde{z}+\partial_{xx}m_{z}^{\uparrow}+\partial_{xx}m_{z}^{\downarrow}\right)\partial_{\tilde{x}}x+\left((\partial_{xy}m_{z}^{\uparrow}-\partial_{xy}m_{z}^{\downarrow})\tilde{z}+\partial_{xy}m_{z}^{\uparrow}+\partial_{xy}m_{z}^{\downarrow}\right)\partial_{\tilde{x}}y,$
(120) $\displaystyle M_{\tilde{x}yz}$
$\displaystyle=\partial_{\tilde{x}}M_{yz}(x,y,\tilde{z})$
$\displaystyle=\left((\partial_{xy}m_{z}^{\uparrow}-\partial_{xy}m_{z}^{\downarrow})\tilde{z}+\partial_{xy}m_{z}^{\uparrow}+\partial_{xy}m_{z}^{\downarrow}\right)\partial_{\tilde{x}}x+\left((\partial_{yy}m_{z}^{\uparrow}-\partial_{yy}m_{z}^{\downarrow})\tilde{z}+\partial_{yy}m_{z}^{\uparrow}+\partial_{yy}m_{z}^{\downarrow}\right)\partial_{\tilde{x}}y,$
(121) $\displaystyle M_{\tilde{y}yz}$
$\displaystyle=\partial_{\tilde{y}}M_{yz}(x,y,\tilde{z})$
$\displaystyle=\left((\partial_{xy}m_{z}^{\uparrow}-\partial_{xy}m_{z}^{\downarrow})\tilde{z}+\partial_{xy}m_{z}^{\uparrow}+\partial_{xy}m_{z}^{\downarrow}\right)\partial_{\tilde{y}}x+\left((\partial_{yy}m_{z}^{\uparrow}-\partial_{yy}m_{z}^{\downarrow})\tilde{z}+\partial_{yy}m_{z}^{\uparrow}+\partial_{yy}m_{z}^{\downarrow}\right)\partial_{\tilde{y}}y,$
(122) $\displaystyle M_{\tilde{x}\tilde{z}z}$
$\displaystyle=\partial_{\tilde{x}\tilde{z}}z(x,y,\tilde{z})$
$\displaystyle=(\partial_{x}m_{z}^{\uparrow}-\partial_{x}m_{z}^{\downarrow})\partial_{\tilde{x}}x+(\partial_{y}m_{z}^{\uparrow}-\partial_{y}m_{z}^{\downarrow})\partial_{\tilde{x}}y,$
(123)

contain the second partial derivatives, $\partial_{xx}m_{z}^{\uparrow},$
$\partial_{xx}m_{z}^{\downarrow}$, $\partial_{xy}m_{z}^{\uparrow}$ and
$\partial_{xy}m_{z}^{\downarrow}$, of the height functions.

The second derivatives needed in the hessians of $A$ are obtained after once
more differentiating and rearranging Equation 56 and Equation 58:

$\displaystyle\partial_{xx}a_{y}$
$\displaystyle=-\frac{\partial_{xx}\alpha+2\partial_{xy}\alpha\,\partial_{x}a_{y}+\partial_{yy}\alpha\,(\partial_{x}a_{y})^{2}}{\partial_{y}\alpha},\quad\partial_{xx}a_{z}=-\frac{\partial_{xx}\alpha+2\partial_{xz}\alpha\,\partial_{x}a_{z}+\partial_{zz}\alpha\,(\partial_{x}a_{z})^{2}}{\partial_{z}\alpha},$
(124) $\displaystyle\partial_{yy}a_{z}$
$\displaystyle=-\frac{\partial_{yy}\alpha+2\partial_{yz}\alpha\,\partial_{y}a_{z}+\partial_{zz}\alpha\,(\partial_{y}a_{z})^{2}}{\partial_{z}\alpha},\quad\partial_{xy}a_{z}=-\frac{\partial_{xy}\alpha+\partial_{xz}\alpha\,\partial_{y}a_{z}+\partial_{yz}\alpha\,\partial_{x}a_{z}+\partial_{zz}\alpha\,\partial_{x}a_{z}\partial_{y}a_{z}}{\partial_{z}\alpha}.$
(125)


* Yuan et al. (2023a) Yuan, S., Zamora, A., & Abel, T. 2023a, MNRAS, 522, 3935, doi: 10.1093/mnras/stad1275
  * Yuan et al. (2023b) Yuan, S., Zhang, H., Ross, A. J., et al. 2023b, arXiv e-prints, arXiv:2306.06314, doi: 10.48550/arXiv.2306.06314
  * Zennaro et al. (2022) Zennaro, M., Angulo, R. E., Contreras, S., Pellejero-Ibáñez, M., & Maion, F. 2022, MNRAS, 514, 5443, doi: 10.1093/mnras/stac1673
  * Zennaro et al. (2021) Zennaro, M., Angulo, R. E., Pellejero-Ibáñez, M., et al. 2021, arXiv e-prints, arXiv:2101.12187.  https://arxiv.org/abs/2101.12187
  * Zhai et al. (2023a) Zhai, Z., Percival, W. J., & Guo, H. 2023a, MNRAS, 523, 5538, doi: 10.1093/mnras/stad1793
  * Zhai et al. (2023b) Zhai, Z., Tinker, J. L., Banerjee, A., et al. 2023b, ApJ, 948, 99, doi: 10.3847/1538-4357/acc65b
  * Zheng et al. (2007) Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760, doi: 10.1086/521074
  * Zheng et al. (2005) Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791, doi: 10.1086/466510