## 2 Regularization in the Level Set Context  
Let the volume \(V\) be divided into \(V_{1}\) and \(V_{2}\) by the interface \(\Gamma\). At any point \(\mathbf{x}\in\Gamma\), the interface is characterized by the normal vector \(\mathbf{n}(\mathbf{x})\). The position of \(\Gamma\) is implicitly defined by means of the signed distance function with respect to \(\Gamma\)[20]  
\[d(\mathbf{x})=\begin{cases}\text{dist}(\mathbf{x},\Gamma)&\text{if}\qquad \mathbf{x}\in V_{1}\\ -\text{dist}(\mathbf{x},\Gamma)&\text{if}\qquad\mathbf{x}\in V_{2}\end{cases} \tag{1}\]  
We introduce the Heaviside function \(\mathcal{H}(d(\mathbf{x}))\) such that \(\mathcal{H}(d(\mathbf{x}))=1\) for \(d(\mathbf{x})>0\) and \(\mathcal{H}(d(\mathbf{x}))=-1\) for \(d(\mathbf{x})\leq 0\), and focus on discontinuous functions of the type  
\[\mathbf{u}(\mathbf{x})=\mathbf{v}(\mathbf{x})+\mathcal{H}(d(\mathbf{x}))\mathbf{ j}(\mathbf{x}) \tag{2}\]  
where \(\mathbf{v}(\mathbf{x})\) is the regular part of the field \(\mathbf{u}\), while \(\mathbf{j}(\mathbf{x})\) denotes the amplitude of the jump of \(\mathbf{u}(\mathbf{x})\) across \(\Gamma\). Differentiation of the discontinuous vector function (2) leads to a singular function of the type  
\[\boldsymbol{\varepsilon}(\mathbf{x})=\
abla\mathbf{u}(\mathbf{x})=\
abla \mathbf{v}(\mathbf{x})+\mathcal{H}(\mathbf{x})\
abla\mathbf{j}(\mathbf{x})+ \delta(d(\mathbf{x}))(\mathbf{j}(\mathbf{x})\otimes\mathbf{n}(\mathbf{x})) \tag{3}\]  
where \(\
abla\) denotes differentiation. The delta function \(\delta(d(\mathbf{x}))\) is introduced through the generalized gradient \(\
abla\mathcal{H}(d(\mathbf{x}))=\delta(d(\mathbf{x}))\mathbf{n}(\mathbf{x})\).  
Our purpose is to numerically integrate terms containing either discontinuous functions such as Equation (2), or singular functions such as Equation (3). In order to avoid quadrature subdomains, we want to replace the Heaviside function \(\mathcal{H}\) by a regularized Heaviside function \(\mathcal{H}_{\rho}\) and the delta function \(\delta\) by a regularized version \(\delta_{\rho}\). In the one-dimensional case, denoting with \(t\) the scalar distance from the discontinuity, regularization can be carried out by replacing the delta function \(\boldsymbol{\delta}(d(\mathbf{x}))\) with [18, 23]  
\[\delta_{\rho}(t)=\frac{1}{V_{\rho}}\varphi(\frac{t}{\rho})\,,\qquad\lim_{\rho \to 0}\delta_{\rho}(t)=\delta(t) \tag{4}\]  
where the weight function \(\varphi\) is a smooth function depending on a length \(\rho\) such that  
\[\int_{-\infty}^{\infty}\delta_{\rho}(\xi)\,d\xi=\frac{1}{V_{\rho}}\int_{- \infty}^{\infty}\varphi(\xi)\,d\xi=1 \tag{5}\]The normalization parameter \(V_{\rho}\) depends on the type of the weight function and the dimension of the space [17]. Examples of weight functions are shown in Table 1; other examples can be found in [24]. The one-dimensional regularized Heaviside function is then obtained as  
\[\mathcal{H}_{\rho}(t)=\int_{-\infty}^{t}\delta_{\rho}(t_{1})\,dt_{1} \tag{6}\]  
A straightforward extension to the multi-dimensional case can be carried out by assuming the regularized delta function to be a function of the signed distance function \(d\) (1)  
\[\delta_{\rho}(d)=\frac{1}{V_{\rho}}\varphi(\frac{d}{\rho}) \tag{7}\]  
In this paper, the level set definition (7) is adopted.