# A penalized $\phi$-FEM scheme for the Poisson Dirichlet problem

Michel Duprez, Killian Vuillemot, Raphaël Bulle, Vanessa Lleras

<figure><img src="image_1.png" /></figure>

## To cite this version:

Michel Duprez, Killian Vuillemot, Raphaël Bulle, Vanessa Lleras. A penalized $\phi$-FEM scheme for the Poisson Dirichlet problem. 2026. hal-05495474

HAL Id: hal-05495474

https://hal.science/hal-05495474v1

Preprint submitted on 5 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A penalized $\varphi$-FEM scheme for the Poisson Dirichlet problem

Raphaël Bulle*

Michel Duprez†

Vanessa Lleras‡

Killian Vuillemot§

## Abstract

In this work, we analyze a penalized variant of the $\varphi$-FEM scheme for the Poisson equation with Dirichlet boundary conditions. The $\varphi$-FEM is a recently introduced unfitted finite element method based on a level-set description of the geometry, which avoids the need for boundary-fitted meshes. Unlike the original $\varphi$-FEM formulation, the method proposed here enforces boundary conditions through a penalization term. This approach has the advantage that the level-set function is required only on the cells adjacent to the boundary in the variational formulation. The scheme is stabilized using a ghost penalty technique. We derive a priori error estimates, showing optimal convergence in the $H^1$ semi-norm and quasi-optimal convergence in the $L^2$ norm under suitable regularity assumptions. Numerical experiments are presented to validate the theoretical results and to compare the proposed method with both the original $\varphi$-FEM and the standard fitted finite element method.

**Keywords.** Immersed boundary method; complex geometries; partial differential equations; level-set methods.

## 1 Introduction

The finite element method (FEM) usually requires meshes that conform to the domain boundary, which can be restrictive for complex or evolving geometries. To address this limitation, Immersed Boundary Methods (IBM) [14] allow the use of meshes that do not fit the physical boundary by extending the problem to a larger domain. More recently, CutFEM [3, 4] has been proposed as a robust alternative, based on integrating the weak formulation only over the parts of elements cut by the boundary and adding stabilization terms to ensure numerical stability and coercivity. Another alternative, called the Shifted Boundary Method [13], uses a Taylor expansion of the boundary condition. The $\varphi$-FEM paradigm introduced in [12] is based on a level-set description of the domain, where the solution is written as $u = \varphi w$, automatically enforcing boundary conditions and achieving optimal accuracy and good conditioning without requiring boundary or cut-cell integration. This method has been extended to various PDEs and interface problems [5, 6, 7, 11]. More recently, this approach has been adapted to a finite difference scheme [9], combined with a neural network [10] and an a posteriori estimator, and has been analysed to estimate the error [2].

<footnote>* Université de Strasbourg, CNRS, Inria, ICube, F-67000 Strasbourg, France. raphael.bulle@inria.fr</footnote>

<footnote>† Université de Strasbourg, CNRS, Inria, ICube, F-67000 Strasbourg, France. michel.duprez@inria.fr</footnote>

<footnote>‡IMAG, Univ Montpellier, CNRS UMR 5149, 34090 Montpellier, France. vanessa.lleras@umontpellier.fr</footnote>

<footnote>§IMAG, Univ Montpellier, CNRS UMR 5149, 34090 Montpellier, France. killian.vuillemot@umontpellier.fr</footnote>
In this paper, we consider the Poisson equation with Dirichlet boundary conditions,

$$ \begin{cases} -\Delta u = f, & \text{in } \Omega, \\ u = 0, & \text{on } \Gamma, \end{cases} \qquad (1) $$

The present contribution is devoted to the numerical analysis of a $\varphi$-FEM scheme introduced in [5] for this problem (see Section 2). This $\varphi$-FEM scheme enforces boundary conditions weakly by introducing a penalization on the unfitted boundary while still using a level-set description of the domain. So this method combines the robustness of mixed formulations with the flexibility of immersed finite element methods. We define $\Omega \subset \mathbb{R}^d$ (here $d = 2, 3$) a domain of boundary $\Gamma$, $f \in L^2(\Omega)$ and $n$ the outward unit normal to $\Omega$. We assume that the domain and its boundary are given by a level-set function $\varphi$ such that

$$ \Omega := \{\varphi < 0\} \quad \text{and} \quad \Gamma := \{\varphi = 0\}. \qquad (2) $$

The present paper is organized as follows: in Section 2, we present the scheme and the main theoretical results. In the second and third Section, we demonstrate these results. Section 4 analyzes the conditioning number of associated matrix. Finally, Section 5 is devoted to numerical illustrations.

## 2 Proposed scheme and theoretical results

Let $\mathcal{T}_h^O$ be a Cartesian grid covering the domain $\mathcal{O}$, with triangular cells of size $h$. Let $\varphi_h = I_{h,\mathcal{O}}^{(l)}\varphi$ the standard Lagrange continuous interpolation of $\varphi$ (of degree $l > 0$) on $\mathcal{T}_h^O$. This interpolation is then used to construct a submesh $\mathcal{T}_h$ of $\mathcal{T}_h^O$ collecting all the cells of the grid intersecting the domain $\{\varphi_h < 0\}$, i.e.

$$ \mathcal{T}_h := \{T \in \mathcal{T}_h^O : T \cap \{\varphi_h < 0\} \neq \emptyset\}. $$

Let us finally introduce another submesh of $\mathcal{T}_h$, collecting the cells intersecting the boundary, i.e. $\{\varphi_h = 0\}$ and a set containing the facets of this submesh, given by

$$ \mathcal{T}_h^\Gamma := \{T \in \mathcal{T}_h \ T \cap \{\varphi_h = 0\} \neq \emptyset\}, \quad \text{and} \quad \mathcal{F}_h^\Gamma := \{F \in \mathcal{T}_h^\Gamma \setminus \partial\Omega_h\}. \qquad (3) $$

Let us now consider the finite element spaces: for $k \ge 1$,

$$ V_h^{(k)} := \{v_h \in H^1(\Omega_h) \ v_h|_T \in \mathbb{P}_k(T) \ \forall T \in \mathcal{T}_h\}, \qquad (4) $$

and a second one, the local version of the previous space, on $\mathcal{T}_h^\Gamma$,

$$ Q_h^{(k)} := \{q_h \in L^2(\Omega_h^\Gamma) \ v_h|_T \in \mathbb{P}_k(T) \ \forall T \in \mathcal{T}_h^\Gamma\}. $$

We will now present the $\varphi$-FEM penalized scheme analyzed here, that relies on the $\varphi$-FEM paradigm used only in the cells of $\Omega_h^\Gamma$. For this scheme, we thus introduce an auxiliary variable $p$ on $\Omega_h^\Gamma$ such that $u = \varphi p$; equation that will be enforced by penalization in the following scheme.

The penalized $\varphi$-FEM scheme to solve (1) is given by: find $u_h \in V_h^{(k)}$ and $p_h \in Q_h^{(k)}$ such that

$$ \int_{\Omega_h} \nabla u_h \cdot \nabla v_h - \int_{\partial\Omega_h} \frac{\partial u_h}{\partial n} v_h + \frac{\gamma}{h^2} \int_{\Omega_h^\Gamma} (u_h - \frac{1}{h}\varphi_h p_h)(v_h - \frac{1}{h}\varphi_h q_h) \\ + G_h^{lhs}(u_h, v_h) = \int_{\Omega_h} f v_h + G_h^{rhs}(v_h), \qquad (5) $$
for all $v_h \in V_h^{(k)}$, $q_h \in Q_h^{(k)}$ with

$$G_h^{lhs}(u, v) = \sigma_D h \sum_{E \in \mathcal{F}_h^\Gamma} \int_E \left[ \frac{\partial u}{\partial n} \right] \left[ \frac{\partial v}{\partial n} \right] + \sigma_D h^2 \sum_{T \in \mathcal{T}_h^\Gamma} \int_T \Delta u \Delta v, \quad (6)$$

and

$$G_h^{rhs}(v) = -\sigma_D h^2 \sum_{T \in \mathcal{T}_h^\Gamma} \int_T f \Delta v \quad (7)$$

for some $\gamma$, $\sigma_D > 0$. As said in the introduction, this scheme has been introduced in [5], but not analyzed.

Remark 1. We can remark three important points:

* In the case of $\mathbb{P}^1$ finite elements, the second-order terms (the ones multiplied by $\sigma_D h^2$) are not necessary.

* In the case of non-homogeneous Dirichlet boundary conditions $u = u_D$, all you need to do is the modification of the penalization form, given by $u = \varphi p + u_D$. It then remains to adapt the scheme to this equation.

* A final observation concerns the comparison with the direct $\varphi$-FEM scheme [12], relying on the assumption $u = \varphi w$ over the entire mesh, thereby redefining the problem's unknown throughout the domain. This has two major consequences:

- the computation relies on the level-set only near the boundary, which makes the scheme more robust to the level-set and, in particular, does not require the gradient or the second derivatives of the level-set function, which could present a singularity (for example, if we consider the signed distance function).

- the scheme is compatible with the Neumann approach presented in [6].

Let us now recall some important hypothesis on the mesh and lemmas from [12] and [8] and introduce our main theoretical results.

**Assumption 1.** The boundary $\Gamma$ can be covered by open sets $O_i$, $i = 1, \dots, I$ on which we can introduce local coordinates $\xi_1, \dots, \xi_d$ with $\xi_d = \varphi$ such that, up to the order $k+1$, all the partial derivatives $\partial^\alpha \xi_i / \partial x^\alpha$ and $\partial x^\alpha / \partial \xi_i$ are bound by a constant $C_0 > 0$. Moreover, on $O$, $\varphi$ is of class $C^{k+1}$ and there exists $m > 0$ such that $|\varphi| \ge m$ on $O \setminus \cup_{i=1, \dots, I} O_i$.

**Assumption 2.** The approximated boundary, defined by $\Gamma_h = \{\varphi_h = 0\}$, can be covered by patches of elements $\{\Pi_r\}_{r=1, \dots, N_\Pi}$ such that:

* Each patch $\Pi_r$ can be written as $\Pi_r = \Pi_r^\Gamma \cup T_r$ where $\Pi_r^\Gamma \subset \mathcal{T}_h^\Gamma$ and $T_r \in \mathcal{T}_h \setminus \mathcal{T}_h^\Gamma$. Moreover $\Pi_r$, contains at most $M$ connected elements, with $M$ independent of $h$;

* The mesh $\mathcal{T}_h^\Gamma$ satisfies $\mathcal{T}_h^\Gamma = \cup_{r=1, \dots, N_\Pi} \Pi_r^\Gamma$;

* Two patches $\Pi_r$ and $\Pi_s$ are disjoints if $r \neq s$.

Let us now present the main theorem of this paper:

**Theorem 1.** We suppose that assumptions 1 and 2 on $\Gamma$ and $\mathcal{T}_h$ are satisfied, $f \in H^{k-1}(\Omega_h)$, $\Omega \subset \Omega_h$ and $\gamma$, $\sigma_D$ large enough. Consider $u \in H^{k+1}(\Omega)$ the solution of (1) and $u_h \in V_h$ the solution of the discrete scheme (5). It holds

$$|u - u_h|_{1,\Omega} \le C h^k \|f\|_{k-1,\Omega_h}$$

with $C > 0$ a constant independent on $h$.
We note that our analysis requires less regularity on $u$ compared to [12]. Moreover, the results of [12] can be adapted to show that the condition number of the associated matrix scales like $O(h^{-2})$, while the $L^2$-error is of order $O(h^{k+1/2})$. In the final section of this paper, optimal convergence of order $O(h^{k+1})$ is observed numerically.

### 3 Coercivity of the bilinear form

This section is devoted to the proof of the coercivity needed for Theorem 1.

**Proposition 1.** Consider the bilinear form given by

$$a_h(u, p; v, q) = \int_{\Omega_h} \nabla u \cdot \nabla v - \int_{\partial \Omega_h} \frac{\partial u}{\partial n} v + \frac{\gamma}{h^2} \int_{\Omega_h^\Gamma} (u - \frac{1}{h} \varphi_h p)(v - \frac{1}{h} \varphi_h q) \\ + \sigma_D h \sum_{F \in \mathcal{F}_h^\Gamma} \int_F \left[ \frac{\partial u}{\partial n} \right] \left[ \frac{\partial v}{\partial n} \right] + \sigma_D h^2 \int_{\Omega_h^\Gamma} \Delta u \Delta v. \quad (8)$$

Under assumptions 1 and 2, for $\gamma$ and $\sigma_D$ sufficiently big, $a_h$ is coercive on $V_h^{(k)} \times Q_h^{(k)}$ for the norm

$$\|(u, p)\|_h^2 = |u|_{1, \Omega_h}^2 + \frac{1}{h^2} \left\| u - \frac{1}{h} \varphi_h p \right\|_{0, \Omega_h^\Gamma}^2 + h \sum_{F \in \mathcal{F}_h^\Gamma} \left\| \left[ \frac{\partial u}{\partial n} \right]_{0, F} \right\|_{0, \Omega_h^\Gamma}^2 + h^2 \| \Delta u \|_{0, \Omega_h^\Gamma}^2, \quad (9)$$

i.e. $\|(u_h, p_h)\|_h^2 \ge C a_h(u_h, p_h; u_h, p_h)$ for each $(u_h, p_h) \in V_h^{(k)} \times Q_h^{(k)}$.

Before proving Proposition 1, we need the following Lemma :

**Lemma 1.** For all $u_h \in V_h^{(k)}$ and $p_h \in Q_h^{(k)}$, there exists $C > 0$ such that

$$\|\varphi_h p_h\|_{0, \Omega_h^\Gamma} \le C(h\|\nabla u_h\|_{0, \Omega_h^\Gamma} + \|u_h - \varphi_h p_h\|_{0, \Omega_h^\Gamma}).$$

*Proof.* Using the Poincaré inequality, combined to the triangular inequality and an inverse inequality, we obtain

$$\begin{aligned} \|\varphi_h p_h\|_{0, \Omega_h^\Gamma} &\le Ch\|\nabla\varphi_h p_h\|_{0, \Omega_h^\Gamma} \le Ch(\|\nabla u_h\|_{0, \Omega_h^\Gamma} + \|\nabla(u_h - \varphi_h p_h)\|_{0, \Omega_h^\Gamma}) \\ &\le C(h\|\nabla u_h\|_{0, \Omega_h^\Gamma} + \|u_h - \varphi_h p_h\|_{0, \Omega_h^\Gamma}). \end{aligned} \quad \square$$

We also need the following Lemma :

**Lemma 2** (cf. [12, Lemma 3.3]). Under the Assumption 2, for all $\beta > 0$, one can choose $\alpha \in ]0, 1[$ depending only on the regularity of the mesh, such that for all $v_h \in V_h^{(k)}$

$$|v_h|_{1, \Omega_h^\Gamma}^2 \le \alpha |v_h|_{1, \Omega_h}^2 + \beta h \sum_{F \in \mathcal{F}_h^\Gamma} \left\| \left[ \frac{\partial v_h}{\partial n} \right]_{0, F} \right\|_{0, \Omega_h^\Gamma}^2 + \beta h^2 \| \Delta v_h \|_{0, \Omega_h^\Gamma}^2.$$

*Proof of Proposition 1.* Let $B_h$ the strip between $\partial \Omega_h$ and $\Gamma_h$, i.e. $B_h = \{\varphi_h > 0\} \cap \Omega_h$. Since $\varphi_h = 0$ on $\Gamma_h$, the boundary term of (8) can be written for all $(u_h, p_h) \in V_h^{(k)} \times Q_h^{(k)}$
as :

$$
\int_{\partial \Omega_h} \frac{\partial u_h}{\partial n} u_h = \overbrace{\int_{B_h} |\nabla u_h|^2}^{I} - \overbrace{\int_{\Gamma_h} \frac{\partial u_h}{\partial n} (u_h - \frac{1}{h} \varphi_h p_h)}^{II} - \underbrace{\sum_{F \in \mathcal{F}_h^\Gamma} \int_{F \cap B_h} \left[ \frac{\partial u_h}{\partial n} \right] u_h}_{III} + \underbrace{\int_{B_h} \Delta(u_h) u_h}_{IV}. \quad (10)
$$

Using Lemma 2, for all $\beta > 0$, there exists $\alpha \in (0, 1)$ such that

$$
I \le |u_h|_{1, \Omega_h^\Gamma}^2 \le \alpha |u_h|_{1, \Omega_h}^2 + \beta h \sum_{F \in \mathcal{F}_h^\Gamma} \left\| \left[ \frac{\partial u_h}{\partial n} \right] \right\|_{0, F}^2 + \beta h^2 \| \Delta u_h \|_{0, \Omega_h^\Gamma}^2.
$$

Moreover, thanks to the trace inequality and an inverse inequality, for all $\varepsilon > 0$,

$$
\begin{align*}
II &\le C \left( \frac{1}{\sqrt{h}} \|\nabla u_h\|_{0,\Omega_h^\Gamma} + \sqrt{h} |\nabla u_h|_{1,\Omega_h^\Gamma} \right) \\
&\quad \times \left( \frac{1}{\sqrt{h}} \left\| u_h - \frac{1}{h} \varphi_h p_h \right\|_{0,\Omega_h^\Gamma} + \sqrt{h} \left| u_h - \frac{1}{h} \varphi_h p_h \right|_{1,\Omega_h^\Gamma} \right) \\
&\le \frac{C}{h} |u_h|_{1,\Omega_h^\Gamma} \left\| u_h - \frac{1}{h} \varphi_h p_h \right\|_{0,\Omega_h^\Gamma} \le C\varepsilon |u_h|_{1,\Omega_h^\Gamma}^2 + \frac{C}{\varepsilon h^2} \left\| u_h - \frac{1}{h} \varphi_h p_h \right\|_{0,\Omega_h^\Gamma}^2.
\end{align*}
$$

Cauchy-Schwarz, Young and trace inequalities as well Lemma 1 lead to

$$
III \le \frac{C}{\varepsilon} h \sum_{F \in \mathcal{F}_h^\Gamma} \left\| \left[ \frac{\partial u_h}{\partial n} \right] \right\|_{0,F}^2 + C\varepsilon |u_h|_{1,\Omega_h^\Gamma}^2 + \frac{C\varepsilon}{h^2} \left\| u_h - \frac{1}{h} \varphi_h p_h \right\|_{0,\Omega_h^\Gamma}^2.
$$

Finally, Lemma 1 gives

$$
IV \le \frac{Ch^2}{\varepsilon} \|\Delta u_h\|_{0,\Omega_h^\Gamma}^2 + C\varepsilon|u_h|_{1,\Omega_h^\Gamma}^2 + \frac{C\varepsilon}{h^2} \left\| u_h - \frac{1}{h}\varphi_h p_h \right\|_{0,\Omega_h^\Gamma}^2.
$$

Hence,

$$
\begin{align*}
a_h(u_h, p_h; u_h, p_h) &\ge (1 - \alpha - C\varepsilon) |u_h|_{1,\Omega_h}^2 \\
&\quad + \left(\sigma_D - \frac{C}{\varepsilon} - \beta\right) \left( h \sum_{F \in \mathcal{F}_h^\Gamma} \left\| \left[ \frac{\partial u_h}{\partial n} \right] \right\|_{0,F}^2 + h^2 \| \Delta u_h \|_{0,\Omega_h^\Gamma}^2 \right) \\
&\quad + \left( \gamma - C\varepsilon - \frac{C}{\varepsilon} \right) \frac{1}{h^2} \left\| u_h - \frac{1}{h} \varphi_h p_h \right\|_{0,\Omega_h^\Gamma}^2.
\end{align*}
$$

Choosing $\varepsilon$ sufficiently small, $\sigma_D$ and $\gamma$ sufficiently big, we obtain the desired result. $\square$

4 Proof of the H$^1$ estimate

To prove Theorem 1, we will need the two following results:
**Lemma 3** (cf. [6, Lemma 3.4]). Under the assumption 1, all function $v \in H^s(\Omega_h)$ vanishing on $\Omega$, with $1 \le s \le k + 1$, satisfies

$$ \|v\|_{0,\Omega_h \setminus \Omega} \le Ch^s \|v\|_{s,\Omega_h \setminus \Omega}. $$

**Lemma 4.** We assume that the domain $\Omega$ is given by the level-set $\phi$, cf. (2), and satisfies Assumption 1. Then, for any $u \in H^{k+1}(\mathcal{O})$ vanishing on $\Gamma$,

$$ \left\| \frac{u}{\phi} \right\|_{k, \mathcal{O}} \le C \|u\|_{k+1, \mathcal{O}} $$

with $C > 0$ depending only on the constants in Assumption 1.

*Proof of Theorem 1, H¹ estimate.* Let $\tilde{u} \in H^{k+1}(\Omega_h)$ an extension of $u$ from $\Omega$ to $\Omega_h$ such that $\tilde{u} = u$ on $\Omega$ and

$$ \|\tilde{u}\|_{k+1, \Omega_h} \le C \|u\|_{k+1, \Omega} \le C \|f\|_{k-1, \Omega}. $$

We consider $\tilde{f} := -\Delta \tilde{u}$ and $p := \frac{h}{\phi} \tilde{u}$ on $\Omega_h$. Then,

$$ a_h(\tilde{u}, p; v_h, q_h) = \int_{\Omega_h} \tilde{f} v_h - \sigma_D h^2 \int_{\Omega_h^\Gamma} \tilde{f} \Delta v_h + \frac{\gamma}{h^2} \int_{\Omega_h^\Gamma} (\tilde{u} - \frac{1}{h} \varphi_h p)(v_h - \frac{1}{h} \varphi_h q_h) $$

for each $(v_h, q_h) \in V_h^{(k)} \times Q_h^{(k)}$. We obtain the following Galerkin orthogonality

$$ a_h(\tilde{u} - u_h, p - p_h; v_h, q_h) = \int_{\Omega_h} (\tilde{f} - f)v_h - \sigma_D h^2 \int_{\Omega_h^\Gamma} (\tilde{f} - f)\Delta v_h + \frac{\gamma}{h^2} \int_{\Omega_h^\Gamma} (\tilde{u} - \frac{1}{h}\varphi_h p)(v_h - \frac{1}{h}\varphi_h q_h). \quad (11) $$

Thus, using the coercivity of the bilinear form (c.f. Proposition 1),

$$ c \|(u_h - I_h \tilde{u}, p_h - I_h p)\|_h \le \sup_{(v_h, q_h)} \frac{a_h(u_h - I_h \tilde{u}, p_h - I_h p; v_h, q_h)}{\|(v_h, q_h)\|_h} \\ \le \sup_{(v_h, q_h)} \frac{I - II - III}{\|(v_h, q_h)\|_h}, $$

with

$$ I = a_h(e_u, e_p; v_h, q_h), \quad II = \int_{\Omega_h} (\tilde{f} - f)v_h - \sigma_D h^2 \int_{\Omega_h^\Gamma} (\tilde{f} - f)\Delta v_h, \\ III = \frac{\gamma}{h^2} \int_{\Omega_h^\Gamma} (\tilde{u} - \frac{1}{h}\varphi_h p)(v_h - \frac{1}{h}\varphi_h q_h), $$

where $e_u = \tilde{u} - I_h \tilde{u}$ and $e_p = p - I_h p$. It now remains to estimate each of the terms. For the term $I$, using the expression (10), the trace inequality [6, Lemma 3.5], interpolation inequalities and the expression (9),

$$ I \le Ch^k \|\tilde{u}\|_{k+1, \Omega_h} \|(v_h, q_h)\|_h \le Ch^k \|f\|_{k-1, \Omega} \|(v_h, q_h)\|_h. $$

Since $\tilde{f} = f$ on $\Omega$ and $\Omega \subset \Omega_h$, using Lemma 3, it holds

$$ \|\tilde{f} - f\|_{0, \Omega_h \setminus \Omega} \le Ch^{k-1} \|\tilde{f} - f\|_{k-1, \Omega_h \setminus \Omega}. $$
Hence, using Lemma 1 and the interpolation inequalities,

$$ II \leq C \|\tilde{f} - f\|_{0,\Omega_h \setminus \Omega} (\|v_h\|_{0,\Omega_h \setminus \Omega} + \sigma h^2 \|\Delta v_h\|_{0,\Omega_h \setminus \Omega}) \\ \leq Ch^k \|f\|_{k-1,\Omega_h} \|(v_h, q_h)\|_h. $$

Let us now estimate term III: Since $\tilde{u} = \frac{1}{h}p\varphi$ on $\Omega_h^\Gamma$,

$$ III \leq \frac{C}{h} \|\tilde{u} - \frac{1}{h}\varphi_h p\|_{0,\Omega_h^\Gamma} \|(v_h, q_h)\|_h \leq \frac{C}{h} \|\varphi - \varphi_h\|_{\infty,\Omega_h^\Gamma} \frac{p}{h} \|v_h, q_h\|_h. $$

Using Lemma 4 and interpolation inequalities,

$$ III \leq Ch^k \|\varphi\|_{W_{k+1}^\infty} \|\tilde{u}\|_{1,\Omega_h^\Gamma} \|(v_h, q_h)\|_h \\ \leq Ch^k \|f\|_{k-1,\Omega} \|(v_h, q_h)\|_h. $$

Using the previous estimates and by definition of $\|\cdot\|_h$,

$$ |u_h - I_h \tilde{u}|_{1,\Omega_h} \leq \|(u_h - I_h \tilde{u}, p_h - I_h p)\|_h \leq Ch^k \|f\|_{k-1,\Omega}. $$

Finally, using the triangle inequality and some interpolation estimates,

$$ |u_h - u|_{1,\Omega} \leq |u_h - I_h \tilde{u}|_{1,\Omega_h} + |I_h \tilde{u} - \tilde{u}|_{1,\Omega_h} \leq Ch^k \|f\|_{k-1,\Omega}. $$

<figure><img src="image_3.png" /></figure>

# 5 Numerical illustration

We conclude this study by numerical test cases implemented in Python, using the Dolfinx [1] library. For these two test cases, we compare our scheme to a standard finite element method and to the first $\varphi$-FEM scheme presented in [12], respectively referred as "Dual $\varphi$-FEM", "Standard-FEM" and "Direct $\varphi$-FEM" in the figures.

## 5.1 First test case: a 2D complex geometry

We first consider a complex geometry (see Fig. 1) defined by

$$ \varphi(x, y) = (-1)^5 \prod_{j=1}^{5} \left( -1 + \exp \left( -\frac{x_j^2}{2l_{x,j}^2} - \frac{y_j^2}{2l_{y,j}^2} \right) \right) - 0.5, \quad (12) $$

with

$$ x_j = \cos(\theta_j)(x - x_{0,j}) - \sin(\theta_j)(y - y_{0,j}) \text{ and } y_j = \sin(\theta_j)(x - x_{0,j}) + \cos(\theta_j)(y - y_{0,j}), $$

for some given parameters $(x_{0,j}, y_{0,j}, l_{x,j}, l_{y,j}, \theta_j)$ given in Table 1. This expression was chosen to fit a liver geometry in an other study.

The source term is given by $f(x, y) = \cos(x) \exp(y)$, and homogeneous Dirichlet boundary conditions are imposed on $\Gamma$. The parameters of the penalized $\varphi$-FEM scheme are set to $\gamma = 100$ and $\sigma_D = 0.1$. The reference solution and the corresponding absolute errors are displayed in Fig. 1.

In Fig. 2, we compare the relative $L^2$ and semi-$H^1$ errors with respect to a fine standard FEM solution computed on a mesh with $h \approx 0.002$. The three methods yield comparable results, each exhibiting optimal convergence rates in both the $L^2$ and $H^1$ norms.
<table><caption>Table 1: Test case 1. Parameters of $\varphi$.</caption><thead><tr><th>j</th><th>x<sub>0,j</sub></th><th>y<sub>0,j</sub></th><th>l<sub>x,j</sub></th><th>l<sub>y,j</sub></th><th>θ<sub>j</sub></th></tr></thead><tbody><tr><td>1</td><td>0.356</td><td>0.507</td><td>0.145</td><td>0.171</td><td>0.000</td></tr><tr><td>2</td><td>0.588</td><td>0.589</td><td>0.153</td><td>0.090</td><td>0.000</td></tr><tr><td>3</td><td>0.569</td><td>0.588</td><td>0.008</td><td>0.008</td><td>0.006</td></tr><tr><td>4</td><td>0.308</td><td>0.443</td><td>0.055</td><td>0.116</td><td>0.622</td></tr><tr><td>5</td><td>0.741</td><td>0.643</td><td>0.058</td><td>0.035</td><td>0.000</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 1: Test case 1. From left to right: reference solution (standard FEM with a fine mesh); difference between the reference solution and the projection of each approximated solution (Standard FEM, Direct $\varphi$-FEM, Dual $\varphi$-FEM).</figcaption></figure>

Furthermore, Fig. 3 (left) compares the total computational time of each method, including mesh generation, cell selection, and assembly and solution of the finite element system. These results show that both $\varphi$-FEM approaches are faster than the standard FEM. However, since the finite element system associated with the penalized scheme is larger due to the introduction of an additional variable, higher computation times are observed compared to the direct $\varphi$-FEM. Finally, Fig. 3 (right) illustrates that optimal conditioning is achieved.

## 5.2 Second test case: a 3D geometry

For this test case, we consider a sphere described by the level-set function

$$ \varphi(x, y, z) = -0.3125^2 + (x - 0.5)^2 + (y - 0.5)^2 + (z - 0.5)^2, $$

and a manufactured solution given by $u_{ex}(x, y, z) = 1 - \exp(\varphi(x, y, z)^2)$ such that $u_{ex}$ vanishes on $\Gamma$. The parameters of the penalized $\varphi$-FEM scheme are set to $\gamma = 100$ and $\sigma_D = 0.01$. The reference solution and the absolute errors are given in Fig. 4. The relative $L^2$ and $H^1$ errors are given in Fig 5, verifying the optimal convergence in both norms.

## Acknowledgment

This work was supported by the Agence Nationale de la Recherche, Project PhiFEM, under grant ANR-22-CE46-0003-01.

## References

[1] I. A. Baratta, J. P. Dean, J. S. Dokken, M. Habera, J. S. Hale, C. N. Richardson, M. E. Rognes, M. W. Scroggs, N. Sime, and G. N. Wells. DOLFINx: the next generation FEniCS problem solving environment. preprint, 2023.
<figure><img src="image_5.png" /><figcaption>Figure 2: **Test case 1.** $L^2$ (left) and semi-$H^1$ (right) relative errors with respect to the mesh size.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3: **Test case 1.** Left: computation time with respect to the mesh size. Right: condition number with respect to the mesh size.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: **Test case 2.** From left to right: reference solution (standard FEM with a fine mesh); difference between the reference solution and the projection of each approximated solution (Standard FEM, Direct $\varphi$-FEM, Dual $\varphi$-FEM).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5: **Test case 2.** $L^2$ (left) and semi-$H^1$ (right) relative errors with respect to the mesh size.</figcaption></figure>

[2] Roland Becker, Raphaël Bulle, Michel Duprez, and Vanessa Lleras. Residual-based a posteriori error estimates with boundary correction for Phi-FEM. working paper or preprint, February 2025.

[3] E. Burman, S. Claus, P. Hansbo, M. G. Larson, and A. Massing. CutFEM: Discretizing geometry and partial differential equations. *International Journal for Numerical Methods in Engineering*, 104(7):472–501, 2015.

[4] E. Burman, P. Hansbo, and M. Larson. Low regularity estimates for cutfem approximations of an elliptic problem with mixed boundary conditions. *Mathematics of Computation*, 93(345):35–54, 2024.

[5] S. Cotin, M. Duprez, V. Lleras, A. Lozinski, and K. Vuillemot. $\varphi$-FEM: An efficient simulation tool using simple meshes for problems in structure mechanics and heat transfer. *Partition of Unity Methods*, pages 191–216, 2023.

[6] M. Duprez, V. Lleras, and A. Lozinski. A new $\phi$-FEM approach for problems with natural boundary conditions. *Numer. Methods Partial Differential Equations*, 39(1):281–303, 2023.

[7] M. Duprez, V. Lleras, and A. Lozinski. $\varphi$-FEM: an optimally convergent and easily implementable immersed boundary method for particulate flows and Stokes equations. *ESAIM: Mathematical Modelling and Numerical Analysis*, 57(3):1111–1142, 2023.

[8] M. Duprez, V. Lleras, and A. Lozinski. $\varphi$-FEM: an optimally convergent and easily implementable immersed boundary method for particulate flows and Stokes equations. *ESAIM Math. Model. Numer. Anal.*, 57(3):1111–1142, 2023.

[9] M. Duprez, V. Lleras, A. Lozinski, V. Vigon, and K. Vuillemot. $\varphi$-FD : A well-conditioned finite difference method inspired by $\varphi$-FEM for general geometries on elliptic PDEs. *Journal of Scientific Computing*, 104(1):1–27, 2025.

[10] M. Duprez, V. Lleras, A. Lozinski, V. Vigon, and K. Vuillemot. $\varphi$-fem-fno: A new approach to train a neural operator as a fast pde solver for variable geometries. *Communications in Nonlinear Science and Numerical Simulation*, 152:109131, 2026.
[11] M. Duprez, V. Lleras, A. Lozinski, and K. Vuillemot. $\varphi$-FEM for the heat equation: optimal convergence on unfitted meshes in space. *Comptes Rendus. Mathématique*, 361(G11):1699-1710, 2023.

[12] M. Duprez and A. Lozinski. $\varphi$-FEM: a finite element method on domains defined by level-sets. *SIAM J. Numer. Anal.*, 58(2):1008-1028, 2020.

[13] A. Main and G. Scovazzi. The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems. *J. Comput. Phys.*, 372:972-995, 2018.

[14] R. Mittal and G. Iaccarino. Immersed boundary methods. *Annu. Rev. Fluid Mech.*, 37:239-261, 2005.