# An equilibrated flux a posteriori error estimator for defeaturing problems

Annalisa Buffa, Ondine Chanon, Denise Grappein, Rafael Vázquez, Martin Vohralík

► To cite this version:

Annalisa Buffa, Ondine Chanon, Denise Grappein, Rafael Vázquez, Martin Vohralík. An equilibrated flux a posteriori error estimator for defeaturing problems. SIAM Journal on Numerical Analysis, 2024, 62 (6), pp.2439-2458. 10.1137/23M1627195. hal-04379236

HAL Id: hal-04379236

https://inria.hal.science/hal-04379236v1

Submitted on 8 Jan 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
An equilibrated flux *a posteriori* error estimator for
defeaturing problems

Annalisa Buffa*

Ondine Chanon†

Denise Grappein‡

Rafael Vázquez§

Martin Vohralík ¶

December 26, 2023

Abstract

An *a posteriori* error estimator based on an equilibrated flux reconstruction is proposed for defeaturing problems in the context of finite element discretizations. Defeaturing consists in the simplification of a geometry by removing features that are considered not relevant for the approximation of the solution of a given PDE. In this work, the focus is on Poisson equation with Neumann boundary conditions on the feature boundary. The estimator accounts both for the so-called *defeaturing error* and for the numerical error committed by approximating the solution on the defeatured domain. Unlike other estimators that were previously proposed for defeaturing problems, the use of the equilibrated flux reconstruction allows to obtain a sharp bound for the numerical component of the error. Furthermore, it does not require the evaluation of the normal trace of the numerical flux on the feature boundary: this makes the estimator well-suited for finite element discretizations, in which the normal trace of the numerical flux is typically discontinuous across elements. The reliability of the estimator is proven and verified on several numerical examples. Its capability to identify the most relevant features is also shown, in anticipation of a future application to an adaptive strategy.

**Keywords:** Geometric defeaturing problems, *a posteriori* error estimation, equilibrated flux

**MSC codes:** 65N15, 65N30

# 1 Introduction

The need of solving problems on complex domains, characterized by the presence of geometrical features of different scales and shapes, arises in many practical applications. In particular, in the process of simulation-based manufacturing, repeated simulations are often to be performed,

<footnote>*MNS, Institute of Mathematics, École Polytechnique Fédérale de Lausanne, Switzerland & IMATI CNR - Via Ferrata 5, 27100 Pavia (annalisa.buffa@epfl.ch)</footnote>

<footnote>†Institute of Analysis and Scientific Computing, Technische Universität Wien, Austria (ondine.chanon@asc.tuwien.ac.at)</footnote>

<footnote>‡Dipartimento di Scienze Matematiche G. L. Lagrange, Politecnico di Torino, Italy, Member of GNCS INdAM Group (denise.grappein@polito.it)</footnote>

<footnote>§Departamento de Matemática Aplicada, Universidade de Santiago de Compostela, Spain (rafael.vazquez@usc.es)</footnote>

<footnote>¶Inria, 2 rue Simone Iff, 75589 Paris, France & CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée, France (martin.vohralik@inria.fr)</footnote>
in order to analyze the impact of design changes or to adjust geometric parameters. In many
cases, before even solving the problem at hand, the first issue to overcome is the definition of the
features themselves and the construction of a suitable computational mesh. For this reason it can
be fundamental to simplify the geometry as much as possible, in order to avoid the definition of
those features which may not have an actual impact on the accuracy of the solution. This process is
commonly called *defeaturing*. Some criteria based on some *a priori* knowledge of the computational
domain and of the properties of the materials have been used in the past (see, e.g., [10, 11, 22]).
However, in order to fully automatize the process, an *a posteriori* criterion is necessary and many
different proposals can be found in literature (see [6, 9, 12–17, 20, 21, 23]).

In this paper we start from the work presented in [4, 5], which proposes an *a posteriori* error
estimator for analysis-aware defeaturing, in the context of the Poisson equation with Neumann
boundary conditions on the feature boundary. In particular in [4], an estimator is designed to
control the overall error between the exact solution of the PDE defined in the exact domain, and
the numerical approximation of the solution of the corresponding PDE defined in the defeatured
domain. This estimator is made by two components, one accounting for the defeaturing error, i.e.
the error committed by neglecting the features, and the other accounting for the numerical error
committed when solving the problem on the defeatured geometry. The first component has the big
advantage of being explicit with respect to the size of the geometrical features, and in [5] the authors
prove that it is a reliable and efficient bound for the energy norm of the defeaturing error. The
second component is instead built as a residual-based estimator of the numerical error. The overall
estimator is defined up to two positive parameters, related to the unknown constants appearing in
the bounds of the defeaturing and of the numerical errors. Such parameters need to be tuned in
order to correctly weight the two components.

In order to partially overcome this issue, in this work we propose a novel *a posteriori* error
estimator that is strongly based on [5] for what concerns its defeaturing component, but which
resorts to an *equilibrated flux reconstruction* (see, among others, [1, 3, 7, 18]). Indeed, one of
the main drawbacks of residual-based error estimators is that the reliability constants are usually
unknown and problem dependent. On the contrary, the difference between the numerical and
the equilibrated flux provides an upper bound for the energy norm of the numerical error having
reliability constant equal to 1. Although we do not get rid of the unknown constant related to
the defeaturing component, the use of the equilibrated flux reconstruction also allows to avoid the
computation of the normal trace of the numerical flux on the feature boundary. This makes the
estimator well-suited for finite element discretizations, in which the normal trace of the numerical
flux is typically not continuous. On the contrary, the estimator proposed in [4] was designed to be
applied along with an IGA discretization.

The equilibrated flux reconstruction is built following the steps in [3, 8], solving mixed local problems on patches of elements and leading to a discrete reconstructed flux in a Raviart—Thomas finite element space.

The paper is organized in four Sections. In Section 2 we introduce some notation and the
defeaturing model problem, while in Section 3 and 4 we derive and analyze an *a posteriori* error
estimator resorting to a generic equilibrated flux reconstruction and providing a bound for the
overall error. Section 5 describes a practical way to build the equilibrated flux reconstruction and,
finally, in Section 6 the proposed estimator is validated by some numerical experiments.
<figure><img src="image_2.png" /><figcaption>(a) Domain with negative feature. (b) Domain with positive feature. (c) Defeatured geometry.</figcaption></figure>

<caption>Figure 1: Domain with a negative feature, domain with a positive feature and corresponding defeatured domain for both configurations.</caption>

## 2 Notation and model problem

In the following we adopt the notation introduced in [5], which is here recalled for the sake of clarity. Let $\omega$ be any open $k$-dimensional manifold in $\mathbb{R}^d$, $d = 2, 3$ and $k \le d$. We denote by $|\omega|$ the measure of $\omega$, and for any function $\varphi$ defined on $\omega$, we denote by $\bar{\varphi}^\omega$ its average over $\omega$. We will denote by $(\cdot, \cdot)_\omega$ the $L^2$-inner product on $\omega$ and by $||\cdot||_\omega$ the corresponding norm. If $k < d$, then $\langle \cdot, \cdot \rangle_\omega$ stands for a duality paring on $\omega$. For future use, let us define the quantity

$$c_{\omega} := \begin{cases} \max(-\log(|\omega|), \zeta)^{\frac{1}{2}} & \text{if } k=1, d=2 \\ 1 & \text{if } k=2, d=3 \end{cases} \qquad (1)$$

where $\zeta \in \mathbb{R}$ is the unique solution of $\zeta = -\log(\zeta)$.

Let us consider an open Lipschitz domain $\Omega \subset \mathbb{R}^d$ and let us denote by $\partial\Omega$ its boundary. We suppose that $\Omega$ contains one feature $F \subset \mathbb{R}^d$, i.e. a geometrical detail of smaller scale, which is assumed to be an open Lipschitz domain, as well. The boundary of $F$ is denoted by $\partial F$. We consider two main types of features. In particular, a feature $F$ is said to be

* negative, if $(\bar{F} \cap \bar{\Omega}) \subset \partial\Omega$;

* positive, if $F \subset \Omega$.

In the following we will refer to $\Omega$ as the *exact* or *original* geometry. For the sake of simplicity we restrict ourselves to the case of an exact geometry with a single feature, but the generalization to the multiple feature case easily follows from [2, 4].

Let us now define the so called *defeatured geometry*, i.e. $\Omega_0 \subset \mathbb{R}^d$ such that

$$\Omega_0 := \begin{cases} \text{int}(\bar{\Omega} \cup \bar{F}) & \text{if } F \text{ is negative} \\ \Omega \setminus \bar{F} & \text{if } F \text{ is positive.} \end{cases}$$

Hence, if the feature is negative, $\Omega \subset \Omega_0$ (Figure 1a), while if the feature is positive, $\Omega_0 \subset \Omega$ (Figure 1b). In the following, the boundary of $\Omega_0$ is denoted by $\partial\Omega_0$.

We denote by $n$, $n_0$ and $n_F$ the unitary outward normals respectively of $\Omega$, $\Omega_0$ and $F$. Let $\partial\Omega = \bar{\Gamma}_D \cup \bar{\Gamma}_N$, with $\Gamma_D \cap \Gamma_N = \emptyset$ and $\Gamma_D \neq \emptyset$, and we assume that $\partial F \cap \Gamma_D = \emptyset$. Let $\gamma_0 := \partial F \setminus \bar{\Gamma}_N \subset \partial\Omega_0$ and, finally, let $\gamma := \partial F \setminus \bar{\gamma}_0 \subset \partial\Omega$, so that $\partial F = \bar{\gamma}_0 \cup \bar{\gamma}$ and $\gamma_0 \cap \gamma = \emptyset$. Let us
observe that, if $\gamma_0 = \emptyset$, then we are in the case of a negative internal feature, i.e. $\Omega$ is a perforated domain (see, for an example, Figures 3 and 11 in Section 6).

On the exact geometry $\Omega$ we use the Poisson problem as a model problem:

$$ egin{cases} -\Delta u = f & \text{in } \Omega \\ u = g_D & \text{on } \Gamma_D \\ \nabla u \cdot \mathbf{n} = g & \text{on } \Gamma_N, \end{cases} \qquad (2) $$

to which we will also refer as the original problem. Defining

$$ H_{0,\Gamma_D}^1(\Omega) = \{v \in H^1(\Omega) : v|_{\Gamma_D} = 0\}, \quad H_{g_D,\Gamma_D}^1(\Omega) = \{v \in H^1(\Omega) : v|_{\Gamma_D} = g_D\}, $$

the variational formulation of Problem (2) reads: find $u \in H_{g_D,\Gamma_D}^1(\Omega)$ which satisfies

$$ (\nabla u, \nabla v)_{\Omega} = (f, v)_{\Omega} + \langle g, v \rangle_{\Gamma_N} \quad \forall v \in H_{0,\Gamma_D}^1(\Omega). \qquad (3) $$

On the defeatured geometry $\Omega_0$ we consider instead the problem

$$ egin{cases} -\Delta u_0 = f & \text{in } \Omega_0 \\ u_0 = g_D & \text{on } \Gamma_D \\ \nabla u_0 \cdot \mathbf{n}_0 = g & \text{on } \Gamma_N \setminus \gamma \\ \nabla u_0 \cdot \mathbf{n}_0 = g_0 & \text{on } \gamma_0 \end{cases} \qquad (4) $$

to which we will also refer to as defeatured problem. With an abuse of notation, in the negative feature case, we denote by $f \in L^2(\Omega_0)$ a suitable $L^2$-extension of $f \in L^2(\Omega)$ to $F$, while the Neumann datum $g_0$ has to be chosen. The variational formulation of problem (4) reads: find $u_0 \in H_{g_D,\Gamma_D}^1(\Omega_0)$ which satisfies, $\forall v \in H_{0,\Gamma_D}^1(\Omega_0)$

$$ (\nabla u_0, \nabla v)_{\Omega_0} = (f, v)_{\Omega_0} + \langle g, v \rangle_{\Gamma_N \setminus \gamma} + \langle g_0, v \rangle_{\gamma_0}. \qquad (5) $$

Let us consider a partition $\mathcal{T}_h$ of $\Omega_0$ consisting of closed triangles $K$ for $d=2$, or tetrahedrons for $d=3$, such that $\overline{\Omega}_0 = \bigcup_{K \in \mathcal{T}_h} K$. Hereby, we suppose that the mesh faces match with the boundaries $\Gamma_D, \Gamma_N \setminus \gamma$ and $\gamma_0$. Let us then introduce the set

$$ Q_h = \mathcal{P}_p(\mathcal{T}_h) := \{q_h \in L^2(\Omega_0) : q_h|_K \in \mathcal{P}_p(K), \forall K \in \mathcal{T}_h\}, \qquad (6) $$

with $\mathcal{P}_p(K)$ denoting the set of polynomials of degree at most $p \ge 1$ on $K \in \mathcal{T}_h$, and

$$ V_h^0 = \{q_h \in C^0(\overline{\Omega}_0) \cap Q_h : q_h|_{\Gamma_D} = 0\}, \quad V_h := \{q_h \in C^0(\overline{\Omega}_0) : q_h|_{\Gamma_D} = g_D\}. $$

In the following, for the sake of simplicity, we assume $f \in Q_h$. Similarly, let us consider the partition of $\partial\Omega_0$ induced by the elements of $\mathcal{T}_h$ and let us denote its restriction to $(\Gamma_N \setminus \gamma) \cup \gamma_0$ by $\partial\Omega_{0,h}^N$. Introducing

$$ g_N = \begin{cases} g & \text{on } \Gamma_N \setminus \gamma \\ g_0 & \text{on } \gamma_0 \end{cases} $$

we assume $g_N$ to be an element of the broken space $\mathcal{P}_p(\partial\Omega_{0,h}^N)$, defined in the same manner as (6). Hence, the finite element approximation of (5) reads as: find $u_0^h \in V_h$ which satisfies, $\forall v_h \in V_h^0$

$$ (\nabla u_0^h, \nabla v_h)_{\Omega_0} = (f, v_h)_{\Omega_0} + \langle g, v_h \rangle_{\Gamma_N \setminus \gamma} + \langle g_0, v_h \rangle_{\gamma_0}. \qquad (7) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Example of geometry with a positive feature F and extension of F to $\tilde{F}$.</figcaption></figure>

Let us remark that our aim is to never solve Problem (2), but to design a proper a posteriori error estimator capable to control the energy norm of the error committed by approximating the exact solution of (2) by $u_0^h$. We will refer to this error as the overall error, as it accounts both for the error introduced by defeaturing and for the error introduced by the numerical approximation of $u_0$. In particular, we aim at designing an estimator based on an equilibrated flux reconstruction, which has the advantage of bounding the numerical error with a sharp reliability constant equal to 1. The flux reconstruction will be used to bound also the defeaturing error even if, in this case, we will not get rid of the unknown constant. In the following we provide the definition of the overall error for the negative and positive feature case, referring again to [5].

**Negative feature:** in this case $\Omega \subset \Omega_0$, hence we restrict $u_0$ to $\Omega$ and we define the overall error as $||\nabla(u - u_0^h)|_\Omega||_\Omega$.

**Positive feature:** this case is slightly more complicated, since $u_0$ and its finite element approximation are defined only on $\Omega_0$ and $\Omega_0 \subset \Omega$. Hence, in order to define the overall error, we need to extend $u_0$ to the feature F.

However meshing F and solving a problem on it may be non trivial, in particular if F has a complex boundary. Hence we follow the steps in [5]: we consider a suitable extension $\tilde{F} \subset \mathbb{R}^d$ of F, being as simple as possible, in particular $F \subset \tilde{F}$ and $\gamma_0 \subset (\partial\tilde{F} \cap \partial F)$, as reported in Figure 2. Let $\gamma$ be decomposed as $\gamma = \text{int}(\bar{\gamma}_s \cup \bar{\gamma}_r)$ where $\gamma_s = \gamma \cap \partial\tilde{F}$ is the portion of $\gamma$ shared by $\partial F$ and $\partial\tilde{F}$ and $\gamma_r = \gamma \setminus \bar{\gamma}_s$ (Figure 2b). We denote by $\tilde{n}$ the unitary outward normal of $\tilde{F}$ and we set $\tilde{\gamma} = \partial\tilde{F} \setminus \partial F$. On $\tilde{F}$ we solve the problem

$$ \left\{ \begin{array}{ll} -\Delta \tilde{u}_0 = f & \text{in } \tilde{F} \\ \tilde{u}_0 = u_0 & \text{on } \gamma_0 \\ \nabla \tilde{u}_0 \cdot \tilde{n} = \tilde{g} & \text{on } \tilde{\gamma} \\ \nabla \tilde{u}_0 \cdot \tilde{n} = g & \text{on } \gamma_s \end{array} \right. \qquad (8) $$

where, with an abuse of notation, we still denote by $f$ any $L^2$-extension of the forcing term to $\tilde{F}$ and the Neumann datum $\tilde{g}$ on $\tilde{\gamma}$ has to be chosen. Introducing

$$ H_{u_0, \gamma_0}^1(\tilde{F}) = \{ v \in H^1(\tilde{F}) : v|_{\gamma_0} = u_0|_{\gamma_0} \}, $$
the variational formulation of (8) is: find $\tilde{u}_0 \in H_{u_0, \gamma_0}^1(\tilde{F})$ which satisfies, $\forall v \in H_{0, \gamma_0}^1(\tilde{F})$

$$ (\nabla \tilde{u}_0, \nabla v)_{\tilde{F}} = (f, v)_{\tilde{F}} + \langle \tilde{g}, v \rangle_{\tilde{\gamma}} + \langle g, v \rangle_{\gamma_s}. \qquad (9) $$

We denote by $\tilde{u}_0^h$ the finite element approximation of $\tilde{u}_0$ on a partition $\tilde{T}_h$ of $\tilde{F}$. Note that this partition does not need to be conforming to $\gamma$. We suppose, however, that the mesh faces match with $\gamma_0$, $\tilde{\gamma}$ and $\gamma_s$ and that $\tilde{T}_h$ matches with $\tilde{T}_h$ on $\gamma_0$. Let

$$ \tilde{Q}_h = \mathcal{P}_p(\tilde{T}_h) := \{ q_h \in L^2(\tilde{F}) : q_h|_{K} \in \mathcal{P}_p(K), \forall K \in \tilde{T}_h \}, \qquad (10) $$

and let us introduce

$$ \tilde{V}_h^0 = \{ q_h \in C^0(\tilde{F}) \cap \tilde{Q}_h : q_h|_{\gamma_0} = 0 \}, \quad \tilde{V}_h := \{ q_h \in C^0(\tilde{F}) \cap \tilde{Q}_h : q_h|_{\gamma_0} = u_0^h|_{\gamma_0} \}. $$

We assume for simplicity that $f|_{\tilde{F}} \in \tilde{Q}_h$. Considering the partition of $\partial \tilde{F}$ induced by the elements of $\tilde{T}_h$ and denoting its restriction to $\gamma_s \cup \tilde{\gamma}$ as $\partial \tilde{F}_h^N$, we also assume that

$$ \tilde{g}_N = \begin{cases} g & \text{on } \gamma_s \\ \tilde{g} & \text{on } \tilde{\gamma} \end{cases} $$

is an element of the broken space $\mathcal{P}_p(\partial \tilde{F}_h^N)$. The finite element approximation of Problem (9) is hence: find $\tilde{u}_0^h \in \tilde{V}_h$ which satisfies, $\forall v \in \tilde{V}_h^0$

$$ (\nabla \tilde{u}_0^h, \nabla v_h)_{\tilde{F}} = (f, v_h)_{\tilde{F}} + \langle \tilde{g}, v_h \rangle_{\tilde{\gamma}} + \langle g, v_h \rangle_{\gamma_s}. \qquad (11) $$

Finally, we define the extended defeatured solution and its numerical approximation as

$$ u_d := \begin{cases} u_0 & \text{in } \Omega_0 \\ \tilde{u}_0 & \text{in } \tilde{F}, \end{cases} \quad u_d^h := \begin{cases} u_0^h & \text{in } \tilde{T}_h \\ \tilde{u}_0^h & \text{in } \tilde{T}_h, \end{cases} \qquad (12) $$

while the overall error is $||\nabla(u - u_d^h)||_\Omega$.

### 3 Negative feature a posteriori error estimator

In this section we propose a reliable estimator for the overall error $||\nabla(u - u_0^h|_\Omega)||_\Omega$, in the case of a single negative feature. To simplify the notation, in the following we omit the explicit restriction of $u_0^h$ (and $u_0$) to $\Omega$.

Let us consider the solution to problem (5): introducing the flux $\sigma = -\nabla u_0$, we have that $\sigma \in \mathcal{H}(\text{div}, \Omega_0)$, $\nabla \cdot \sigma = f$, $-\sigma \cdot n_0 = g$ on $\Gamma_N \setminus \gamma$ and $-\sigma \cdot n_0 = g_0$ on $\gamma_0$. At discrete level, a suitable definition of flux is more involved. Indeed $\nabla u_0^h \notin \mathcal{H}(\text{div}, \Omega_0)$ and hence the divergence equation and the Neumann boundary condition are not exactly satisfied. The idea behind the equilibrated flux reconstruction is to use the discrete solution $u_0^h$ to build a discrete flux $\sigma_h$ such that $\sigma_h \in \mathcal{H}(\text{div}, \Omega_0)$ is an approximation of $\sigma$ satisfying

$$ \left\{ \begin{array}{ll} \nabla \cdot \sigma_h = f & \text{in } \Omega_0 \\ \sigma_h \cdot n_0 = -g & \text{on } \Gamma_N \setminus \gamma \\ \sigma_h \cdot n_0 = -g_0 & \text{on } \gamma_0. \end{array} \right. \qquad (13) $$
We will give more details about how an equilibrated flux reconstruction can actually be computed in Section 5, following [3, 8]. For the time being we assume that we have $\sigma_h$ at our disposal.

Referring to the notation introduced in Section 2, let us introduce, on $\gamma$, the quantity

$$d_{\gamma}^{h} := g + \boldsymbol{\sigma}_{h} \cdot \mathbf{n} \quad \text{on } \gamma$$

which is the error between the Neumann datum $g$ on $\gamma$ and the normal trace of the equilibrated flux reconstruction. Following [5], denoting by $\overline{d}_{\gamma}^{h}$ the average of $d_{\gamma}^{h}$ over $\gamma$, let us define

$$\mathcal{E}_{\gamma} := \left( |\gamma|^{\frac{1}{d-1}} \left\| d_{\gamma}^{h} - \overline{d}_{\gamma}^{h} \right\|_{\gamma}^{2} + c_{\gamma}^{2} |\gamma|^{\frac{d}{d-1}} \left\| \overline{d}_{\gamma}^{h} \right\|^{2} \right)^{\frac{1}{2}}, \quad (14)$$

$$\mathcal{E}_0 := ||\boldsymbol{\sigma}_h + \nabla u_0^h||_{\Omega_0} \quad (15)$$

where $c_{\gamma}$ is defined as in (1). Let us remark that, unlike [5], the quantity $d_{\gamma}^{h}$ does not depend on the normal trace of the numerical flux, but on the normal trace of the equilibrated flux reconstruction, which is continuous across the elements of the mesh $\mathcal{T}_h$. The following proposition establishes our a posteriori bound:

**Proposition 1.** Let $u$ be the solution of (3) and $u_0^h$ the solution of (7). Then

$$||\nabla(u - u_0^h)||_{\Omega} \le C_D \mathcal{E}_{\gamma} + \mathcal{E}_0, \quad (16)$$

with $C_D > 0$ being a constant independent of the size of feature $F$.

*Proof.* Let $v \in H_{0,\Gamma_D}^1(\Omega)$. Adding and subtracting $(\boldsymbol{\sigma}_h, \nabla v)_{\Omega}$, exploiting (3), applying Green's theorem and using the characterization of $\boldsymbol{\sigma}_h$ provided in (13), we have

$$\begin{aligned}
(\nabla(u - u_0^h), \nabla v)_{\Omega} &= (\nabla u + \boldsymbol{\sigma}_h, \nabla v)_{\Omega} - (\boldsymbol{\sigma}_h + \nabla u_0^h, \nabla v)_{\Omega} \\
&= (f - \nabla \cdot \boldsymbol{\sigma}_h, v)_{\Omega} + \langle g + \boldsymbol{\sigma}_h \cdot \mathbf{n}, v \rangle_{\Gamma_N} - (\boldsymbol{\sigma}_h + \nabla u_0^h, \nabla v)_{\Omega} \\
&= \langle g + \boldsymbol{\sigma}_h \cdot \mathbf{n}, v \rangle_{\gamma} - (\boldsymbol{\sigma}_h + \nabla u_0^h, \nabla v)_{\Omega} \\
&= \langle d_{\gamma}^h, v \rangle_{\gamma} - (\boldsymbol{\sigma}_h + \nabla u_0^h, \nabla v)_{\Omega}.
\end{aligned} \quad (17)$$

Referring the reader to the steps reported in [5], with the difference that the numerical flux is here substituted by the equilibrated flux reconstruction, it is possible to prove that

$$\langle d_{\gamma}^{h}, v \rangle_{\gamma} \le C_D \mathcal{E}_{\gamma} ||\nabla v||_{\Omega} \quad (18)$$

with $C_D > 0$ being a constant independent of the size of feature $F$ (see Theorem 4.3 in [5]). If we choose $v = u - u_0^h \in H_{0,\Gamma_D}^1(\Omega)$ in (17), we apply (18) and the Cauchy-Schwarz inequality, we have

$$\begin{aligned}
||\nabla(u - u_0^h)||_{\Omega}^2 &\le C_D \mathcal{E}_{\gamma} ||\nabla(u - u_0^h)||_{\Omega} + ||\boldsymbol{\sigma}_h + \nabla u_0^h||_{\Omega} ||\nabla(u - u_0^h)||_{\Omega} \\
&\le C_D \mathcal{E}_{\gamma} ||\nabla(u - u_0^h)||_{\Omega} + ||\boldsymbol{\sigma}_h + \nabla u_0^h||_{\Omega_0} ||\nabla(u - u_0^h)||_{\Omega} \\
&= (C_D \mathcal{E}_{\gamma} + \mathcal{E}_0) ||\nabla(u - u_0^h)||_{\Omega},
\end{aligned}$$

where we have also exploited the fact that, in the negative feature case, $\Omega \subset \Omega_0$. Simplifying on both sides yields (16). $\square$
**Remark 1.** It is well known from literature (see, among others, [1, 3, 7, 18]) that the quantity $\mathcal{E}_0$ provides a sharp upper bound for the numerical error $||\nabla(u_0 - u_0^h)||_{\Omega_0}$. Let us remark that, if no feature is present, the same result is provided also by (16). Indeed, if $\gamma = \emptyset$, then $u = u_0$, $\Omega = \Omega_0$ and (16) reduces to

$$||\nabla(u_0 - u_0^h)||_{\Omega_0} \le ||\boldsymbol{\sigma}_h + \nabla u_0^h||_{\Omega_0}.$$

For this reason we will refer to $\mathcal{E}_0$ as the numerical component of the estimator and to $\mathcal{E}_\gamma$ as the defecturing component.

# 4 Positive feature a posteriori error estimator

In this section we propose a reliable estimator for the overall error $||\nabla(u - u_d^h)||_{\Omega}$, in the case of a single positive feature $F$. For the sake of generality, we consider the case in which $F$ is embedded in a smooth extension $\tilde{F}$, as detailed in Section 2.

Let us introduce an equilibrated flux reconstruction on $\tilde{F}$, i.e. a discrete flux $\tilde{\boldsymbol{\sigma}}_h \in \mathcal{H}(\text{div}, \tilde{F})$ built somehow from $\tilde{u}_0^h$ such that

$$\left\{ \begin{array}{ll} \nabla \cdot \tilde{\boldsymbol{\sigma}}_h = f & \text{in } \tilde{F} \\ \tilde{\boldsymbol{\sigma}}_h \cdot \tilde{n} = -\tilde{g} & \text{on } \tilde{\gamma} \\ \tilde{\boldsymbol{\sigma}}_h \cdot \tilde{n} = -g & \text{on } \gamma_s. \end{array} \right. \qquad (19)$$

Again the details on the construction of this flux will be provided in Section 5 and, for the time being, we assume we have $\tilde{\boldsymbol{\sigma}}_h$. In this case we define on $\gamma_0$ the quantity

$$d_{\gamma_0}^h := \tilde{\boldsymbol{\sigma}}_h \cdot \boldsymbol{n}_F - g_0 \quad \text{on } \gamma_0$$

which approximates the jump in the normal derivative of $u_d$ on $\gamma_0$, while on $\gamma_r$ we define

$$d_{\gamma_r}^h := \tilde{\boldsymbol{\sigma}}_h \cdot \boldsymbol{n}_F + g \quad \text{on } \gamma_r$$

which is the error between the Neumann datum $g$ on $\gamma_r$ and the normal trace of the equilibrated flux reconstruction computed on $\tilde{F}$. Again we observe how, unlike [5], the normal trace of the numerical flux is not involved in the definition of these quantities. Denoting by $d_{\gamma_0}^h$ the average of $d_{\gamma_0}^h$ on $\gamma_0$ and by $d_{\gamma_r}^{h, \gamma_r}$ the average of $d_{\gamma_r}^h$ on $\gamma_r$ and following [5], let us introduce

$$\tilde{\mathcal{E}}_{\gamma_0} := \left( |\gamma_0|^{\frac{1}{d-1}} \left\| d_{\gamma_0}^h - \overline{d_{\gamma_0}^h}_{\gamma_0} \right\|_{\gamma_0}^2 + c_{\gamma_0}^2 |\gamma_0|^{\frac{d}{d-1}} \overline{d_{\gamma_0}^h}_{\gamma_0} \right)^{\frac{1}{2}} \qquad (20)$$

$$\tilde{\mathcal{E}}_{\gamma_r} := \left( |\gamma_r|^{\frac{1}{d-1}} \left\| d_{\gamma_r}^h - \overline{d_{\gamma_r}^h}_{\gamma_r} \right\|_{\gamma_r}^2 + c_{\gamma_r}^2 |\gamma_r|^{\frac{d}{d-1}} \overline{d_{\gamma_r}^h}_{\gamma_r} \right)^{\frac{1}{2}}, \qquad (21)$$

where $c_{\gamma_0}$ and $c_{\gamma_r}$ are defined as in (1). Let us also define

$$\tilde{\mathcal{E}}_0 := ||\tilde{\boldsymbol{\sigma}}_h + \nabla \tilde{u}_0^h||_{\tilde{F}}, \qquad (22)$$

and let us recall that $\mathcal{E}_0 = ||\boldsymbol{\sigma}_h + \nabla u_0^h||_{\Omega_0}$, where $\boldsymbol{\sigma}_h$ is an equilibrated flux reconstruction defined in $\Omega_0$ as in (13).
**Proposition 2.** Let $u$ be the solution of (3) and $u_d^h$ be defined as in (12). Then

$$ ||\nabla(u - u_d^h)||_{\Omega} \le \tilde{C}_D (\tilde{\mathcal{E}}_{\gamma_0}^2 + \tilde{\mathcal{E}}_{\gamma_r}^2)^{\frac{1}{2}} + (\tilde{\mathcal{E}}_0^2 + \mathcal{E}_0^2)^{\frac{1}{2}} \quad (23) $$

with $\tilde{C}_D$ being a constant independent of the size of feature F.

*Proof.* Let us consider the restriction to $\Omega_0$ of the solution $u$ of problem (2), which satisfies

$$ \left\{ \begin{array}{ll} -\Delta u|_{\Omega_0} = f & \text{in } \Omega_0 \\ u|_{\Omega_0} = g_D & \text{on } \Gamma_D \\ \nabla u|_{\Omega_0} \cdot \mathbf{n} = g & \text{on } \Gamma_N \setminus \gamma \\ \nabla u|_{\Omega_0} \cdot \mathbf{n}_0 = \nabla u \cdot \mathbf{n}_0 & \text{on } \gamma_0. \end{array} \right. \qquad (24) $$

Omitting the explicit restriction on $u$ to $\Omega_0$ the variational formulation of problem (24) reads: find $u \in H_{g_D, \Gamma_D}^1(\Omega)$ which satisfies, $\forall v \in H_{0, \Gamma_D}^1(\Omega)$

$$ (\nabla u, \nabla v)_{\Omega_0} = (f, v)_{\Omega_0} + \langle g, v \rangle_{\Gamma_N \setminus \gamma} + \langle \nabla u \cdot \mathbf{n}_0, v \rangle_{\gamma_0}. \quad (25) $$

Let $v \in H_{0, \Gamma_D}^1(\Omega_0)$. Adding and subtracting $(\sigma_h, \nabla v)_{\Omega_0}$, exploiting (25), applying Green's theorem and using the characterization of $\sigma_h$ provided in (13), we have

$$ \begin{aligned} (\nabla(u - u_0^h), \nabla v)_{\Omega_0} &= (\nabla u + \sigma_h, \nabla v)_{\Omega_0} - (\sigma_h + \nabla u_0^h, \nabla v)_{\Omega_0} \\ &= (f - \nabla \cdot \sigma_h, v)_{\Omega_0} + \langle g + \sigma_h \cdot \mathbf{n}, v \rangle_{\Gamma_N \setminus \gamma} \\ &\quad + \langle \nabla u \cdot \mathbf{n}_0 + \sigma_h \cdot \mathbf{n}_0, v \rangle_{\gamma_0} - (\sigma_h + \nabla u_0^h, \nabla v)_{\Omega_0} \\ &= \langle \nabla u \cdot \mathbf{n}_0 - g_0, v \rangle_{\gamma_0} - (\sigma_h + \nabla u_0^h, \nabla v)_{\Omega_0}. \end{aligned} \quad (26) $$

In order to obtain an actual error indicator we need to estimate the quantity $\langle \nabla u \cdot \mathbf{n}_0 - g_0, v \rangle_{\gamma_0}$ on the right hand side, and for this reason we must consider the error committed on the feature as well. Hence, let us consider the restriction to the positive feature $F$ of the solution $u$ of (2), satisfying

$$ \left\{ \begin{array}{ll} -\Delta u|_F = f & \text{in } F \\ \nabla u|_F \cdot \mathbf{n}_F = g & \text{on } \gamma \\ \nabla u|_F \cdot \mathbf{n}_F = \nabla u \cdot \mathbf{n}_F & \text{on } \gamma_0, \end{array} \right. \qquad (27) $$

so that, omitting the explicit restriction of $u$ to $F$, $u \in H^1(F)$ is one of the infinitely many solutions, defined up to a constant, of

$$ (\nabla u, \nabla v)_F = (f, v)_F + \langle g, v \rangle_\gamma + \langle \nabla u \cdot \mathbf{n}_F, v \rangle_{\gamma_0} \quad \forall v \in H^1(F). \quad (28) $$

Let $v \in H^1(F)$. Adding and subtracting $(\tilde{\sigma}_h, \nabla v)_F$, exploiting (28), applying Green's theorem and using the characterization of $\tilde{\sigma}_h$ provided in (19), we have

$$ \begin{aligned} (\nabla(u - \tilde{u}_0^h), \nabla v)_F &= (\nabla u + \tilde{\sigma}_h, \nabla v)_F - (\tilde{\sigma}_h + \nabla \tilde{u}_0^h, \nabla v)_F \\ &= (f - \nabla \cdot \tilde{\sigma}_h, v)_F + \langle g + \tilde{\sigma}_h \cdot \mathbf{n}_F, v \rangle_{\gamma_r \cup \gamma_s} \\ &\quad + \langle \nabla u \cdot \mathbf{n}_F + \tilde{\sigma}_h \cdot \mathbf{n}_F, v \rangle_{\gamma_0} - (\tilde{\sigma}_h + \nabla \tilde{u}_0^h, \nabla v)_F \\ &= \langle g + \tilde{\sigma}_h \cdot \mathbf{n}_F, v \rangle_{\gamma_r} + \langle \nabla u \cdot \mathbf{n}_F + \tilde{\sigma}_h \cdot \mathbf{n}_F, v \rangle_{\gamma_0} - (\tilde{\sigma}_h + \nabla \tilde{u}_0^h, \nabla v)_F. \end{aligned} \quad (29) $$
Choosing $v \in H_{0,\Gamma_D}(\Omega)$, observing that $v|_{\Omega_0} \in H_{0,\Gamma_D}^1(\Omega_0)$ and $v|_F \in H^1(F)$, and summing (26) and (29), we obtain

$$ 
\begin{aligned}
(\nabla(u - u_d^h), \nabla v)_\Omega &= (\nabla(u - u_0^h), \nabla v)_{\Omega_0} + (\nabla(u - \tilde{u}_0^h), \nabla v)_F \\
&= \langle g + \tilde{\sigma}_h \cdot n_F, v \rangle_{\gamma_r} + \langle \tilde{\sigma}_h \cdot n_F - g_0, v \rangle_{\gamma_0} \\
&\quad - (\sigma_h + \nabla u_0^h, \nabla v)_{\Omega_0} - (\tilde{\sigma}_h + \nabla \tilde{u}_0^h, \nabla v)_F \\
&= \langle d_{\gamma_r}^h, v \rangle_{\gamma_r} + \langle d_{\gamma_0}^h, v \rangle_{\gamma_0} - (\sigma_h + \nabla u_0^h, \nabla v)_{\Omega_0} - (\tilde{\sigma}_h + \nabla \tilde{u}_0^h, \nabla v)_F,
\end{aligned} 
\quad (30) $$

where we have used that $n_F = -n_0$ on $\gamma_0$.

For the terms involving $d_{\gamma_r}^h$ and $d_{\gamma_0}^h$ we proceed similarly to the negative feature case: referring the reader to Theorem 5.5 in [5], it is possible to prove that there exists a constant $\tilde{C}_D > 0$ such that

$$ \langle d_{\gamma_0}^h, v \rangle_{\gamma_r} + \langle d_{\gamma_r}^h, v \rangle_{\gamma_0} \le \tilde{C}_D (\tilde{\mathcal{E}}_{\gamma_r}^2 + \tilde{\mathcal{E}}_{\gamma_0}^2)^{\frac{1}{2}} ||\nabla v||_{\Omega}. \quad (31) $$

If we choose $v = u - u_d^h \in H_{0,\Gamma_D}(\Omega)$ in (30), we use (31) and the Cauchy–Schwarz inequality we obtain

$$ 
\begin{aligned}
||\nabla(u - u_d^h)||_{\Omega}^2 &\le \tilde{C}_D (\tilde{\mathcal{E}}_{\gamma_r}^2 + \tilde{\mathcal{E}}_{\gamma_0}^2)^{\frac{1}{2}} ||\nabla(u - u_d^h)||_{\Omega} + ||\sigma_h + \nabla u_0^h||_{\Omega_0} ||\nabla(u - u_d^h)||_{\Omega_0} \\
&\quad + ||\tilde{\sigma}_h + \nabla \tilde{u}_0^h||_F ||\nabla(u - u_d^h)||_F \\
&\le \tilde{C}_D (\tilde{\mathcal{E}}_{\gamma_r}^2 + \tilde{\mathcal{E}}_{\gamma_0}^2)^{\frac{1}{2}} ||\nabla(u - u_d^h)||_{\Omega} \\
&\quad + (||\sigma_h + \nabla u_0^h||_{\Omega_0}^2 + ||\tilde{\sigma}_h + \nabla \tilde{u}_0^h||_F^2)^{\frac{1}{2}} (||\nabla(u - u_d^h)||_{\Omega_0}^2 + ||\nabla(u - u_d^h)||_F^2)^{\frac{1}{2}} \\
&\le (\tilde{C}_D (\tilde{\mathcal{E}}_{\gamma_r}^2 + \tilde{\mathcal{E}}_{\gamma_0}^2)^{\frac{1}{2}} + (\mathcal{E}_0^2 + \tilde{\mathcal{E}}_0^2)^{\frac{1}{2}}) ||\nabla(u - u_d^h)||_{\Omega},
\end{aligned} 
$$

where, in the last step, we have exploited the fact that $F \subseteq \tilde{F}$. The thesis follows by simplifying on both sides. $\square$

**Remark 2.** If the feature $F$ is relatively simple, there is no need to use an extension, and problem (9) is solved directly in $F$. In this case, maintaining the tilde notation (since $F = \tilde{F}$), expression (23) simplifies into

$$ ||\nabla(u - u_d^h)||_{\Omega} \le \tilde{C}_D \tilde{\mathcal{E}}_{\gamma_0} + (\tilde{\mathcal{E}}_0^2 + \mathcal{E}_0^2)^{\frac{1}{2}}. \quad (32) $$

# 5 Equilibrated fluxes reconstruction

In this section we describe how to build, in practice, an equilibrated flux starting from the discrete solution of the defeatured problem $u_0^h$ or from $u_d^h$ in the positive feature case. The proposed procedure is directly adapted from [3, 8] and resorts to a local reconstruction of the fluxes. Given the triangular/tetrahedral mesh $\mathcal{T}_h$ built on the defeatured geometry $\Omega_0$, let us denote by $\mathcal{N}_h$ the set of its vertices and let us divide it into interior vertices $\mathcal{N}_h^{\text{int}}$ and boundary vertices $\mathcal{N}_h^{\text{ext}}$. We aim at reconstructing the flux in the Raviart–Thomas space of order $p \ge 1$, namely in

$$ M_h := \{v_h \in \mathbf{H}(\text{div}, \Omega_0) : v_h|_K \in [\mathcal{P}_p(K)]^d + x\mathcal{P}_p(K), \forall K \in \mathcal{T}_h\}. $$
The best choice for the equilibrated flux reconstruction would then be

$$ \sigma_h = \arg \min_{v_h \in \mathcal{M}_h} \|v_h + \nabla u_0^h\|_{\Omega_0} \quad \text{subject to} \quad \left\{ \begin{array}{ll} \nabla \cdot v_h = f & \text{in } \Omega_0 \\ v_h \cdot n = -g & \text{on } \Gamma_N \setminus \gamma \\ v_h \cdot n_0 = -g_0 & \text{on } \gamma_0. \end{array} \right. \quad (33) $$

However, finding $\sigma_h$ through (33) implies solving a global optimization problem in the domain $\Omega_0$.

Following [3, 8] we adopt instead a different strategy, in which local flux reconstructions are built on patches $\omega_a$ of elements sharing a vertex $a \in \mathcal{N}_h$. Let us denote by $\psi_a$ the hat function in $P_1(\tilde{T}_h) \cap H^1(\Omega_0)$ taking value 1 in vertex $a$ and 0 on all the other vertices. Let us denote by $\partial\omega_a$ the boundary of the patch $\omega_a$ and let $\partial\omega_a^0 \subseteq \partial\omega_a$ be defined as

$$ \partial\omega_a^0 = \{x \in \partial\omega_a : \psi_a(x) = 0\}, $$

and $\partial\omega_a^\psi = \partial\omega_a \setminus \partial\omega_a^0$. Let us remark that, if $a \in \mathcal{N}_h^{\text{int}}$ then $\partial\omega_a^0 = \partial\omega_a$. Let $\Gamma_N^0 = (\Gamma_N \setminus \gamma) \cup \gamma_0$ and let us introduce

$$ \mathcal{M}_h^{a,0} = \{v_h \in \mathcal{M}_h(\omega_a) : v_h \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a^0 \cup (\partial\omega_a^\psi \cap \Gamma_N^0)\} $$

and

$$ \mathcal{M}_h^a := \left\{ \begin{array}{ll} M_h^{a,0} & \text{if } a \in \mathcal{N}_h^{\text{int}} \\ \{v_h \in \mathcal{M}_h(\omega_a) : \begin{array}{l} v_h \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a^0, \\ v_h \cdot n_{\omega_a} = -\psi_a g \text{ on } \partial\omega_a^\psi \cap (\Gamma_N \setminus \gamma), \\ v_h \cdot n_{\omega_a} = -\psi_a g_0 \text{ on } \partial\omega_a^\psi \cap \gamma_0 \end{array} \} & \text{if } a \in \mathcal{N}_h^{\text{ext}}, \end{array} \right. $$

$$ Q_h^a := \left\{ \begin{array}{ll} \{q_h \in Q_h(\omega_a) : (q_h, 1)_{\omega_a} = 0\} & \text{if } a \in \mathcal{N}_h^{\text{int}} \text{ or } a \in \text{int}(\overline{\Gamma_N^0}) \\ Q_h(\omega_a) & \text{if } a \in \mathcal{N}_h^{\text{ext}} \text{ and } a \notin \text{int}(\overline{\Gamma_N^0}), \end{array} \right. $$

where $\mathcal{M}_h(\omega_a)$ and $Q_h(\omega_a)$ are respectively the restrictions of $\mathcal{M}_h$ and $Q_h$ to the patch $\omega_a$ and $Q_h$ is defined as in (6). We then look for local equilibrated flux reconstructions as

$$ \sigma_h^a = \arg \min_{v_h \in \mathcal{M}_h^a} \|v_h + \psi_a \nabla u_0^h\|_{0,\omega_a} \quad \text{subject to} \quad \nabla \cdot v_h = \psi_a f - \nabla \psi_a \cdot \nabla u_0^h \quad (34) $$

and then we set $\sigma_h = \sum_{a \in \mathcal{N}_h} \sigma_h^a$.

The optimization problem (34) is equivalent to look for $\sigma_h^a \in \mathcal{M}_h^a$ and $\lambda_h^a \in Q_h^a$ such that,

$$ (\sigma_h^a, v_h)_{\omega_a} - (\lambda_h^a, \nabla \cdot v_h)_{\omega_a} = -(\psi_a \nabla u_0^h, v_h)_{\omega_a} \quad \forall v_h \in \mathcal{M}_h^{a,0} \quad (35) $$

$$ (\nabla \cdot \sigma_h^a, q_h)_{\omega_a} = (\psi_a f, q_h)_{\omega_a} - (\nabla \psi_a \cdot \nabla u_0^h, q_h)_{\omega_a} \quad \forall q_h \in Q_h^a, \quad (36) $$

which is the strategy that we are actually adopting in practice.

The equilibrated flux on the extension $\tilde{F}$ of a positive feature $F$ is reconstructed exactly in the same manner. Denoting by $\tilde{\mathcal{N}}_h^{\text{int}}$ and $\tilde{\mathcal{N}}_h^{\text{ext}}$ respectively the internal and the boundary vertices of the mesh $\tilde{T}_h$ of $\tilde{F}$, introducing

$$ \tilde{M}_h(\tilde{T}_h) := \{v_h \in \mathcal{H}(\text{div}, \tilde{F}) : v_h|_K \in [\mathcal{P}_p(K)]^d + x\mathcal{P}_p(K), \forall K \in \tilde{T}_h\}, $$
and recalling the definition of $\tilde{Q}_h$ given in (10), we look for the couple $(\tilde{\sigma}_h^a, \tilde{\lambda}_h^a)$ in the sets and spaces

$$ \tilde{M}_h^a := \left\{ \begin{array}{ll} \{ v_h \in \tilde{M}_h(\omega_a) : v_h \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a^0 \cup (\partial\omega_a^\psi \cap (\tilde{\gamma} \cup \gamma_s)) \} & \text{if } a \in \tilde{\mathcal{N}}_h^{\text{int}} \\ \\ \{ v_h \in \tilde{M}_h(\omega_a) : \begin{aligned}[t] v_h \cdot n_{\omega_a} &= 0 \text{ on } \partial\omega_a^0, \\ v_h \cdot n_{\omega_a} &= -\psi_a \tilde{g} \text{ on } \partial\omega_a^\psi \cap \tilde{\gamma}, \\ v_h \cdot n_{\omega_a} &= -\psi_a g \text{ on } \partial\omega_a^\psi \cap \gamma_s \end{aligned} \} & \text{if } a \in \tilde{\mathcal{N}}_h^{\text{ext}}, \end{array} \right. \\ \\ \tilde{Q}_h^a := \left\{ \begin{array}{ll} \{ q_h \in \tilde{Q}_h(\omega_a) : (q_h, 1)_{\omega_a} = 0 \} & \text{if } a \in \tilde{\mathcal{N}}_h^{\text{int}} \text{ or } a \in \text{int}(\tilde{\gamma} \cup \gamma_s) \\ \\ \tilde{Q}_h(\omega_a) & \text{if } a \in \tilde{\mathcal{N}}_h^{\text{ext}} \text{ and } a \notin \text{int}(\tilde{\gamma} \cup \gamma_s) \end{array} \right. $$

solving a problem analogous to (35)-(36).

# 6 Numerical experiments

In this section we propose some numerical experiments to validate the proposed estimator. We here focus on the case $d = 2$ and $p = 1$. All the simulations were performed in Matlab and meshes were built using the Triangle mesh generator [19]. For each element $K$ of a mesh $T_h$ we denote by $h_K$ the diameter of the element, and we choose as a mesh parameter $h = \max_{K \in \mathcal{T}_h} h_K$.

In the following we define the *effectivity index* as the ratio between the total estimator and the overall error. i.e.

$$ \eta = \frac{\mathcal{E}_{\text{tot}}}{\|\nabla(u - u_0^h)\|_{\Omega}}. \qquad (37) $$

In case of a domain $\Omega$ characterized by a single negative feature, following Proposition 1, the total estimator is defined as

$$ \mathcal{E}_{\text{tot}} = C_D \mathcal{E}_{\gamma} + \mathcal{E}_0, \qquad (38) $$

while in the case of a single positive feature it is

$$ \mathcal{E}_{\text{tot}} = \tilde{C}_D \tilde{\mathcal{E}}_{\gamma_0} + (\tilde{\mathcal{E}}_0^2 + \mathcal{E}_0^2)^{\frac{1}{2}}, \qquad (39) $$

assuming $\tilde{F} = F$ (see Remark 2). For all the proposed experiments, a reference solution is built by solving the problem on the original geometry $\Omega$ by linear finite elements on a very fine mesh. With an abuse of notation this reference solution is still denoted by $u$.

Three numerical examples are proposed. In Test 1 we consider the case of a single negative internal feature, analyzing the convergence of the total estimator and of the overall error under mesh refinement and feature size reduction. Test 2 deals instead with the case of positive and negative boundary features and the convergence of the estimator and of the error are again analyzed. Finally, in Test 3 we consider the presence of multiple internal negative features, showing how the proposed estimator allows to point out which features have the greatest impact on the error.

## 6.1 Test 1: negative internal feature

For this first numerical example we consider a square domain characterized by a single negative (see Figure 3a). We denote by $\epsilon$ the characteristic size of the feature, i.e. the radius of the circle circumscribing the feature itself.
<figure><img src="image_4.png" /><figcaption>Figure 3: Test 1: reference solution $u_\epsilon$ ($\epsilon = 0.14$), and defeatured geometry with an example of computational mesh used to compute the defeatured solution and the equilibrated flux reconstruction.</figcaption></figure>

Setting $\Omega_0 = (0, 1)^2$ and $\Omega_\epsilon = \Omega_0 \setminus \overline{F_\epsilon}$ we consider on the exact geometry $\Omega_\epsilon$ the problem

$$ \left\{ \begin{array}{ll} -\Delta u_{\epsilon} = f & \text{in } \Omega_{\epsilon} \\ u_{\epsilon} = 0 & \text{on } \partial\Omega_{\epsilon} \setminus \gamma_{\epsilon} \\ \nabla u_{\epsilon} \cdot \mathbf{n} = 0 & \text{on } \gamma_{\epsilon}, \end{array} \right. \qquad (40) $$

with $f(x, y) = x$. For what concerns the defeatured problem, $\gamma_0 = \emptyset$, since the feature is internal.

Figure 3b reports an example of the computational mesh $T_h$ used to solve the defeatured problem on $\Omega_0$. We remark how the mesh has no need to be conforming to the feature boundary since the equilibrated flux $\sigma_h$ is reconstructed on the defeatured geometry, which is blind to the feature, and $\mathcal{E}_\gamma$ is computed by simply defining a proper quadrature rule on the feature boundary itself and evaluating the normal trace of $\sigma_h$ in the chosen quadrature nodes.

In the numerical experiments that follow we choose $C_D = 1$ in (38). Figure 4 shows the convergence of the estimator $\mathcal{E}_{\text{tot}}$ and of the energy norm of the overall error $||\nabla(u_\epsilon - u_0^h)||_\Omega$, under mesh refinement. The values of $\mathcal{E}_\gamma$ and $\mathcal{E}_0$ are also reported. Three fixed values of $\epsilon$ are considered, namely $\epsilon = 7.00 \cdot 10^{-2}, 1.75 \cdot 10^{-2}, 4.83 \cdot 10^{-3}$. As expected, going from a coarse to a fine mesh, the error reaches a plateau when the defeaturing error becomes more relevant than the numerical one. The bigger the feature is, the earlier the plateau is reached. The same behavior is captured also by the estimator. For a fixed feature size, the value of $\mathcal{E}_\gamma$ remains constant while $\mathcal{E}_0$ converges as $O(h)$, as expected.

The trend of the effectivity index $\eta$ under mesh refinement and for the three considered feature sizes is reported in Figure 6a. As expected, since the numerical source of the error is sharply bounded by $\mathcal{E}_0$, when $\mathcal{E}_0 \gg \mathcal{E}_\gamma$ we have $\eta \sim 1$. The effectivity index appears instead to be around 2.5 when the defeaturing component is dominating. The highest values of $\eta$, namely $\eta \sim 3$, are registered when $\mathcal{E}_\gamma > \mathcal{E}_0$ but the two components have still a comparable magnitude.

Figure 5 focuses instead on the convergence of the estimator and of the error under the reduction of the feature size, for three fixed mesh sizes, namely $h = 1.25 \cdot 10^{-1}, 3.13 \cdot 10^{-2}, 7.81 \cdot 10^{-3}$. As expected, both the error and the estimator reach a plateau when the numerical error dominates over the defeaturing one. The value of the effectivity index $\eta$ is reported in Figure 6b, with the same considerations done for the Figure 6a still holding.
<figure><img src="image_5.png" /><figcaption>Figure 4: Test 1: energy norm of the error (full red line), total estimator (full blue line) and its components (dashed lines) under mesh refinement and for three fixed feature sizes</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Test 1: energy norm of the error, total estimator and its components under the reduction of the feature size and for three different fixed mesh sizes.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Test 1: Effectivity index under mesh refinement and under feature size reduction.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Test 2: Exact geometry with a negative and a positive feature.</figcaption></figure>

## 6.2 Test 2: positive and negative boundary features

In this second numerical example we consider the case of a positive and a negative boundary feature. As in the previous test case we choose $\Omega_0 = (0, 1)^2$, while we define

$$F_n = \left( \frac{1-\epsilon}{2}, \frac{1+\epsilon}{2} \right) \times (1-\epsilon, 1),$$

$$F_p = \left( \frac{1-\epsilon}{2}, \frac{1+\epsilon}{2} \right) \times (-\epsilon, 0).$$

For the negative feature case we chose as exact geometry $\Omega_n = \Omega_0 \setminus \overline{F_n}$, while for the positive feature case we choose $\Omega_n = \text{int}(\overline{\Omega_0} \cup \overline{F_p})$, as reported in Figure 7. In both cases we consider the problem

$$\left\{ \begin{array}{ll} -\Delta u = f & \text{in } \Omega_{\star} \\ u = 0 & \text{on } \Gamma_D \\ \nabla u \cdot n = 0 & \text{on } \Gamma_N \end{array} \right. \qquad (41)$$

with $\star = \{n, p\}$, $f = 1$, $\Gamma_D = \{(x, y) : x = 0 \lor x = 1\}$ and $\Gamma_N = \partial\Omega_\epsilon \setminus \Gamma_D$. We recall that, according to its definition, $\Gamma_N$ includes also the feature boundary.

Figure 8 compares the convergence under mesh refinement of the total estimator and error in the negative and in the positive feature case and for two different values of $\epsilon$. Let us recall how, in presence of a negative feature the total estimator is defined as in (38), with $\mathcal{E}_\gamma = \mathcal{E}_{\gamma_n}$ computed on $\gamma_n$ (see Figure 7a) from the equilibrated flux reconstructed on $\Omega_0$; in the case of a positive feature, instead, the definition is provided by (39) with $\tilde{\mathcal{E}}_{\gamma_0} = \tilde{\mathcal{E}}_{\gamma_{0,p}}$ computed on $\gamma_{0,p}$ (see Figure 7b) from the equilibrated flux reconstructed on the feature itself. In the following we choose $C_D = \tilde{C}_D = 1$ and, for the sake of clarity, we denote respectively by $\mathcal{E}_{\text{tot}}^n$ and $\mathcal{E}_{\text{tot}}^p$ the total estimator in the negative and in the positive feature case. As in the previous test case, fixing the feature size and refining the mesh, we can observe how the overall error reaches a plateau, and how this behavior is captured also by the total estimator, both in the negative and in the positive feature case. As expected, a bigger feature (Figure 8a) produces a stagnation of the error and of the estimator already for coarse
<figure><img src="image_9.png" /><figcaption>Figure 8: Test 2: energy norm of the error and total estimator under mesh refinement, for the negative and the positive feature cases and for two different feature sizes.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Test 2: effectivity index for the negative and the positive feature cases for two different feature sizes.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Test 2, simultaneous presence of positive and negative feature: energy norm of the error, total estimator with its components and effectivity index. For both features, $\epsilon = 0.2$.</figcaption></figure>
meshes, while if the feature is smaller (Figure 8b) the defeaturing source of error becomes relevant only for finer meshes.

Figure 9 shows the trend of the effectivity index related to the curves reported in Figure 8. Both for the negative and the positive feature case we observe that, as expected, $\eta \sim 1$ when the numerical component is dominating (coarse meshes in Figure 8b). We can instead observe how $\eta \sim 1.5$ when the defeaturing component dominates (fine meshes in Figure 8a), and how the effectivity index is in general lower with respect to the internal negative feature case (Test 1), with $\eta < 2$ even when both the defeaturing and the numerical component have a significant impact.

Finally, Figure 10 refers to the case in which the negative and the positive features are simultaneously present, i.e. $\Omega = \text{int}(\Omega_0 \cup \overline{F_p}) \setminus \overline{F_n}$. For both features we choose $\epsilon = 0.2$. Figure 10a reports the convergence of the error and of the estimator under mesh refinement. The total estimator is, in this case, defined as

$$ \mathcal{E}_{\text{tot}} = C(\tilde{\mathcal{E}}_{\gamma_{0,p}}^2 + \mathcal{E}_{\gamma_n}^2)^{\frac{1}{2}} + (\mathcal{E}_0^2 + \tilde{\mathcal{E}}_0^2)^{\frac{1}{2}}, $$

with $C > 0$ being a constant independent of the size of both features. In particular, we choose here $C = 1$. As in the previous test cases, the estimator appears to correctly capture the behavior of the overall error. The corresponding effectivity index is reported in Figure 10b.

## 6.3 Test 3: multiple internal features

For this last numerical example we consider a case with multiple internal features, similar to the one proposed in [2, 4]. Our aim is to show the capability of the proposed estimator to identify the most relevant features and to provide a criterion to decide whether a feature should be added or not, according also to the magnitude of the numerical source of the error.

Let us define the defeatured geometry again as $\Omega_0 = (0, 1)^2$ and let us consider a set of $I$ features $\mathcal{F} = \{F_i\}_{i \in \mathcal{I}}$, $\mathcal{I} = \{1, ..., I\}$, each of which is a polygon of 16 faces, inscribed in a circle of radius $\epsilon_i$ and centered in $\mathbf{x}_C^i$. In particular we choose $I = 5$ and

$$ \begin{aligned} \mathbf{x}_C^1 &= (0.12, 0.12), \quad \epsilon_1 = 0.02; & \mathbf{x}_C^2 &= (0.35, 0.35), \quad \epsilon_2 = 0.05 \\ \mathbf{x}_C^3 &= (0.65, 0.65), \quad \epsilon_3 = 0.10; & \mathbf{x}_C^4 &= (0.20, 0.68), \quad \epsilon_4 = 0.05; \\ \mathbf{x}_C^5 &= (0.65, 0.16), \quad \epsilon_5 = 0.05. \end{aligned} $$

The boundary of the $i$-th feature is denoted by $\gamma_i$. On the exact geometry $\Omega = \Omega_0 \setminus \overline{\cup_{i \in \mathcal{I}} F_i}$, which is reported in Figure 11a, we consider the problem

$$ \left\{ \begin{array}{ll} \Delta u = 0 & \text{in } \Omega \\ u = g_D & \text{on } \Gamma_D \\ \nabla u \cdot \mathbf{n} = 0 & \text{on } \Gamma_N \end{array} \right. \qquad (42) $$

with $\Gamma_D = \{(x, y) : x = 0 \lor y = 0\}$, $g_D(x, y) = e^{-8(x+y)}$ and $\Gamma_N = \partial\Omega \setminus \Gamma_D$, including also the feature boundaries. The reference solution $u$ is also reported in Figure 11a.

Let $\mathcal{M} \subseteq \mathcal{F}$ be a subset of features indexed by $j \in \mathcal{I}^* \subseteq \mathcal{I}$ and let us denote by $\Omega_0^\mathcal{M}$ a generic partially defeatured geometry, obtained by including the features in $\mathcal{M}$ to the defeatured geometry, i.e. $\Omega_0^\mathcal{M} = \Omega_0 \setminus \overline{\cup_{j \in \mathcal{I}^*} F_j}$. If $\mathcal{M} = \emptyset$, then $\Omega_0^\mathcal{M} = \Omega_0$. An example of a computational mesh defined on $\Omega_0$ is reported in Figure 11b. We will use $u_0^h$ to refer both to the numeric solution computed
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Test 3: reference solution $u$ and totally defeatured geometry $\Omega_0$ with an example of computational mesh.</figcaption></figure>

on $\Omega_0$ and to the one computed on $\Omega_0^M$, the meaning being clear from the context. In presence of multiple negative features the total estimator is defined as $\mathcal{E}_{\text{tot}} = \alpha_D \mathcal{E}_\gamma + \mathcal{E}_0$, where

$$ \mathcal{E}_{\gamma} = \left( \sum_{k \in \mathcal{I} \setminus \mathcal{I}_{*}} \mathcal{E}_{\gamma_k}^2 \right)^{\frac{1}{2}}. $$

Table 1 reports the value of the components of the total estimator for differently refined meshes and for different choices of $\mathcal{M}$, i.e. of the partially defeatured geometry on which $u_0^h$ is computed. The rows of the table are divided into three sets, corresponding to three differently refined meshes. The variation in the number of degrees of freedom which can be observed when a feature is included into the geometry is related to the adaptation of the mesh to the feature boundary and to the deletion of the degrees of freedom lying inside the feature itself. Looking at the columns from 5 to 9 we can see how feature $F_1$ is clearly the most relevant, since $\mathcal{E}_{\gamma_1} > \mathcal{E}_{\gamma_i}$ for all $i > 1$. This is expected since, although being the smallest feature, it is located in a region in which the gradient of the solution is very steep. Feature $F_3$ is, instead, almost irrelevant: despite being the biggest one it is located in a region in which the solution is rather flat and hence its impact on the solution accuracy tends to be negligible. As expected, the values of $\mathcal{E}_{\gamma_i}$ are independent from the mesh size, meaning that the relevance of the features can be evaluated even on a coarse mesh. However, the choice of including the $i$-th feature into the geometry should be taken by comparing $\mathcal{E}_{\gamma_i}$ with $\mathcal{E}_0$, which is a sharp indicator of the numerical source of error. In particular, a value of $\mathcal{E}_{\gamma_i}$ considerably bigger than $\mathcal{E}_0$ means that we will not be able to significantly reduce the error by mesh refinement, unless the feature is added. This is true, for example, for feature $F_1$ with the second considered mesh and for features $F_1$ and $F_2$ for the finest mesh.

Table 2 focuses exactly on these cases, reporting the values of the energy norm of the overall error and of the total estimator, along with the corresponding effectivity index. In particular, $\mathcal{E}_{\text{tot}}$ is computed with $\alpha_D = 1$. In Table 1, we can observe that, if $\mathcal{M} = \emptyset$, the reduction of $\mathcal{E}_0$ when going from $h = 3.13 \cdot 10^{-2}$ ($\sim 5k$ degrees of freedom) to $h = 1.56 \cdot 10^{-2}$ ($\sim 20k$ degrees of freedom) is of about 50%. This is expected, since $\mathcal{E}_0$ should converge as $\mathcal{O}(h)$. However, looking at rows 1 and 3 in Table 2, we see that the drop in the total estimator, and hence in the error, is under 20%.
<table><caption>Table 1: Test 3: components of the total estimator for differently refined mesh and for different choices of the (partially) defatured geometry.</caption><thead><tr><th>h</th><th>M</th><th>N<sub>dof</sub></th><th>E<sub>0</sub></th><th>E<sub>γ<sub>1</sub></sub></th><th>E<sub>γ<sub>2</sub></sub></th><th>E<sub>γ<sub>3</sub></sub></th><th>E<sub>γ<sub>4</sub></sub></th><th>E<sub>γ<sub>5</sub></sub></th><th>E<sub>γ</sub></th></tr></thead><tbody><tr><td rowspan="2">6.25e-2</td><td>Ø</td><td>1240</td><td>0.101</td><td>0.147</td><td>0.050</td><td>0.008</td><td>0.026</td><td>0.036</td><td>0.162</td></tr><tr><td>F<sub>1</sub></td><td>1380</td><td>0.096</td><td>-</td><td>0.048</td><td>0.008</td><td>0.025</td><td>0.035</td><td>0.065</td></tr><tr><td rowspan="3">3.13e-2</td><td>Ø</td><td>4960</td><td>0.052</td><td>0.146</td><td>0.050</td><td>0.008</td><td>0.025</td><td>0.036</td><td>0.162</td></tr><tr><td>F<sub>1</sub></td><td>5501</td><td>0.050</td><td>-</td><td>0.048</td><td>0.008</td><td>0.025</td><td>0.035</td><td>0.065</td></tr><tr><td>F<sub>1</sub>, F<sub>2</sub></td><td>5684</td><td>0.053</td><td>-</td><td>-</td><td>0.007</td><td>0.024</td><td>0.035</td><td>0.043</td></tr><tr><td rowspan="5">1.56e-2</td><td>Ø</td><td>19840</td><td>0.027</td><td>0.146</td><td>0.050</td><td>0.008</td><td>0.025</td><td>0.036</td><td>0.161</td></tr><tr><td>F<sub>1</sub></td><td>21963</td><td>0.026</td><td>0</td><td>0.048</td><td>0.008</td><td>0.025</td><td>0.035</td><td>0.065</td></tr><tr><td>F<sub>1</sub>, F<sub>2</sub></td><td>22618</td><td>0.027</td><td>-</td><td>-</td><td>0.007</td><td>0.024</td><td>0.034</td><td>0.042</td></tr><tr><td>F<sub>1</sub>, F<sub>2</sub>, F<sub>5</sub></td><td>23953</td><td>0.027</td><td>-</td><td>-</td><td>0.007</td><td>0.024</td><td>-</td><td>0.025</td></tr><tr><td>F<sub>1</sub>, F<sub>2</sub>, F<sub>4</sub>, F<sub>5</sub></td><td>24372</td><td>0.027</td><td>-</td><td>-</td><td>0.007</td><td>-</td><td>-</td><td>0.007</td></tr></tbody></table>

Adding feature $F_1$ and refining the mesh at the same time the drop is instead of about 60%, as it can be seen by comparing rows 1 and 4 in Table 2. Let us remark that, for the finest considered mesh, also feature $F_2$ becomes rather relevant. However, $E_{\gamma_2}$ is closer to $E_0$ and hence adding it to the geometry has a smaller impact on the solution accuracy.

This experiment is to be intended as a preliminary test for the use of the proposed estimator in an adaptive strategy, involving both geometrical adaptation (i.e. feature inclusion) and local mesh refinement. We decide to leave this to a forthcoming work: indeed the procedure adopted for the computation of the equilibrated flux requires the mesh to be conforming to the domain boundaries, and hence also to the boundary of the features which are actually included in the partially defatured geometry. However, to build an efficient and flexible adaptive algorithm we do not want to remesh the geometry each time a feature is added, and for this reason a generalization of the equilibrated flux reconstruction to trimmed meshes needs to be considered.

<table><caption>Table 2: Test 3: energy norm of the error, total estimator and effectivity index for differently refined mesh and for different choices of the (partially) defatured geometry.</caption><thead><tr><th>h</th><th>M</th><th>||∇(u - u<sub>0</sub><sup>h</sup>)||<sub>Ω</sub></th><th>E<sub>tot</sub></th><th>η</th></tr></thead><tbody><tr><td rowspan="2">3.13e-2</td><td>Ø</td><td>0.079</td><td>0.214</td><td>2.71</td></tr><tr><td>F<sub>1</sub></td><td>0.052</td><td>0.115</td><td>2.21</td></tr><tr><td rowspan="3">1.56e-2</td><td>Ø</td><td>0.066</td><td>0.188</td><td>2.85</td></tr><tr><td>F<sub>1</sub></td><td>0.033</td><td>0.091</td><td>2.76</td></tr><tr><td>F<sub>1</sub>, F<sub>2</sub></td><td>0.028</td><td>0.069</td><td>2.46</td></tr></tbody></table>
Conclusions

In this work we have proposed a new *a posteriori* error estimator for defeaturing problems based
on an equilibrated flux reconstruction and designed for finite elements. The Poisson equation with
Neumann boundary conditions on the feature boundary was taken as a model problem. The relia-
bility of the estimator has been proven both in the negative and in the positive feature case, and
tested with several numerical examples. The choice of using an equilibrated flux reconstruction leads
to an estimator which is able to bound sharply the numerical component of the error and which
never requires to evaluate the normal trace of the numerical flux, which is typically discontinuous
on element edges in a standard finite element discretization.

This work is to be intended as a preliminary analysis for the use of the proposed estimator in
an adaptive strategy, allowing not only for mesh refinement, but also for an automatic inclusion
of those features whose absence causes most of the accuracy loss. The proposed estimator does
not require the mesh to be conforming to the feature boundary until the feature is included in the
computational domain itself. Indeed, computing the integral of the normal trace of the equilibrated
flux reconstruction on a generic curve is always possible, regardless of the intersections with the mesh
elements. However, the procedure which was adopted to reconstruct the equilibrated flux is designed
for meshes which are conforming to the computational domain boundary and this would require to
remesh the domain each time a feature is added by the adaptive procedure, hence increasing the
complexity of the algorithm. For this reason, an extension of the equilibrated flux reconstruction
to the case of trimmed meshes needs to be considered, so that the geometry never needs to be
remeshed. This generalization is left to a forthcoming work, that is currently under preparation.
Although the proof of the reliability of the estimator holds in $\mathbb{R}^d$, $d = 2, 3$, and for any polynomial
order $p \ge 1$, we decided to propose numerical experiments only in $\mathbb{R}^2$ and for $p = 1$. The application
of the estimator on more complex, realistic and tridimensional geometries, or the use of a higher
order finite element approximation, are left to a forthcoming work as well, both extensions having
an impact only on implementation aspects.

Acknowledgments

The authors are grateful to Zhaonan Dong (Inria, Paris) for sharing his code on equilibrated fluxes
and to Paolo Bardella (DET, Politecnico di Torino) for the MeshToolbox library.

The authors acknowledge the support of the Swiss National Science Foundation (via project MINT
n. 200021_215099, PDE tools for analysis-aware geometry processing in simulation science) and of
European Union Horizon 2020 FET program (under grant agreement No 862025 (ADAM2)).

Author O. Chanon acknowledges the support of the Swiss National Science Foundation through
the project n.P500PT 210974.

References

[1] M. Ainsworth and J.T. Oden. "A posteriori error estimation in finite element analysis". In: *Computer Methods in Applied Mechanics and Engineering* 142.1 (1997), pp. 1-88.

[2] P. Antolín and O. Chanon. “Analysis-aware defeaturing of complex geometries with Neumann features”. In: *Int. J. Numer. Methods Eng.* (2023).
[3] D. Braess and J. Schöberl. "Equilibrated Residual Error Estimator for Edge Elements". In: *Mathematics of Computation* 77.262 (2008), pp. 651-672.

[4] A. Buffa, O. Chanon, and R. Vázquez. *Adaptive analysis-aware defeaturing*. 2022. arXiv: 2212.05183.

[5] A. Buffa, O. Chanon, and R. Vázquez. "Analysis-aware defeaturing: Problem setting and a posteriori estimation." In: *Math. Models Methods Appl. Sci.* 32.2 (2022), pp. 359-402.

[6] K. K. Choi and N.-H. Kim. "Structural Sensitivity Analysis and Optimization 1: Linear Systems". In: *Springer Science and Business Media* (2005).

[7] P. Destuynder and B. Métivet. "Explicit Error Bounds in a Conforming Finite Element Method". In: *Mathematics of Computation* 68.228 (1999), pp. 1379-1396.

[8] A. Ern and M. Vohralík. "Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations". In: *SIAM Journal on Numerical Analysis* 53.2 (2015), pp. 1058-1081.

[9] R. Ferrandes et al. "A posteriori evaluation of simplification details for finite element model preparation". In: *Computers and Structures* 87.1 (2009), pp. 73-80.

[10] L. Fine, L. Remondini, and J.-C. Leon. "Automated generation of FEA models through idealization operators". In: *International Journal for Numerical Methods in Engineering* 49.1-2 (2000), pp. 83-108.

[11] G. Foucault, P. M. Marin, and J.-C Léon. "Mechanical Criteria for the Preparation of Finite Element Models". In: *International Meshing Roundtable Conference*. 2004.

[12] S. H. Gopalakrishnan and K. Suresh. "A formal theory for estimating defeaturing-induced engineering analysis errors". In: *Computer-Aided Design* 39.1 (2007), 60-68.

[13] S. H. Gopalakrishnan and K. Suresh. "Feature sensitivity: a generalization of topological sensitivity". In: *Finite Elements in Analysis and Design* 44.11 (2008), 696-704.

[14] M. Li and S. Gao. "Estimating defeaturing-induced engineering analysis errors for arbitrary 3D features". In: *Computer-Aided Design* 43.12 (2013), 1587-1597.

[15] M. Li, S. Gao, and R. R. Martin. "Estimating the effects of removing negative features on engineering analysis". In: *Computer-Aided Design* 43.11 (2011). Solid and Physical Modeling, 1402-1412.

[16] M. Li, S. Gao, and K. Zhang. "A goal-oriented error estimator for the analysis of simplified designs". In: *Computer Methods in Applied Mechanics and Engineering* 255 (2013), 89-103.

[17] S. Gao Li and R. R. Martin. "Engineering analysis error estimation when removing finite-sized features in nonlinear elliptic problem". In: *Computer-Aided Design* 45.2 (2013), 361-372.

[18] R. Luce and B. I. Wohlmuth. "A Local a Posteriori Error Estimator Based on Equilibrated Fluxes". In: *SIAM Journal on Numerical Analysis* 42.4 (2005), pp. 1394-1414.

[19] J. R. Shewchuk. "Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator". In: *Applied Computational Geometry: Towards Geometric Engineering*. Ed. by Ming C. Lin and Dinesh Manocha. Vol. 1148. Lecture Notes in Computer Science. From the First ACM Workshop on Applied Computational Geometry. Springer-Verlag, May 1996, pp. 203-222.

[20] J. Sokolowski and A. Zochowski. "On the topological derivative in shape optimization". In: *SIAM journal on control and optimization* 37.4 (1999), 1251-1272.
[21] J. Tang, S. Gao, and M. Li. "Evaluating defeaturing-induced impact on model analysis". In: *Mathematical and Computer Modelling* 57.3 (2013), pp. 413-424.

[22] A. Thakur, A. G. Banerjee, and S. K. Gupta. "A survey of CAD model simplification techniques for physics-based simulation applications". In: *Computer-Aided Design* 41.2 (2009), pp. 65-80.

[23] I. Turevsky, S. H. Gopalakrishnan, and K. Suresh. "An efficient numerical method for computing the topological sensitivity of arbitrary-shaped features in plate bending". In: *International journal for numerical methods in engineering* 79.13 (2011), 1683-1702.