Analysis-aware defeaturing of complex geometries with Neumann features

Local modifications of a computational domain are often performed in order to simplify the meshing process and to reduce computational costs and memory requirements. However, removing geometrical features of a domain often introduces a non-negligible error in the solution of a differential problem in which it is defined. In this work, we extend the results from [1] by studying the case of domains containing an arbitrary number of distinct Neumann features, and by performing an analysis on Poisson's, linear elasticity, and Stokes' equations. We introduce a simple, computationally cheap, reliable, and efficient a posteriori estimator of the geometrical defeaturing error. Moreover, we also introduce a geometric refinement strategy that accounts for the defeaturing error: Starting from a fully defeatured geometry, the algorithm determines at each iteration step which features need to be added to the geometrical model to reduce the defeaturing error. These important features are then added to the (partially) defeatured geometrical model at the next iteration, until the solution attains a prescribed accuracy. A wide range of two- and three-dimensional numerical experiments are finally reported to illustrate this work.


Introduction
With the advance of engineering knowledge, simulations are performed on objects of increasing geometric complexity, nowadays mainly described by three-dimensional Computer-Aided Design (CAD) models.These models often contain a large number of geometric details of different scales, also called geometric features.Unfortunately, the construction of a finite element mesh on such complex domains may fail, or if it does not, the mesh generation may be very difficult; see for example [2] dealing with the complexity arising from an automatic all-hexahedral mesh generation for complex B-Reps (boundary representations).Moreover, the resulting mesh may require a very large number of elements, therefore leading to simulations which are too costly or even unfeasible.For instance, it has been shown in [3,4] that the cost of the underlying simulation may be increased by up to a factor 10 in the presence of a single geometric feature of relatively small size.
However, depending on the problem at hand, such high model complexity may be unnecessary.That is, the geometric description of the object may require a high number of degrees of freedom, but not all of them are needed to perform an accurate analysis, and taking all of them into account is potentially too costly.To deal with complex geometries and to accelerate the process of analysis-aware geometric design, it is therefore essential to be able to simplify the geometric model, process also known as defeaturing.This is a very common practice among finite element analysts.See, as matter of example, the case illustrated in Figure 1.There, a CAD design with numerous features as holes, rounds, and a carved logo (Figure 1a) is defeatured to create a simpler model (Figure 1b) that is easier to mesh.Each finite element mesh in Figure 1 was generated using Gmsh [5] with the same mesh algorithm and parameters.Nevertheless, the mesh of the original design has 5 times more nodes than the one of the defeatured model.Likely, most of those extra Analysis-aware defeaturing of complex geometries with Neumann features  Both finite element meshes were generated with Gmsh [5] using the same mesh algorithm and parameters.degrees of freedom, required for correctly representing the geometrical details, will improve little the accuracy of the solution obtained with such mesh, but will increase the computational cost and memory requirements significantly.
Nevertheless, it is important to consider how such geometrical simplifications will impact the analysis solution, i.e., to control the error introduced by defeaturing, in order to provide an accurate solution of the problem at hand.The literature on the subject is still relatively scarce, and a lot remains to be done.To estimate the defeaturing error on the solution of a partial differential equation (PDE), some a posteriori criteria have been developed: The one introduced in [6] uses an approximation of the error in energy norm; in [7,8], an estimator is found using topological sensitivity analysis; adjoint theory is used in the series of works [9][10][11] to describe the first order defeaturing error on a quantity of interest; an estimator is introduced in [12] based on the reciprocal theorem stating the conservation of solution flux in the features; and in the series of works [13][14][15][16], defeaturing error is expressed as a modeling error directly on the differential problem, both for negative features (holes) and positive ones (protrusions).In these latter papers, the modeling error is then estimated with the dual weighted residual method [17].Nevertheless, very few of those works come with a sound mathematical theory, and most of them rely on some heuristics.
In the recent paper [1], the authors have tackled this issue: A precise mathematical framework is defined for geometries for which a single feature of very generic shape is removed, and an efficient and reliable a posteriori estimator of the defeaturing error is derived in the context of Poisson's equation in the energy norm.Worth mentioning related research papers include [18][19][20] that study heterogeneous and perforated materials, [21] that are interested in the error introduced by the approximation of boundary conditions, and [22] that more generally studies modeling errors coming from dimension reduction, homogenization and model simplification.
We generalize here the work [1].The original contributions of the presented research are the following: • We consider geometries for which an arbitrary number N of distinct Neumann geometrical features are removed from the computational domain, and we introduce a defeaturing error estimator whose effectivity index is independent from N .Note in particular that the features do not need to be small.
• We do not only consider Poisson's equation, but also linear elasticity and Stokes' problems.While the considered problems are exclusively linear, the proposed estimator could be extended to encompass nonlinear problems as well.
• Based on the proposed defeaturing error estimator, an adaptive strategy for geometric refinement is Analysis-aware defeaturing of complex geometries with Neumann features introduced: The algorithm begins with a completely simplified geometry and, at each iteration, identifies the necessary features to be added to the geometric model in order to reduce the defeaturing error.
• Finally, more complex numerical experiments with engineering interest are presented, in particular the simulation of a three-dimensional CAD design.
The error estimator formulation depends on the problem at hand.Hence, in this work, different expressions will be provided for the Poisson, Stokes, and linear elasticity problems.In particular, the differential problems and their corresponding defeaturing error estimator depend on the material properties and on the boundary conditions.E.g., some features will play more or less predominant roles depending on the chosen boundary conditions.
In general, and in practice, differential problems cannot be solved exactly, and thus a numerical method is used to solve them approximately in a finite dimensional space.However, in this article, we neglect the contribution of the error coming from the numerical approximation of the solution, as we concentrate our efforts on the defeaturing error.Deriving an estimator that also includes the numerical error contribution is the subject of a subsequent work, see [23,24].An early method for defeaturing and coarsening, called composite finite elements and developed by Hackbusch and Sauter in [25,26], deserves mentioning.It assumes that defeaturing arises solely from the limitation of the mesh in resolving geometric details within the computational domain.In contrast, our approach aims to treat defeaturing independently of any discretization, allowing us to separate the strict defeaturing error from the numerical approximation error.
Therefore, in Section 2, we first introduce the problem of defeaturing and the corresponding notation that will be used throughout the article, and we precisely define the defeaturing error that we aim at estimating.Then, we state the main results that are obtained, namely the reliability and the efficiency of an a posteriori estimator of the defeaturing error whose effectivity index is independent from the number of the features and their size.Subsequently, in Sections 3 and 4, we precisely define the exact and defeatured problems when the differential problem at hand is, respectively, Poisson's and the linear elasticity or Stokes' equations, and we propose in each case a defeaturing error estimator.We then introduce in Section 5 an adaptive geometric refinement strategy driven by the defeaturing error estimators previously defined.We perform in Section 6 a validation of the presented theory, of the aforementioned estimators' properties, and of the proposed adaptive strategy, thanks to an extensive set of two-and three-dimensional numerical experiments.To perform these tests, we use isogeometric analysis (IGA) [27,28] on very fine meshes, in order to have a negligible numerical error with respect to the defeaturing error.Since IGA and defeaturing both pursue the scope of reducing the gap between the design and the analysis phases, IGA is a natural method of choice.Nevertheless, it is important to remark that the proposed techniques are completely discretization agnostic and any other numerical method could be used.Some conclusions are finally drawn in Section 7. In the Appendix A, the reliability and efficiency of the proposed defeaturing error estimator are demonstrated in the framework of Stokes' equations, as the proofs for Poisson's and the linear elasticity equations are very similar, easier, and can be found in details in the corresponding thesis [24].
In the following, the operator ≲ is used to mean any inequality which neither depends on the number of features nor on their size, but which can depend on their shape, on their type (i.e., whether they are holes, also called negative features, or protrusions, also called positive features), and on the space dimension n = 2 or n = 3.Moreover, for all D ⊂ R n , and for all Λ ⊂ ∂D, we denote by |D| the n-dimensional Lebesgue measure of D, by |Λ| the (n − 1)-dimensional Hausdorff measure of Λ, by D and Γ the closure of D and of Γ, respectively, and by int(D) and int(Λ) the interior of D and of Γ, respectively.

Presentation of the problem
In this section, we present the considered problem, following and extending the setting in [1].Let us consider a potentially complicated, connected open Lipschitz domain Ω ⊂ R n , n = 2 or n = 3, on which we want to solve a differential problem P(Ω) which contains some boundary conditions.More precisely, let us assume that Ω  contains geometrical details of smaller scale also called features, and assume that the feature information of Ω is known a priori .As illustrated in Figure 2, we say that a feature or complex if it has both negative and positive components.A negative feature corresponds to a part where some material has been removed (e.g., a hole), a positive feature corresponds to the addition of some material (e.g., a protrusion), and a feature is complex in the most general situation in which there has been both the addition and the removal of some material.Note that an internal feature (e.g., an internal hole) is a special case of negative feature.We assume that the positive and the negative part of each feature is a connected open Lipschitz domain of R n .
However, solving the given differential problem P(Ω) can be very complicated due to the complexity of Ω, coming from the presence of the features.Therefore, we solve instead a similar problem but in a defeatured domain Ω 0 , where the features of Ω are removed: Holes are filled with some material and protrusions are cut out of the computational domain.This differential problem on Ω 0 is denoted by P(Ω 0 ) and it is called defeatured (or simplified) problem.
The defeaturing of the computational domain introduces an error in the problem's solution, and we are interested in controlling the energy norm of this so-called defeaturing error.In other words, if u is the solution of the exact problem P(Ω) and u 0 is the solution of the defeatured problem P(Ω 0 ), then we are interested in controlling the error "u − u 0 " in the exact energy norm, which is the energy norm defined by problem P(Ω) and denoted as |||•||| Ω .Since u is defined in Ω and u 0 is defined in Ω 0 , the defeaturing error needs to be more accurately defined.To do so, we need to introduce some further geometric notation.
Generalizing the work from [1], let N f ≥ 1, N f ∈ N, denote the total number of (possibly complex) features of Ω, gathered into the set F := F k N f k=1 .Let us make the following assumption on the features.Assumption 2.1 The features in F are separated, that is, ) one cannot have an increasingly large number of features that are arbitrarily close to one another.Remark 2.2 In the currently considered setting, features are discrete objects.Note that: • It is always possible to satisfy condition (a) of Assumption 2.1 by changing the numbering of the features.Indeed, if there are can be considered as a single feature that replaces the two features F k and F ℓ .However, note that the treatment of a geometry in which the boundary is complicated everywhere (for instance a fractal-like domain) is not considered here, see instead [29][30][31] for first results in this different framework.
• The precise mathematical definition of condition (b) will be given at the end of Section 2.3, as it requires technical geometric definitions that would hinder here the readability of the main results of this manuscript.
Since the features in F are assumed to be generally complex, this means that, for all k = 1, . . ., N f , the feature where F k n and F k p are open domains such that, if we define then In this setting, we define the defeatured geometry by Note that in the case in which the exact domain Ω contains a single feature F , i.e.N f = 1, then 1) and ( 2) and apply, therefore generalizing the definition of Ω 0 given in [1].
Remark 2.3 Given a complicated geometry Ω without any further information, one cannot always easily tell whether the features it contains are negative or positive, see Figure 4. Therefore, this is often a choice that the user needs to make, based on the available geometric information at hand.If one has access to a simplified geometry, for instance thanks to the history of CAD operations from which the exact geometry Ω is built, then it is possible to define the features from Ω and Ω 0 , instead of defining Ω 0 from Ω and the features.The identification of features in a given geometry and the construction of a corresponding simplified geometric model can be complicated tasks, see [32] for a review of possible techniques.However, this goes beyond the scope of this work, which supposes at its roots that the feature information is known a priori.Now, we remark that the solution u 0 of P(Ω 0 ) is not defined in F p , i.e., in the positive component of the features of Ω, since F p ̸ ⊂ Ω 0 but F p ⊂ Ω.Thus, to be able to compare u and u 0 in F p , we need to solve for each feature F k ∈ F a local differential problem P F k p in F k p , whose solution suitably extends u 0 to F k p .However, and as already noted in [1], the domain F k p might be complicated or even non-smooth, thus finding the solution of the extension problem P F k p might be cumbersome.
where we recall the definition of Ω ⋆ from (1).We are now finally able to precisely define the defeaturing error in the energy norm by |||u − u d ||| Ω .

Main results
In this article, we will precisely define the considered differential problems P(Ω), P(Ω 0 ), and P F k p for all k = 1, . . ., N f in the context of Poisson's, linear elasticity, and Stokes' equations.Note that when the linear elasticity equations are considered, then u, u 0 , u k , and u d are vector-valued functions that we will also, respectively, write u, u 0 , u k , and u d .Similarly, for Stokes equations, u, u 0 , u k , and u d also include the pressure, i.e., we will write them (u, p), (u 0 , p 0 ), (u k , p k ), and (u d , p d ), respectively.
For each case, we will then define an a posteriori estimator E (u d ) of the defeaturing error in the energy norm (Sections 3 and 4), in the case in which Neumann boundary conditions are imposed on the boundaries of the features.This assumption is formalized as follows: Then for each one of the differential problems P(Ω) that will be considered in this article, the introduced a posteriori defeaturing error estimator E (u d ) verifies the two following Theorems 2.5 and 2.7.
Theorem 2.5 (Reliability) Under Assumptions 2.1 and 2.4, the defeaturing error estimator E (u d ) is reliable, meaning that it is an upper bound for the defeaturing error in the energy norm: Before stating Theorem 2.7, let us just introduce the following definition.Definition 2.6 Let Λ be an (n − 1)-dimensional subset of R n .We say that Λ is regular if Λ is piecewise smooth and shape regular, that is, if there exists N Λ ∈ N \ {0} such that for all ℓ 1 , ℓ 2 = 1, . . ., N Λ with where the hidden constant may depend on N Λ but not on the measure of each component Λ ℓ1 , and Λ ℓ1 is smooth.
Theorem 2.7 (Efficiency) Under Assumptions 2.1 and 2.4, if the boundaries of the features are shape regular as in Definition 2.6, then the defeaturing error estimator E (u d ) is efficient up to oscillations, meaning that it is a lower bound for the defeaturing error in the energy norm: where osc(u d ) is a higher order term with respect to the size of the features.
We remark again that the inequalities in Theorems 2.5 and 2.7 do not depend neither on the size nor on the number of features.In addition, as it will be shown in Sections 3 and 4 for Poisson, Stokes, and linear elasticity problems, the estimator E (u d ) is simple and computationally cheap, and it is not only driven by geometrical considerations, but also by the PDE at hand.
The explicit expression of the oscillation term osc(u d ) will be provided in Appendix A, together with the proofs of Theorems 2.5 and 2.7, in the framework of Stokes' equations, as the oscillations corresponding to Poisson's and linear elasticity equations are very similar and their details can be found in [24].The key issue in the analysis is to track the dependence of all constants from the sizes of the features and from their number.

Some further geometric notation
To be able to correctly define the defeaturing error estimator E (u d ), we need to introduce some further notation identifying specific pieces of boundaries of the features.In the following, we use the upper index k to refer to the feature F k , with k = 1, . . ., N f , and the lower indices n and p to refer to negative and positive components of the features, respectively.Now in particular, we let γ be the union of the pieces of boundaries of Ω that are removed by defeaturing, and we let γ 0 be the union of the pieces of boundaries of Ω 0 replacing them, that is, and as illustrated in Figure 3, Using this notation and following the discussion in Section 2.1, we can now precisely determine which simple extension F k p of the positive component F k p can be chosen, for all k = 1, . . ., N f (see Figure 5).More precisely, for the solution of P F k p restricted to F k p to correctly define an extension of the solution of Note that it is possible to have F k p ∩ Ω 0 ̸ = ∅.Let us also define We remark that for all k = 1, . . ., N f , one can look at F k p as the defeatured geometry of the positive component F k p , that is, as a geometry simplified from the exact geometry F k p , for which G k p is a negative feature (see Figure 5c).To simplify the following exposition, and even if this hypothesis could easily be removed, let us make the following assumption regarding these domain extensions.Assumption 2.8 Let us assume that In addition, let γk := ∂ F k p \ ∂F k p , and let γ k p be decomposed as This notation is illustrated in Figure 5.Then, similarly to the previously introduced notation, let In the sequel, we will see that the boundaries γ n , γ 0,p , and γ r will play an important role in the definition of the defeaturing error estimators E (u d ).Therefore, let us also introduce the following notation: and the sets Finally, let n and n 0 respectively denote the unit outward normal vectors to Ω and Ω 0 .Moreover, for all k = 1, . . ., N f , let ñk be the unitary outward normal to F k p , and let n k ≡ n F k denote the unitary outward normal vector to F k n and to F k p .Note that the vectors n k may not be uniquely defined if the outward normal to F k n is of opposite sign of the outward normal to F k p (see Figure 5 for instance), but we allow this abuse of notation since the context will always make it clear.
With this technical notation at hand, we can now precisely restate condition (b) of Assumption 2.1.

Assumption 2.1 (b) There exist sub-domains Ω
where the hidden constant is independent of the size of the features, i.e., the measure of Ω k is comparable with the measure of Ω, not with the measure of the feature F k , domains Ω k is limited and notably smaller than the total number of features N f .This condition is illustrated in Figure 6.Note that if N f = 1, one can take Ω 1 := Ω.

Defeaturing in Poisson's equation
In this section, we consider the Poisson equation and precisely define the exact problem P(Ω), the defeatured problem P(Ω 0 ), and the extension problems P F k p for all k = 1, . . ., N f .Then in this context, we give a precise definition of the proposed reliable and efficient a posteriori estimator E (u d ) of the energy norm of the defeaturing error |||u − u d ||| Ω .
In the following, we denote by H s (D) the Sobolev space of order s ∈ R in a domain D ⊂ R n , whose classical norm and semi-norm are written ∥ • ∥ s,D and | • | s,D , respectively.We will also denote by L 2 (D) := H 0 (D).Moreover, if we let φ ⊂ ∂D and z ∈ H

Exact and defeatured problems
Let us first introduce Poisson's problem P(Ω) in the exact geometry Ω.To do so, let Note that one order of regularity more than usual is required.As we will further see, to define the proposed defeaturing error estimator, we need a solution whose normal derivative along the features' boundaries belongs to L 2 .Then, the problem reads: find u ∈ H 1 g D ,Γ D (Ω), the weak solution of If H 1 0,Γ D (Ω) is equipped with the L 2 (Ω)-norm of the gradient ∥∇•∥ 0,Ω , then we know from Riesz representation theorem that this problem is well-posed.However, note that this problem is usually not solved in practice, as it is assumed to be computationally expensive.
Let us therefore introduce the corresponding Poisson problem P(Ω 0 ) in the defeatured geometry Ω 0 .To do so, we need to consider an L 2 -extension of the restriction f | Ω⋆ in all the negative components F k n , k = 1, . . ., N f , that we still write f ∈ L 2 (Ω 0 ) by abuse of notation.Note that such an extension is not needed in the positive components of the features.Then instead of ( 7), we solve the following defeatured Poisson problem: after choosing for all v 0 ∈ H 1 0,Γ D (Ω 0 ).Let us recall that, according to Assumption 2.4, Γ D ∩ γ = ∅.From the Riesz representation theorem, we know that problem (8) is well-posed.
Remark 3.1 The best possible choices for the defeatured problem data f in F k n for k = 1, . . ., N f and g 0 on γ 0 will be guided by Remark 3.2.However, remark that in applications, the defeatured problem data are rarely chosen, being rather determined by the problem at hand, and being often taken as the natural extension of the exact problem data.For instance, if the right hand side f corresponds to the gravity field, then its only possible physical extension to the negative components of the features is still the gravity field itself.Moreover, one of the aims of defeaturing relies in avoiding the local meshing of the features when the differential problem is solved numerically.Therefore, when the geometry is defeatured, one does not want to have to do a local meshing to capture the behavior of the defeaturing problem data, as this would be equivalent to meshing the feature.We only anticipate here that the best possible defeatured problem data allow the conservation of the solution flux in the positive and negative components of the features.That is, the defeaturing error will be smaller if the defeatured problem data verify for all k = 1, . . ., N f , ˆγk Now, we need to extend the solution u 0 of (8) to F k p , for k = 1, . . ., N f , as discussed in Section 2.1, where F k p satisfies the properties given in (4).To do so, we need to consider an L 2 -extension of the restriction p , that we still write f by abuse of notation.Then, we solve the following extension As before, from Riesz representation theorem, we know that problem ( 10) is wellposed.Once again, the choice of the defeatured problem data f in G k p and gk on γk , for k = 1, . . ., N f , will be guided by Remark 3.2.Thus, as in ( 8), the best possible choices will allow for the conservation of the solution flux in every feature extension G k p .That is, the best possible choices verify for all k = 1, . . ., N f , Then, the defeatured solution u d ∈ H g D ,Γ D (Ω) is defined from u 0 and u k for k = 1, . . ., N f as in (3), and the energy norm of the defeaturing error is defined as

Defeaturing error estimator
In this section, we generalize the defeaturing error estimator introduced in [1], in the case of a geometry that presents multiple features.As for the single feature case, the derived estimator is an upper bound and a lower bound (up to oscillations) of the energy norm of the defeaturing error.
Let us recall the definition of the defeatured solution u d from (3), and the definitions of Σ, Σ k , Σ n , Σ 0,p , and Σ r from (6).Then for all γ ∈ Σ, let k γ ≡ k if γ ∈ Σ k for some k = 1, . . ., N f , and let In other words, d γ represents the Neumann error on the boundaries γ n and γ r due to the defeaturing process, and the jump of the normal derivative of the defeatured solution on the boundaries γ 0,p .Then, if we let η ∈ R be the unique solution of η = − log(η), and if for all γ ∈ Σ, we let we define the a posteriori defeaturing error estimator as: where, for all γ ∈ Σ, Analysis-aware defeaturing of complex geometries with Neumann features where d γ γ denotes the average value of d γ over γ.Note that we can rewrite the estimator feature-wise as follows: where for all k = 1, . . ., N f , we define E k (u k ) as the defeaturing error estimator for feature F k , that is, .
The proposed estimator indicates that the whole information on the error introduced by defeaturing multiple features, in the energy norm, is encoded in the boundary of the features, and can be accounted by suitably evaluating the error made on the normal derivative of the solution.This result generalizes the one from [1].
As a consequence, if these terms dominate, this means that the defeatured problem data should be more accurately chosen, namely g 0 , gk , and the extension of f to G k p .Moreover, under the reasonable flux conservation assumptions (9) and (11), the defeaturing error estimator (15) rewrites However, when n = 2, and under the flux conservation conditions ( 9) and (11), Ẽ (u d ) is sub-optimal since in this case, Ẽ (u d ) ≲ max γ∈Σ c γ E (u d ).Indeed, no lower bound can be provided for Ẽ (u d ).

Defeaturing in linear elasticity and in Stokes' equations
In this section, we consider the linear elasticity and the Stokes equations, and, following the same structure as for the Poisson's problem in Section 3, we define precisely the exact P(Ω), defeatured P(Ω 0 ), and extension problems P F k p for all k = 1, . . ., N f .Then in this context, we give a precise definition of the proposed a posteriori estimator of the energy norm of the defeaturing error, together with the proof of its reliability and efficiency (up to oscillations).
In the following, we denote by H s (D) := [H s (D)] n the vector-valued Sobolev space of order s ∈ R in a domain D ⊂ R n , whose classical norm and semi-norm are again written ∥ • ∥ s,D and | • | s,D , respectively.We will also denote by L 2 (D) := H 0 (D).Moreover, if we let φ ⊂ ∂D and z ∈ H 1 2 (φ), and to be able to deal with boundary conditions, we will denote by where tr φ (y) denotes the trace of y on φ.

Exact and defeatured problems
Let us first introduce the Stokes problem P(Ω) in the exact geometry Ω.To do so, considering a function v : Ω → R n , let ε(v) := 1 2 ∇v + ∇v T be the linearized strain rate tensor in Ω and let be the viscous stress tensor of the considered Newtonian fluid, where the constant µ > 0 is the dynamic viscosity.Note that σ(v) is the viscous stress tensor and not the total Cauchy stress tensor that would be defined by σ(v, q) := σ(v) − q I I I n for some function q : Ω → R in the space of pressures.Assumption 4.1 Hereinafter, for the sake of simplicity in the exposition, we assume that λ is constant everywhere, and it is therefore naturally extended to the defeatured geometry Ω 0 .
, and f ∈ L 2 (Ω).We are interested in the following Stokes problem defined in the exact geometry Ω: find (u, p) ∈ H 1 g D ,Γ D (Ω) × L 2 (Ω), the weak solution of , where for all v, w ∈ H 1 (Ω) and all q ∈ L 2 (Ω), If we equip H 1 0,Γ D (Ω) with the norm ∥∇ • ∥ 0,Ω , it is possible to show by Ladyzhenskaya-Babuška-Brezzi theorem [33] that problem (18) is well-posed.However, note that this problem is never solved in practice, as it is assumed to be computationally too expensive.Remark 4.2 In a more general setting, the viscous stress tensor writes instead of (17), where µ > 0 and λ ≥ 0 are the dynamic and bulk viscosities, respectively, and I I I n is the identity tensor in R n×n .However, note that in this case, • if f c ≡ 0, then (18) is the system of equations describing a linear elastic problem in the incompressible limit λ → ∞, and the constitutive relation (20) then simplifies to (17); we use the change of variables, p ′ := p − λ∇u and σ ′ (u) := 2µε(u) as in (17), then problem (18) remains identical if we replace p by p ′ and the general expression σ(u) from ( 20) by σ ′ (u) from (17).
Analysis-aware defeaturing of complex geometries with Neumann features Hence, instead of (20), we choose the simpler expression (17) for σ, without loss of generality.

Remark 4.3
The linear elasticity equations can be obtained from Stokes' equations ( 18) by removing the pressure terms and the divergence condition.In this case, σ(v) is the Cauchy stress tensor in the medium, satisfying Hooke's law, and in the general form (20), λ and µ denote the first and second Lamé coefficients, respectively, with µ > 0 and λ + 2 3 µ > 0. Since the corresponding defeatured problem and the derivation of the defeaturing error estimator for the linear elasticity equations can simply be obtained by removing the pressure terms everywhere, then we will not detail this case and the interested reader is referred to [24].Now, let us introduce the corresponding Stokes problem P(Ω 0 ) in the defeatured geometry Ω 0 .To do so, we need to choose an L 2 -extension of f and an L 2 -extension of f c in the negative components of the features F n , that we still write f and f c by abuse of notation.Moreover, we assume that the viscous stress tensor σ also satisfies (17) on functions defined everywhere in Ω 0 .Then instead of the exact problem (18), the following defeatured problem is solved: after choosing that is, (u By Ladyzhenskaya-Babuška-Brezzi theorem, problem ( 22) is well-posed.
Remark 4.4 As for the Poisson problem, and even though this is rarely a choice in applications, the best possible choices for the defeatured problem data f and f c in F k n for k = 1, . . ., N f and g 0 on γ 0 will be guided by Remark 4.5.We anticipate that the defeaturing error will be smaller if the defeatured problem data satisfy a conservation assumption of the solution flux in the positive and negative components of the features.I.e., if they verify for all k = 1, . . ., N f , ˆγk Now, for all k = 1, . . ., N f , we need to extend the solution (u 0 , p 0 ) of ( 22) to F k p as discussed in Section 2.1, where F k p satisfies the properties given in (4).To do so, let us choose an L 2 -extension of the restriction p , that we still write f and f c by abuse of notation.Moreover, we assume that the viscous stress tensor σ also satisfies (17) on functions defined in F k p .Then we define for all k = 1, . . ., N f the following extension of the solution (u 0 , p 0 ) of (21) Analysis-aware defeaturing of complex geometries with Neumann features By Ladyzhenskaya-Babuška-Brezzi theorem, problem (25) is well-posed.The choice of the defeatured problem data f and f c in G k p , and gk on γk for k = 1, . . ., N f will be guided by Remark 4.5.We anticipate here that as before, the best possible choices satisfy a conservation assumption of the solution flux in every feature's extension G k p .I.e., they verify for all k = 1, . . ., N f , Then, the defeatured solution (u is defined from (u 0 , p 0 ) and (u k , p k ) for k = 1, . . ., N f as in (3) and we define the defeaturing error as

Defeaturing error estimator
Let us recall the definition of the defeatured solution (u d , p d ) from ( 3), and the definitions of Σ, Σ k , Σ n , Σ 0,p , and Σ r from (6).Then for all γ ∈ Σ, let k γ ≡ k if γ ∈ Σ k for some k = 1, . . ., N f , and let us redefine d γ for all γ ∈ Σ in the context of Stokes equations.That is, let In other words, and as for Poisson's equations, d γ represents the Neumann error on the boundaries γ n and γ r due to the defeaturing process, and the jump of the normal derivative of the defeatured solution on the boundaries γ 0,p .Then, we define the a posteriori defeaturing error estimator as: where for all γ ∈ Σ, where c γ is defined in ( 14), d γ γ denotes the dimension-wise average value of d γ over γ, and ∥ • ∥ ℓ 2 denotes the discrete ℓ 2 -norm as d γ γ ∈ R n .Note that, as for Poisson's problem, we can rewrite the estimator feature-wise as where for all k = 1, . . ., N f , we define E k (u k , p k ) as the defeaturing error estimator for feature F k , that is, .
Analysis-aware defeaturing of complex geometries with Neumann features However, when n = 2 and under the flux conservation conditions (23) and (26) Remark 4.7 If the linear elasticity equations are considered instead of Stokes' equations, with λ and µ denoting the Lamé coefficients, then the pressure terms in the Neumann error (28) should be removed, and the weight µ − 1 2 in equation ( 30) should be replaced by ρ − 1 2 , where ρ = µ if λ ≥ 0, and ρ = min µ, 3 2 λ + µ otherwise.The coefficient ρ corresponds to the coercivity constant of the bilinear form a(•, •), see Appendix A. Then equations ( 29) and (31) defining the estimator remain the same for the linear elasticity equations.

An adaptive geometric refinement strategy
In this section, we aim at defining an adaptive analysis-aware defeaturing strategy in a geometry Ω containing N f ≥ 1 distinct complex features.More precisely, starting from a fully defeatured geometry Ω 0 , we want to precisely define a strategy that determines when and which geometrical features need to be reinserted in the geometrical model, among those that have been removed by defeaturing.Note that the word defeaturing may be misleading when thinking of an adaptive strategy: The geometry Ω 0 in which the problem is actually solved is (partially) defeatured, but the adaptive algorithm selects the features that need to be added to the geometrical model, in order to solve the differential problem up to a given accuracy.The concept of geometric adaptivity is illustrated in Figure 7.
In the sequel, we elaborate on each of the building blocks which compose one iteration of an iterative process: To do so, let i ∈ N be the current iteration index of the adaptive geometric refinement strategy.For simplicity in this section, let us always write u d the defeatured solution, even in the context of linear elasticity for which it should be u d , or in the context of Stokes equations for which it should be (u d , p d ).To begin the process, let Ω (0) 0 be the fully defeatured geometry defined as in (2).That is, Ω (0) 0 is the domain in which all features of Ω are removed: Their positive component is cut out, and their negative component is filled with material.Since some features will be reinserted during the adaptive process, we denote Ω (i) 0 the simplified geometry at the i-th iteration, and in general, we use the upper index (i) to refer to objects at the same iteration.However, to alleviate the notation, we will drop the index (i) when it is clear from the context.In particular, we will write Ω 0 ≡ Ω (i) 0 .

Solve and Estimate
We first solve the defeatured problem (8)/( 22) defined in the (partially) defeatured geometry Ω 0 .Then, we solve the local extension problem (10)/(25) for each feature having a non-empty positive component.We thus obtain the defeaturing solution u d ≡ u (i) d defined in (3), as an approximation of the exact solution u of (7)/(18) at iteration i.Then, the defeaturing error is estimated by E (u d ) defined in ( 15)/(29).

Mark
Recalling that N f ≡ N (i) f at the current iteration i, we select and mark some features F km km∈Im ⊂ F with I m ⊂ 1, . . ., N (i) f to be added to the (partially) defeatured geometry Ω 0 ≡ Ω (i) 0 .To do so, we employ in the following a maximum strategy, but any other convergent marking technique such as the Dörfler strategy [34] could be used.That is, let us first recall definition ( 16)/(31) of the single feature contributions E k (u d ) of the defeaturing error estimator E (u d ), for k = 1, . . ., N f .Then, after choosing a marking parameter 0 < θ ≤ 1, a feature F km is marked, i.e., k m ∈ I m , if it verifies In other words, the set of marked features are the ones giving the most substantial contribution to the defeaturing error estimator.The smallest is θ, the more features are selected, and viceversa.

Refine
In this step, the defeatured geometry Ω (i) 0 is refined, meaning that the marked features F k k∈Im are inserted in the geometrical model.That is, the new partially defeatured geometrical model Ω (i+1) 0 at the next iteration is built as follows: And thus in particular, , and as in definition (2), .
Once the mesh and the defeatured geometry have been refined, the modules SOLVE and ESTIMATE presented in Section 5.1 can be called again.To do so, we update Ω 0 as Ω k∈Im , and renumber the features from 1 to . The adaptive loop is continued until a certain given tolerance on the error estimator E (u d ) is reached, or until the set F is empty, meaning that all the features have been added to the geometrical model.Remark 5.1 Note that a more precise geometric refinement strategy could be performed since G k p can be seen as a negative feature of F k p whose simplified domain is F k p , for all k = 1, . . ., N f .More precisely, one could consider separately the contributions to E k (u d ) given by • γ k n and γ k 0,p , which indicate whether feature F k should be added to the defeatured geometrical model Ω 0 ; • γ k r , which indicates whether the negative feature G k p of F k p should be removed from the simplified positive component F k p of F k .However, since this adds an extra complexity without introducing new conceptual ideas, this strategy is not further developed in the remainder of this article.Remark 5.2 For a given a design, engineers often consider more than one set of boundary conditions, or in other words, they analyze the same design under different loading scenarios.Then, to decide whether the design is valid, they consider an envelope of the results thanks to a max-like function.The proposed adaptive strategy could be adapted to this context by computing the defeaturing error estimator (15)/(29) for each loading scenario.Then, one marks the union of the sets of features marked in each case, as the high-fidelity model used for the analysis should include the features that are important in all the considered loading scenarios.

Numerical considerations and experiments
In this section, we perform a few numerical experiments to illustrate the validity of the proposed a posteriori estimators of the defeaturing error, introduced in Sections 3 and 4. Thanks to these experiments, we also demonstrate that the adaptive procedure presented in Section 5 ensures the convergence of the defeaturing error in the energy norm.
For the numerical approximation of the differential problems treated in this section, we use isogeometric analysis (IGA) on very fine meshes, and multipatch and unfitted boundary techniques for the geometrical description of the features.More specifically, a code has been developed on top of GeoPDEs [35], an opensource and free Octave/MATLAB package for the resolution of PDEs using IGA.The local meshing process required for the integration of trimmed elements uses the in-house tools presented in [36,37], that have been linked to GeoPDEs.Finally, the 3D numerical experiment of Section 6.3 has been performed using an inhouse C++ library for immersed problems that implements the folded decomposition technique presented in [38].The interested reader is referred to [27] for a presentation of isogeometric analysis and advanced spline technologies.
It is important to remark that, even if in this work we adopted spline based discretizations, the estimator and methodology presented in this work are discretization agnostic, being possible to use other techniques.

Impact of some feature properties on the defeaturing error
In this section, we study the impact of some properties of the geometrical features on the defeaturing error and estimator.In particular, we study the influence of the size and shape of the features, of the distance between them, and of their number.

Size of the features
Let us consider a numerical experiment first studied in [1], in which a computational domain contains a very small but important feature, and a large feature whose presence or absence does not affect much the solution accuracy.More precisely, let Ω 0 := (0, 1) 2 be the fully defeatured geometry, and let Ω : where F 1 is a circular hole of radius 10 −3 centered at (1.1 • 10 −3 , 1.1 • 10 −3 ), and F 2 is a larger circular hole of radius 10 −1 centered at (8.9 • 10 −1 , 8.9 • 10 −1 ).The considered geometry is illustrated in Figure 8a.We consider Poisson's problem (7) in the exact computational domain Ω, and its defeatured version (8) defined in Ω 0 , with f (x, y) := 128e −8(x+y) in Ω 0 , g D (x, y) := e −8(x+y) on Γ D := [0, 1) × {0} ∪ {0} × [0, 1) , the bottom and left sides, g(x, y) := −8e −8(x+y) on ∂Ω 0 \ Γ D , and finally g ≡ 0 on ∂F 1 ∪ ∂F 2 .The exact and defeatured solutions u and u 0 of these Poisson's problems have a high gradient in the region around the small feature F 1 , while they are almost identically equal to zero in the region around the large feature F 2 .We then perform the same test, but with square holes instead of circular ones.That is, we consider the same defeatured geometry Ω 0 = (0, 1) 2 , and the same data to solve Poisson's equation (7), but now Ω := Ω 0 \ (F 1 ∪ F 2 ), where F 1 and F 2 are squares centered at (1.1 • 10 −3 , 1.1 • 10 −3 ) and (8.9 • 10 −1 , 8.9 • 10 −1 ), respectively, and whose sides have length 2 • 10 −3 and 2 • 10 −1 , respectively.The geometry is illustrated in Figure 8b, and as before, the solution has a high gradient close to the bottom left corner where F 1 is located, and is almost constantly equal to zero close to the top right corner where F 2 is located.
The values of the defeaturing error estimator (15) and of the defeaturing error |||u − u d ||| Ω are reported in Table 1.In both geometries, independently of the shape of the features, F 1 is indeed more important than F 2 since the estimator for F 1 is four orders of magnitude larger than the estimator for F 2 .This result was expected because of the solution's very high gradient close to F 1 , and because of the homogeneous boundary conditions imposed on the feature's boundaries.In both cases, the proposed estimator well estimates the defeaturing error since the effectivity index is reasonably low, with values comparable to the single feature experiments performed in [1].These results perfectly agree with the theory developed in  Section 3, and illustrate how the proposed estimator does not only depend on geometrical considerations, but also on the analysis behind, i.e., on the considered differential problem that one wants to solve.Hence the name analysis-aware defeaturing.

Distance between features
The following numerical example is used to show that the separability Assumption 2.1 is very weak, as one can consider features that are arbitrarily close to one another, as soon as the number of close features is bounded.Indeed, consider a geometry with either two square features, one positive and one negative, or one complex feature, as follows.Let Ω 0 := (0, 1) 2 , let δ ∈ (−0.1, 0.8), and let ) and F 2 δ := 0.5 + δ 2 , 0.6 + δ 2 × (0.9, 1) , as illustrated in Figure 9.That is, • if δ ≤ 0, then F 1 δ ∪ F 2 δ needs to be considered as a single feature because of Assumption 2.1, where F 1 δ is the positive component of that feature, and F 2 δ is its negative component.In this case, we let γ 1 0,δ := γ 0,p , γ 2 0,δ := γ 0,n , γ 1 δ := γ p , and γ 2 δ := γ n ; • if δ > 0, then F 1 δ and F 2 δ are two distinct features satisfying Assumption 2.1 and separated by a distance δ, where F 1 δ is positive, and F 2 δ is negative.In this case, we let Let us consider Poisson problem (7) with f ≡ 0 in Ω, g D (x, y) := 40 cos(πx)+10 cos(5πx) on Γ D := (0, 1)×{0}, and g ≡ 0 on Γ N := ∂Ω δ \ Γ D .We solve the defeatured Poisson problem (8) with the same data, and we take g 0 ≡ 0 on Finally, we solve the Dirichlet extension problem (10) in F 1 δ = F 1 δ .We choose different values of δ in order to consider different cases: • with δ = 2 • 10 −1 , the distance between the features and the distance between γ 1 0,δ and γ 2 0,δ are of the same order of magnitude as the measures of γ 1 0,δ and γ 2 0,δ ; δ • with δ = 2 • 10 −4 , the distance between γ 1 0,δ and γ 2 0,δ is several orders of magnitude smaller than the measures of γ 1 0,δ and γ 2 0,δ ; • with δ = 0, the boundaries of the feature components intersect in one single point; • with δ = −1 • 10 −3 , the measure of the intersection between the boundaries of the feature components is several orders of magnitude smaller than the measures of the boundaries of the features; • with δ = −9.9• 10 −2 , the measure of the intersection between the boundaries of the feature components is of the same order of magnitude as the measures of the boundaries of the features.
The results are presented in Table 2, and we indeed see that the defeaturing estimator approximates well the defeaturing error in all the different presented cases.In particular, we observe that the effectivity index is not influenced by the distance separating the positive and negative components of the feature(s).This confirms the fact that Assumption 2.1 in not very restrictive in practice.

Number of features
Under Assumption 2.1, the effectivity index of the defeaturing error estimator should not depend on the number of features that are present in the original geometry Ω.To verify this, let Ω 0 := (0, 1) 2 be the fully defeatured domain, and let Ω := Ω 0 \ N f k=1 F k , where N f = 27, and the features F k are circular holes of random radii in the interval (0, 0.01) which are randomly distributed in Ω 0 , under the condition that Assumption 2.1 is satisfied.For the sake of reproducibility, the values of the radii and centers of the features are reported in Table 3.The exact domain Ω with all the 27 features is represented in Figure 10a.
We want to find a good approximation of the solution of Poisson's problem (7) in Ω, whose exact solution is shown in Figure 10c, being f (x, y) := −18e −3(x+y) in Ω 0 , g D (x, y) := e −3(x+y) on the bottom and left boundaries, i.e., on Γ D := [0, 1) × {0} ∪ {0} × [0, 1) , g(x, y) := −3e −3(x+y) on ∂Ω 0 \ Γ D , and g ≡ 0 on ∂F k for k = 1, . . ., N f .Thus we perform the adaptive algorithm introduced in Section 5 starting from the fully defeatured domain Ω (0) 0 := Ω 0 , with marking parameter θ = 0.95.I.e., at every iteration we include the 5% of the features whose error contributes the most.We recursively solve the partially defeatured problem (8) in Ω (i) 0 at each iteration i ≥ 0, and we call u (i) 0 its solution.The results are presented in Figure 11, and the sets of added features at each iteration are the following: {1}, {2, 6}, {4, 16}, {8}, {3}, {5}, {13}, {12}, {17}, {22}, {11}, {10}, {15}, {23}, {24, 16}, {7}, {20, 27}, {18, 25}, {21}, {14}, {19, 9}.For instance, the error is divided by 10 when 9 out of the 27 features are inserted in the partially defeatured geometrical model, i.e., a third of total number of features; this happens at iteration i = 7, and Ω (7) 0 is represented in Figure 10b.On the other hand, features 19, 20, 24, and 25, that are placed near the right top corner, where the solution is gradient is very small, are only activated in the last iterations.We remark that the iteration index is directly linked to the number N (i) f at each iteration i, that is, to the number of features that are still missing in the simplified geometrical model Ω      Feature index k the 27 features in Ω.Moreover, we can see that the effectivity index is independent of the number of features that are not in the simplified geometrical model in which the problem is solved.Indeed, η eff remains almost constant at each iteration, between 2.1 and 3.8.This result perfectly agrees with the theory developed in this paper, in particular the reliability and efficiency results of Theorems 2.5 and 2.7.

Shape of the features
The presence of a feature in the computational domain may greatly but locally perturb the solution, for instance because of a sharp or even re-entrant corner.When such feature is removed, the defeatured solution becomes smoother.In this section, the presented numerical illustration demonstrates that the proposed estimator is able to capture the correct behavior of the error independently of the shape of the feature, and with a low effectivity index.
In this experiment, we aim at finding the exact solution (u, p) of the Stokes' problem (18) defined in Ω, where f (x, y) := exp 4 (x − 0.5) 2 + (y − 0.5) 2 and f c ≡ 0 in Ω, g D ≡ 0 in ∂Ω 0 , and g ≡ 0 in ∂F k for all k = 1, 2, 3. Instead of solving problem (18), we tackle the approximate one (21) in Ω 0 with the same data; in particular, f is naturally extended in the features, and we obtain the defeatured solution (u 0 , p 0 ).
The magnitude of the velocity fields u and u d is shown in Figure 12.As it can be seen in the zoomin (Figure 12b), in this case, the sharp non-convex feature F 3 introduces a localized perturbation that is similar to ones produced by the sharp convex F 1 and smooth F 2 features, as all of them are located at solution regions with similar gradient magnitudes.Indeed, their contributions E k (u d , p d ), for k = 1, 2, 3, reported in Table 4, are of the same order.From those individual contributions, through (31), the total defeaturing error estimate is computed to be E (u d , p d ) = 1.0895 • 10 −1 , and the corresponding energy error is |||(u − u d , p − p d )||| Ω = 7.5381 • 10 −2 .Therefore, the effectivity index of the proposed estimator is which is notably very low (at the same level as the ones observed in [1] for geometries with a single feature): The estimator is able to estimate the effect of different features, independently of their shape.

Lid-driven cavity
For this next numerical experiment, let us consider Stokes' problem in a lid-driven cavity [39] in which three holes are located.More precisely, we consider an exact domain Ω = Ω 0 \ (F 1 ∪ F 2 ∪ F 3 ) for which Ω 0 := (0, 1) 2 is the fully defeatured domain, and F 1 , F 2 , and F 3 are three circular holes of radius 0.01 centered, respectively, at (0.011, 0.989), (0.5, 0.75), and (0.011, 0.011).For this test, we let Γ D := ∂Ω 0 , Γ N := ∂F 1 ∪ ∂F 2 ∪ ∂F 3 , f ≡ 0 in Ω 0 , g ≡ 0 on Γ N , and (1, 0) on (0, 1) × {1} the top boundary, (0, 0) everywhere else on Γ D .Results are presented in Figure 13, in which the magnitude of the velocity fields u and u 0 is shown, and Table 5 reports each feature contribution's E k (u d , p d ) for k = 1, 2, 3. We observe that the presence of the features changes the fluid velocity inside the cavity, especially for features F 2 and F 3 .This is expected because of the homogeneous Neumann boundary conditions on the features.However, features only change locally the velocity field of the fluid around them.Thus, feature F 1 brings the largest contribution to the estimator while feature F 3 brings the smallest one.This is coherent with the theory since the velocity has a large gradient close to the moving boundary at the top, while it is almost constantly equal to zero around features F 2 and F 3 (note the different scales in Figures 13b-13d).We therefore expect F 1 to contribute the most to the overall defeaturing error, even if at a first glance it is the hole whose absence changes the least the velocity of the fluid (recall Figure 13b).
which is notably very low, as in the previous experiment.Note that to compute the defeaturing error, a very fine mesh around the holes had to be taken in Ω in order to be able to neglect the component of the error coming from the numerical approximation of the problem, as represented in Figure 14.This is not required when the holes are filled as in the defeatured geometry Ω 0 : This shows the potential of defeaturing, in terms of memory and computational time savings.Far from the features and high solution gradients, a coarser grid is considered.

Three-dimensional elastic structure
For this last numerical experiment, let us consider the exact domain Ω and the corresponding defeatured domain Ω 0 represented in Figure 15.More precisely, the base has dimensions 200 × 200 × 20 [mm], and the cylinder has a height of 150 [mm].Moreover, and in particular, the exact domain contains 20 features numbered as illustrated in Figure 15a: • F 1 to F 4 are the four letters of the carved "EPFL" logo, in order (see also Figure 7 in which these features are more clearly visible).
• F 5 to F 8 are the four holes in the stiffeners, counted counter-clockwise beginning from the one on the left of the "EPFL" logo.
• F 9 to F 12 are the four holes in the vertical part of the structure, counted counter-clockwise beginning from the one above the "EPFL" logo.
• F 13 to F 20 are the eight rounds present on the left and right diagonal angles of the stiffeners, counted counter-clockwise beginning from the left round of the stiffener on the left of the "EPFL" logo.
Rounds, holes, and carved logos are three of the most typical features that finite element analysis practitioners encounter in CAD designs.These features are interesting to analyze, since they are usual candidates to be removed before creating a finite element mesh.
Taking the origin at the bottom lower left corner of the structure, let Γ D be the bottom of the structure and Γ and where Γ top is the top face of the cylinder.Then, let u ∈ H 1 0,Γ D (Ω) be the solution of the linear elasticity problem given by (18) in which the pressure terms and the divergence condition are removed, and where the material properties correspond to steel.That is, the Lamé parameters λ and µ are expressed in terms of the Young modulus E = 210 [GPa] and Poisson's ration ν = 0.3 [-] as     Then we compute the defeatured solution u d ≡ u 0 ∈ H 1 0,Γ D (Ω 0 ) given by problem (22) in which again the pressure terms and the divergence condition are removed.Finally, we compute the estimator E (u d ) defined in (31) with p d ≡ 0 by computing each feature contribution E k (u d ) for k = 1, . . ., 20.
A rather fine mesh is used in order to reduce the error derived from numerical approximation.More precisely, the bounding box of Ω 0 is meshed with n el = 128 elements per direction, and B-splines of degree 2 and regularity 1 are used.Results are presented in Table 6, where we report each feature's contribution E k (u d ) for k = 1, . . ., 20.The obtained total error estimator is equal to E (u d ) = 1.716 • 10 −6 [J].Moreover, the magnitude of the solution displacements u and u d , and the corresponding von Mises stress distributions are shown in Figures 16 and 17, respectively.
We can first see that the absence of features F 5 to F 8 in the defeatured geometry significantly affects the solution in the stiffeners.This is indeed reflected in the estimator: The estimator contributions of those four features is very large, corresponding to around half of the total error estimator.On the other hand, the solution is basically constant around the "EPFL" logo, no deformation is observed around it.We can therefore expect that the absence of features F 1 to F 4 in the defeatured geometry is not affecting much the accuracy of the solution.This is indeed observed in the estimator contributions of those features, as E 1 (u d ) to E 4 (u d ) are the lowest contributions of the estimator, corresponding to around 1% − 2% of E (u d ).This is a typical situation that simulation practitioners encounter daily: Carved logos and trademarks are usually defeatured before creating a finite element mesh, since they complicate the meshing process and increase the number of elements (see, e.g., Figure 1a), but they contribute little to the accuracy of the problem's solution.The proposed estimator identifies them straightaway.
Let us now run the adaptive algorithm introduced in Section 5 with θ = 0.99 as marking parameter, until all features are added to the geometrical model.We call u    iteration i.In Table 7, we report the indices of the features that are added to the defeatured geometrical model at each iteration, together with the value of the estimator E u Comparing the values of the estimator at each iteration and the von Mises stress distributions around each feature, we can see that the features that are added to the geometrical model at each iteration seem to be the ones that are affecting the most the solution accuracy, as one would expect.In Table 7 we can also see that to reduce the error estimator by 90%, it is enough to consider 12 out of the total 20 features of Ω (see iteration 4, whose solution is represented in Figure 18).
For instance, the holes F 9 and F 11 are added before the holes F 10 and F 12 during the adaptive process, because of the direction in which the structure is bending due to the applied traction along the x-direction; this is reflected by the variation of the von Mises stresses that are larger in F 9 and F 11 than in F 10 and F 12 .We can also see that larger stresses are present nearer the rounds F 14 , F 15 , F 18 , and F 19 than around the other four rounds.This is again coming from the direction of the bending.And very interestingly, the estimator is able to capture this effect, as rounds F 14 , F 15 , F 18 , and F 19 are introduced in the defeatured geometry after iteration 3, while the other rounds are introduced later, after iterations 4 and 5. See, for instance, the stress distribution in the connection between features F 17 and F 18 and the main cylinder (zoom-in regions in Figures 16,17,and 18).As it can be appreciated, the stress concentration is higher around feature F 18 , fact that also reveals the estimator value in Table 6, and the fact that feature F 18 is activated before than F 17 (see Table 8 and Figure 18).
Rounds are other typical examples of features that are candidates to be removed.However, in this case, the situation is usually less clear.Indeed, on the one hand, rounds complicate the meshing process and increase the number of elements in the model.But on the other hand, depending on the boundary conditions, removing rounds may lead to the creation of singularities in the solution.The proposed estimator is able to determine the impact of removing those rounds.
Finally, the numerical error is not considered in this article.However, it is interesting to note that the  is reported in Figure 19, and the features chosen at each iteration are reported in Table 8.We can observe that except for features whose error contributions are very close to one another, the adaptive algorithm is able to correctly choose the important features, even on the coarsest mesh.

Conclusions
In the context of the Poisson's, linear elasticity, and Stokes equations, we have studied the accuracy impact of removing distinct Neumann features from geometries in which the solution of a PDE is sought.In particular, we have generalized the a posteriori estimator of the energy norm of the defeaturing error from [1] to twoand three-dimensional geometries containing an arbitrary number of negative, positive, or generally complex features.The proposed estimator has the following properties: • it is not only driven by geometrical considerations, but also by the PDE at hand; • it is able to weight the impact of defeaturing in the energy norm, and its effectivity index is independent of the size of the geometrical features and of their number; • it is able to determine whether the defeaturing error comes from the choice of defeaturing data (right hand side and Neumann boundary conditions), or if it comes from the importance of the presence of the feature itself; • it is rigorously proven to be reliable and efficient up to oscillations; • it is naturally decomposed into single feature contributions; • it is simple, computationally cheap, and embarrassingly parallel, as it only requires the evaluation of fluxes through boundary pieces and the resolution of small problems at feature level.
• it has been tested on an extensive set of numerical experiments: In all of them, the estimator acts as an excellent approximation of the defeaturing error.
Note however that our framework does not include the case of a geometry whose boundary is complex everywhere as considered for instance in [29][30][31], because of Assumption 2.1.Then, with the help of the proposed error estimator, we have been able to design an adaptive geometric refinement strategy taking into account the defeaturing errors.More precisely, starting from a fully defeatured geometry, features are iteratively added to the geometrical model when their absence is responsible for most of the solution accuracy loss.That is, the strategy is able to build a (partially) defeatured geometric model containing few features, for which the defeaturing error is below a prescribed tolerance.Presented numerical experiments have demonstrated the convergence of the defeaturing error during the adaptive loop.In a subsequent work [23], the proposed adaptive strategy will be combined with a mesh refinement strategy in the case in which a finite element method (and in particular isogeometric analysis [27,28]) is used to approximately solve the PDE at hand.In a very similar fashion, we can deduce that for all k = 1, . . ., N f and all (v k , q b(e u , q) = − ˆΩ q∇ • e u dx = 0.

A.2 Efficiency
In this section, we prove Theorem 2.7 in the context of Stokes' equations.That is, we prove that the error indicator defined in ( 29) is efficient under Assumptions 2.1 and 2.4, i.e., it is a lower bound for the defeaturing error, up to oscillations.The proof is given in the general case for n = 3, but under the following additional assumption if n = 2: Assumption A.1 If n = 2, assume that the compatibility conditions (23) and (26) are satisfied.Note however that the numerical experiments show that the result also holds without compatibility conditions if n = 2.The same restrictions were considered in the efficiency proofs of [1] in the two-dimensional case.
In the following, for all D ⊂ R n and Λ ⊂ ∂D, we denote by H 1 2 00 (Λ) the space of functions w ∈ L 2 (Λ) such that the zero-extension w ⋆ of w to ∂D belongs to H 1 2 (∂D), and we denote by H − 1 2 00 (Λ) its dual space.For the sake of simplicity, let us strengthen Definition 2.6, even though the proofs could be easily generalized to the weaker setting.Definition A.2 An (n−1)-dimensional subset Λ of R n is regular if Λ is piecewise shape regular and composed of flat elements, that is, if there is N Λ ∈ N such that for all ℓ 1 , ℓ 2 = 1, . . ., N Λ with ℓ 1 ̸ = ℓ 2 , Λ = int That is, we suppose that the boundaries γ ∈ Σ of the features are shape regular as in Definition A.2 instead of Definition 2.6, for simplicity.This allows us to easily define the following Clément operator.Whenever some boundary Λ is regular, then for all m ∈ N, let Q pw m,0 (Λ) be the space of continuous piecewise polynomials of degree at most m on each variable, that vanish at the boundary ∂Λ.Then we define Π m,Λ : L 2 (Λ) → Q pw m,0 (Λ) (46) as the extension of the Clément operator [40] developed in [41] on Λ.
Analysis-aware defeaturing of complex geometries with Neumann features

Figure 1 :
Figure 1: Example of an original (left) and defeatured CAD designs (right).Both finite element meshes were generated with Gmsh [5] using the same mesh algorithm and parameters.

Figure 2 :
Figure 2: Domains with different types of geometrical features F .In each case, the negative component of F is dashed while its positive component is filled in gray.
General complex feature in a domain Ω.

Figure 3 :
Figure 3: Geometries containing different types of features.

Figure 5 :
Figure 5: Domain with one complex feature and illustration of the notation.The different domains are distinguished by different gray intensities.

Assumption 2 . 4
Let us decompose ∂Ω = Γ D ∪ Γ N such that Dirichlet boundary conditions are imposed on Γ D with |Γ D | > 0, Neumann boundary conditions are imposed on Γ N , and Γ D ∩ Γ N = ∅.Then we assume that for all k = 1, .

Remark 3 . 2
As analogously noted in [1, Remarks 4.1 and 5.3], the terms involving the average values of d γ in the estimator E (u d ) only depend on the defeatured problem data since for all k = 1, . . ., N f ,

Remark 4 . 5 . 4 . 6
Similarly to Remark 3.2, the terms involving the component-wise average values of d γ in E (u d , p d ) only depend on the defeatured problem data.As a consequence, if these terms dominate, this means that the defeatured problem data should be more accurately chosen.Moreover, under the reasonable data compatibility conditions (23) and (26) that represent flux conservation assumptions in this context, the defeaturing error estimator (29) rewrites E (u d , p d ) := µ − Remark Similarly to Remark 3.3, note that

F 1 F 2 Ω
(a) Exact domain with two circular features.

F 1 F 2 Ω
(b) Exact domain with two square features.

Figure 8 :
Figure 8: Numerical test 6.1.1 -Exact geometries used for the comparison between features' sizes (not at scale).
Zoom on part of the upper boundary of Ω0 (up) and Ω δ (down) for δ = 0. Zoom on part of the upper boundary of Ω0 (up) and Ω δ (down) for δ = −0.05.

0
with respect to Analysis-aware defeaturing of complex geometries with Neumann features Exact solution in Ω. .

Figure 11 :
Figure 11: Numerical test 6.1.3-Behavior of the defeaturing error and estimator with respect to the number of features in the defeatured geometrical models Ω (i) 0 .Each marker corresponds to the value at one iteration.
(a) Magnitude of u (left) and u d (right) in the domains Ω and Ω0, respectively.(b) Zoom of the central region: Magnitude of u (left) and u d (right) in the domains Ω and Ω0, respectively.

Figure 12 :
Figure 12: Numerical test 6.1.4-Magnitude of the velocity in the exact Ω and fully defeatured Ω 0 domains.

( a )
Magnitude of u (left) and u0 (right) in the domains Ω and Ω0, respectively.(b) Magnitude of u (left) and u0 (right) around F 1 .(c) Magnitude of u (left) and u0 (right) around F 2 .(d) Magnitude of u (left) and u0 (right) around F 3 .

Figure 13 :
Figure 13: Numerical test 6.2 -Magnitude of the velocity in the exact domain Ω and in the corresponding fully defeatured domain Ω 0 .

Figure 14 :
Figure 14: Numerical test 6.2 -Mesh used to compute an overkilled solution of the lid driven cavity Stokes' problem in the exact domain Ω.

d
the solution of the defeatured problem at Analysis-aware defeaturing of complex geometries with Neumann features Exact geometry Ω and numbering of the 20 features (in color).(b) Defeatured geometry Ω0
magnitude of the solution displacement and the corresponding von Mises stress distribution at iteration 4 are represented in Figure18.

Figure 17 :Figure 18 :
Figure 17: Numerical test 6.3 -Exact solution in the exact geometry Ω.The views correspond to those of Figure 15, and the deformed configuration is magnified [×5 • 10 3 ] for visualization purposes.

Table 1 :
Numerical test 6.1.1 -Results of the comparison between features' sizes.

Table 2 :
Numerical test 6.1.2-Values of the defeaturing error and estimator for different values of δ.The cases in which δ > 0 correspond to separate features, while the cases δ < 0 correspond to features with overlapping boundaries.

Table 3 :
Numerical test 6.1.3-Data of the 27 circular features.

Table 4 :
Numerical test 6.1.4-Feature contributions E k (u d , p d ) to the multi-feature estimator E (u d , p d ).

Table 5 :
Numerical test 6.2 -Feature contributions E k (u d , p d ) to the multi-feature estimator E (u d , p d ).The total defeaturing error estimate is equal to E (u d , p d ) = 35.622,and |||(u − u d , p − p d )||| Ω = 31.308.Therefore, the effectivity index of the proposed estimator is equal to

Table 7 :
Numerical test 6.3 -Results of the adaptive defeaturing strategy for n el = 128 elements.

Table 8 :
Numerical test 6.3 -Marked features at each iteration on different mesh refinements.estimator is still able to drive the proposed adaptive strategy on a coarser mesh.Indeed, this algorithm has been performed on multiple meshes containing a different number n el of elements in each space direction of the bounding box of Ω 0 .More precisely, we have considered n el = 8, 32, and 128.In all three cases, the convergence of the estimator E u