# Strict lower bounds with separation of sources of error in non-overlapping domain decomposition methods

Valentine Rey, Pierre Gosselet, Christian Rey

► To cite this version:

Valentine Rey, Pierre Gosselet, Christian Rey. Strict lower bounds with separation of sources of error in non-overlapping domain decomposition methods. International Journal for Numerical Methods in Engineering, 2016, 10.1002/nme.5244. hal-01332674

HAL Id: hal-01332674

https://hal.science/hal-01332674

Submitted on 16 Jun 2016

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>
Strict lower bounds with separation of sources of error in non-overlapping
domain decomposition methods

V. Rey¹, P. Gosselet¹, C. Rey²

¹ LMT-Cachan / Ecole Normale Supérieure de Cachan, CNRS, Université Paris Saclay
61, avenue du président Wilson, 94235 Cachan, France,

² Safran Tech, rue des Jeunes Bois Châteaufort CS 80112,78772 Magny les Hameaux, France

June 16, 2016

Abstract

This article deals with the computation of guaranteed lower bounds of the error in the framework of finite
element (FE) and domain decomposition (DD) methods. In addition to a fully parallel computation, the proposed
lower bounds separate the algebraic error (due to the use of a DD iterative solver) from the discretization error (due
to the FE), which enables the steering of the iterative solver by the discretization error. These lower bounds are
also used to improve the goal-oriented error estimation in a substructured context. Assessments on 2D static linear
mechanic problems illustrate the relevance of the separation of sources of error and the lower bounds' independence
from the substructuring. We also steer the iterative solver by an objective of precision on a quantity of interest.
This strategy consists in a sequence of solvings and takes advantage of adaptive remeshing and recycling of search
directions.

**Keywords:** Verification; Error estimation; Finite element method; Domain decomposition methods; FETI; BDD

# 1 Introduction

Virtual testing is a useful tool for engineers to certify structures without resorting to experimental tests. However, its massive adaption comes with several challenges. Among others, virtual testing requires the capability to solve large problems (several millions degrees of freedom) and to warrant the quality of the results provided by simulations. To tackle these difficulties, we propose to use domain decomposition methods and verification. On the one hand, non-overlapping domain decomposition methods [11, 16, 7] are well-known techniques that enable the solving of large mechanical problems by exploiting parallel computers' performance. On the other hand, verification provides tools to estimate the distance between the unknown exact solution and the computed approximated solution. This distance is the approximation error, it can be estimated by a global energy norm or by local quantities of interest (goal-oriented error estimation).

This paper is the continuation of papers connecting domain decomposition methods and verification. In [22], the authors proposed a parallel error estimator based on the error in constitutive relation [14] in a substructured framework. They described a methodology to construct the required admissible fields for error estimation. Those fields were rebuilt using quantities naturally processed during the solving and preconditioning steps of classical domain decomposition algorithms as inputs for classical equilibration techniques used in parallel on each subdomain. In [28], a new parallel error estimator that separates the discretization error (due to the finite element method) from the algebraic error (due to the iterative solver) was proposed. This new estimator enables the definition of a new stopping criterion for the iterative solver, no longer defined regardless the discretization, which avoids over-solving. Finally, in [27], this work was extended to goal-oriented error estimation. The exact value of a linear quantity of interest defined by an extractor [19, 32, 18] was estimated using global error estimation of a forward problem and an adjoint problem. It was shown that these two problems could be solved simultaneously thanks to a block-Krylov algorithm [30] steered by an objective on the error on the quantity of interest.

Upper bounds are the main concern of verification and the literature on lower bounds is scarce. By exploiting the
residual equation [24] and constructing a continuous error estimation, a lower bound of the error can be computed. The
question of the construction of a continuous error estimation has already been adressed in many papers (for instance
[20, 4, 9]) where the authors benefit the computation of an upper bound for the computation of a lower bound.

The objective of this paper is to investigate the computation of a lower bound of the error in a substructured
framework. In line with the papers associating domain decomposition methods and verification, the lower bounds are
computed in parallel and they separate the two sources of error (discretization error and algebraic error). Based on
the results demonstrated in [20], the computation of the lower bound does not involve significant cost since it exploits
the fields computed during the reconstruction of a statically admissible stress fields.
The paper is organized as follows. In section 2, we define the reference problem and recall the principle of the error in constitutive relation. We also recall the domain decomposition methods' principles and highlight the fields built at each iteration. Finally, we recall the parallel error estimators developed in [21, 28] and give a brief state of the art of the computation of continuous fields for lower bounds of the error. In section 3, two theorems providing lower bounds with and without separation of sources of error are demonstrated. The parallel reconstruction of fields required to compute these bounds is detailed. We also show how to benefit from this global information on the error to better the goal-oriented error estimation. In section 4, we apply these lower bounds on two-dimensional mechanical structures. We compare the lower bounds provided by a sequential approach and primal and dual approaches and also study the independence with respect to the substructuring. We illustrate the convergence of the lower bounds during the iterations and illustrate the separation of sources. Finally, in order to reach an objective of precision on a quantity of interest, we apply an auto-adaptive strategy on one of the structures. In this strategy, we use the separation of sources of error to define the stopping criterion for the iterative solver. Benefiting the informations from a first solving, we process adaptive remeshing to better the FE solution and lower the error bounds and we recycle the search directions generated (Krylov subspace recycling, see [25, 10, 29]) to speed up further solvings. Section 5 concludes the paper.

# 2 Settings

## 2.1 Reference problem

Let $\mathbb{R}^d$ represents the physical space. Let us consider the static equilibrium of a (polyhedral) structure which occupies the open domain $\Omega \subset \mathbb{R}^d$ and which is subjected to given body force $f \in L^2(\Omega)$ within $\Omega$, to given traction force $g \in L^2(\partial_g\Omega)$ on $\partial_g\Omega$ and to given displacement field $u_d$ on the complementary part of the boundary (such that meas($\partial_u\Omega$) $\neq$ 0). We assume that the structure undergoes small perturbations and that the material is linear elastic, characterized by Hooke's elasticity tensor $\mathbf{H}$. Let $u$ be the unknown displacement field, $\underline{\varepsilon}(u)$ the symmetric part of the gradient of $u$, $\underline{\sigma}$ the Cauchy stress tensor. Let $\omega$ be an open subset of $\Omega$.

We introduce two affine subspaces and one positive form:

* Affine subspace of kinematic admissible fields (KA-fields)

$$KA(\omega) = \{\underline{u} \in (\mathbb{H}^1(\omega))^d, \underline{u} = \underline{u}_d \text{ on } \partial\omega \cap \partial_u\Omega\} \quad (1)$$

and we note $KA^0(\omega)$ the following linear subspace:

$$KA^0(\omega) = \{\underline{u} \in (\mathbb{H}^1(\omega))^d, \underline{u} = 0 \text{ on } \partial\omega \cap \partial_u\Omega\} \quad (2)$$

and $KA^{00}(\omega)$ the following linear subspace:

$$KA^{00}(\omega) = \{\underline{u} \in (\mathbb{H}^1(\omega))^d, \underline{u} = 0 \text{ on } \partial\omega \setminus \partial_g\Omega\} \quad (3)$$

*Remark.* Note that if $\omega = \Omega$, $KA^{00}(\omega)$ and $KA^0(\omega)$ are identical.

* Affine subspace of statically admissible fields (SA-fields)

$$SA(\omega) = \left\{ \underline{\tau} \in (L^2(\omega))_{\text{sym}}^{d \times d}; \forall \underline{u} \in KA^{00}(\omega), \int_{\omega} \underline{\tau} : \underline{\varepsilon}(\underline{u}) d\Omega = \int_{\omega} f \cdot \underline{v} d\Omega + \int_{\partial\omega \cap \partial_g\Omega} \underline{g} \cdot \underline{v} dS \right\} \quad (4)$$

* Error in constitutive relation [14]

$$e_{CR_{\omega}}(\underline{u}, \underline{\sigma}) = \|\underline{\sigma} - \mathbb{H} : \underline{\varepsilon}(\underline{u})\|_{H^{-1}, \omega} \quad (5)$$

$$\text{where } \|\underline{x}\|_{H^{-1}, \omega} = \sqrt{\int_{\omega} (\underline{x} : \mathbb{H}^{-1} : \underline{x}) d\Omega}$$

The mechanical problem set on $\Omega$ can be formulated as:

$$\text{Find } (\underline{u}_{ex}, \underline{\sigma}_{ex}) \in KA(\Omega) \times SA(\Omega) \text{ such that } e_{CR_{\Omega}}(\underline{u}_{ex}, \underline{\sigma}_{ex}) = 0 \quad (6)$$

The solution to this problem, named “exact” solution, exists and is unique.
Remark. The formulation (6) is equivalent to the classical following formulation:

$$ \text{Find } \underline{u} \in \text{KA}(\Omega) \text{ such that } \forall \underline{v} \in \text{KA}^{00}(\Omega), a(\underline{u}, \underline{v}) = L(\underline{v}) \quad (7) $$

with

$$ a(\underline{u}, \underline{v}) = \int_{\Omega} \underline{\varepsilon}(\underline{u}) : \mathbb{H} : \underline{\varepsilon}(\underline{v}) d\Omega \quad (8) $$

and

$$ L(\underline{v}) = \int_{\Omega} \underline{f} \cdot \underline{v} d\Omega + \int_{\partial_g \Omega \cap \partial\Omega} \underline{g} \cdot \underline{v} dS \quad (9) $$

### 2.1.1 Finite element approximation

Let us consider a mesh of $\Omega$ to which we associate the finite-dimensional subspace $KA_H(\Omega)$ of $KA(\Omega)$. $\mathcal{T}$ is the set of elements of the mesh and $\mathcal{V}$ is the set of vertexes. The classical finite element displacement approximation consists in searching:

$$ \begin{gathered} \underline{u}_H \in \text{KA}_H(\Omega) \\ \underline{\sigma}_H = \mathbb{H} : \underline{\varepsilon}(\underline{u}_H) \\ \int_{\Omega} \underline{\sigma}_H : \underline{\varepsilon}(\underline{v}_H) d\Omega = \int_{\Omega} \underline{f} \cdot \underline{v}_H d\Omega + \int_{\partial_g \Omega} \underline{g} \cdot \underline{v}_H dS, \quad \forall \underline{v}_H \in \text{KA}_H^0(\Omega) \end{gathered} \quad (10) $$

Of course the approximation is due to the fact that in most cases $\underline{\sigma}_H \notin \text{SA}(\Omega)$.

After introducing the matrix $\underline{\varphi}_H$ of shape functions which form a basis of $KA_H(\Omega)$ (extended to Dirichlet degrees of freedom) and the vector of nodal unknowns $\underline{u}$ so that $\underline{u}_H = \underline{\varphi}_H \underline{u}$, the classical finite element method leads to the linear system:

$$ \begin{pmatrix} \underline{K}_{rr} & \underline{K}_{rd} \\ \underline{K}_{dr} & \underline{K}_{dd} \end{pmatrix} \begin{pmatrix} \underline{u}_r \\ \underline{u}_d \end{pmatrix} = \begin{pmatrix} \underline{f}_r \\ \underline{f}_d \end{pmatrix} + \begin{pmatrix} 0 \\ \underline{\lambda}_d \end{pmatrix} \quad (11) $$

where $\underline{K}$ is the (symmetric semi positive definite) stiffness matrix and $\underline{f}$ is the vector of generalized forces; Subscript $d$ stands for Dirichlet degrees of freedom (where displacements are prescribed) and Subscript $r$ represents the remaining degrees of freedom so that unknowns are $\underline{u}_r$ and $\underline{\lambda}_d$ where Vector $\underline{\lambda}_d$ represents the nodal reactions:

$$ \lambda_d^T = \int_{\Omega} \underline{\sigma}_H : \underline{\varepsilon}(\underline{\varphi}_{Hd}) d\Omega - \int_{\Omega} \underline{f} \cdot \underline{\varphi}_{Hd} d\Omega - \int_{\partial_g \Omega} \underline{g} \cdot \underline{\varphi}_{Hd} dS \quad (12) $$

where $\underline{\varphi}_{Hd}$ is the matrix of shape functions restricted to the Dirichlet nodes and $\underline{n}$ the outer normal vector.

### 2.1.2 A posteriori error estimation

**Upper bound of the discretization error** The estimator we choose is based on the error in constitutive relation, which gives a guaranteed estimator for the discretization error.

The fundamental relation is the following (Prager-Synge theorem, see for instance [15]):

$$ \begin{gathered} \forall (\underline{\hat{u}}, \underline{\hat{\sigma}}) \in \text{KA}(\Omega) \times \text{SA}(\Omega), \\ \|\underline{\varepsilon}(\underline{u}_{ex}) - \underline{\varepsilon}(\underline{\hat{u}})\|_{\mathbb{H},\Omega}^2 + \|\underline{\sigma}_{ex} - \underline{\hat{\sigma}}\|_{\mathbb{H}^{-1},\Omega}^2 = e_{CR\Omega}^2(\underline{\hat{u}}, \underline{\hat{\sigma}}) \end{gathered} \quad (13) $$

We note $\|\underline{u}\|_{\Omega} = \|\underline{\varepsilon}(\underline{u})\|_{\mathbb{H},\Omega}$ the energy norm of the displacement, and since we can choose $\underline{\hat{u}} = \underline{u}_H \in \text{KA}(\Omega)$, we retain the following upper bound for the error $\underline{\varepsilon}_{discr} = \underline{u}_{ex} - \underline{u}_H$:

$$ e_{discr} := \|\underline{\varepsilon}_{discr}\|_{\Omega} \le e_{CR\Omega}(\underline{u}_H, \underline{\hat{u}}) \quad (14) $$

The construction of $\underline{\hat{u}} \in \text{SA}(\Omega)$ is a complex problem solved by various approaches [14, 20, 23, 26].

The techniques [14, 23, 26] are two-steps procedures. The first step consist in building a set of equilibrated tractions or works along the edges of the elements of the mesh. Each method proposes its own strategy to reconstruct such tractions. The second step, common to all methods, is the solving of Neumann problems on each element using the equilibrated tractions as Neumann conditions.

The technique developed in [20] does not require equilibrated fluxes but only the solving local problems on star-patches (a star patch, denoted by $\omega_i$, is composed of the elements sharing the vertex $i$ and corresponds to the support of the shape function associated to the vertex $i$). This is the reason why the technique is sometimes called the flux-free technique.

Local problems (on element or star-patch) are usually solved on a space of finite dimension which is richer than the finite element space restrained to the support of the local problem. The space is enriched either thanks to higher degree polynomial shape functions or thanks to mesh refinement (each element being divided into several smaller elements).
**Lower bound of the discretization error** A lower bound of the true error can be obtained using Cauchy-Schwarz inequality in the following residual equation:

$$ 
\begin{aligned}
& \forall \underline{w} \in \mathrm{KA}^0(\Omega) \\
& \int_{\Omega} \underline{\varepsilon}(\underline{u}_{ex} - \underline{u}_H) : \mathbb{H} : \underline{\varepsilon}(\underline{w}) d\Omega = \int_{\Omega} \underline{f} \cdot \underline{w} d\Omega + \int_{\partial_g \Omega \cap \partial\Omega} \underline{g} \cdot \underline{w} dS - \int_{\Omega} \underline{\sigma}_{\underline{H}} : \underline{\varepsilon}(\underline{w}) d\Omega \\
& := R_H(\underline{w})
\end{aligned} 
\quad (15) $$

Therefore, every displacement field $\underline{w} \in \mathrm{KA}^0(\Omega) \setminus \{\underline{0}\}$ can be used to obtain the following strict lower bound [24]:

$$ \|\underline{u}_{ex} - \underline{u}_H\|_{\Omega} \geq \frac{|R_H(\underline{w})|}{\|\underline{w}\|_{\Omega}} \quad (16) $$

The accuracy of the lower bound depends on the quality of the continuous field $\underline{w}$, which is often called continuous error estimate. Indeed, the lower bound equals the true error for the continuous field $\underline{w} = \underline{u}_{ex} - \underline{u}_H$. The construction of $\underline{w}$ was mainly studied in [4, 20, 9] of which we recall the main results.

In [4], the error is estimated thanks to an implicit residual-based error estimator with local solvings on elements. A continuous field $\underline{w} \in \mathrm{KA}^0(\Omega) \setminus \{\underline{0}\}$ is constructed from the element estimators by averaging on the edges.

In [20], local problems on star-patches are solved to compute an upper bound of the error:

$$ 
\begin{aligned}
& \text{Find } \underline{\varepsilon}^i \in \mathrm{KA}^0(\omega_i) \text{ such that } \forall \underline{v} \in \mathrm{KA}^0(\omega_i) \\
& a(\underline{\varepsilon}^i, \underline{v}) = R_H(\varphi_H^i \underline{v})
\end{aligned} 
\quad (17) $$

where $\varphi_H^i$ is the shape function associated to the central node i of the star-patch. The previous problem is solved on a space $\mathrm{KA}_h^0(\omega_i)$ richer than $\mathrm{KA}_H^0(\omega_i)$.

The proposed continuous displacement field is:

$$ \underline{w} = \Pi_h \left( \sum_{i \in \mathcal{V}} \varphi_H^i \underline{\varepsilon}^i \right) \quad (18) $$

where $\Pi_h$ is the projector on the space $\mathrm{KA}_h^0(\omega_i)$. Discontinuities between star-patches vanish thanks to the multiplication by $\varphi_H^i$. The projector eases the computation of $R_H(\underline{w})$.

Note that in the same article, an enhanced estimate is proposed to better the lower bound:

$$ \|\underline{u}_{ex} - \underline{u}_H\|_{\Omega}^2 \geq \frac{R_H(\underline{w})^2}{\|\underline{w}\|_{\Omega}^2 - \|\underline{\varepsilon}_G\|_{\Omega}^2} \quad (19) $$

where $\underline{\varepsilon}_G$ is the solution of the following global problem:

$$ 
\begin{aligned}
& \text{Find } \underline{\varepsilon}_G \in \mathrm{KA}_H^0(\Omega) \text{ such that } \forall \underline{v} \in \mathrm{KA}_h^0(\omega_i) \\
& a(\underline{\varepsilon}_G, \underline{v}) = -a(\underline{w}, \underline{v})
\end{aligned} 
\quad (20) $$

In [9], the statically admissible stress field is built from a displacement field which is the sum of the solutions of local problems on star-patches with homogeneous boundary Dirichlet conditions:

$$ 
\begin{aligned}
& \text{Find } \underline{w}^i \in \mathrm{KA}^{0, \omega_i}(\omega_i) \text{ such that } \forall \underline{v} \in \mathrm{KA}^{0, \omega_i}(\omega_i) \\
& \int_{\omega_i} \underline{\varepsilon}(\underline{w}^i) : \mathbb{H} : \underline{\varepsilon}(\underline{v}) d\Omega = \int_{\omega_i} (\varphi_i \underline{f} - \underline{\sigma}_{\underline{H}} \mathrm{grad}(\varphi_i)) \underline{v} d\Omega - \int_{\omega_i} \varphi_i \underline{\sigma}_{\underline{H}} : \underline{\varepsilon}(\underline{v}) d\Omega
\end{aligned} 
\quad (21) $$

where $\mathrm{KA}^{0, \omega_i}$ is the space of continuous displacement fields that equal to zero on the boundary of the star-patch $\omega_i$. The continuous field $\underline{w} \in \mathrm{KA}^0(\Omega)$ is the sum of the solutions of the previous problem:

$$ \underline{w} = \sum_{i \in \mathcal{V}} \underline{w}^i \quad (22) $$

To conclude this brief review, there exist various techniques to construct $\underline{w}$. They always take advantage of the field computed during the estimation of an upper bound so that the extra-cost is very limited.

### 2.1.3 Substructured formulation

Let us consider a decomposition of domain $\Omega$ in $N_{sd}$ regular open subsets $(\Omega^{(s)})_s$ such that $\Omega^{(s)} \cap \Omega^{(s')} = \emptyset$ for $s \neq s'$ and $\bar{\Omega} = \bigcup_s \bar{\Omega}^{(s)}$. We note $\partial_g \Omega^{(s)} = \partial \Omega^{(s)} \cap \partial_g \Omega$ the Neumann border of subdomains.
The mechanical problem on the substructured configuration writes :

$$ \forall s \left\{ \begin{array}{l} \underline{u}^{(s)} \in \mathrm{KA}(\Omega^{(s)}) \\ \underline{\sigma}^{(s)} \in \mathrm{SA}(\Omega^{(s)}) \\ e_{CR_{\Omega^{(s)}}}(\underline{u}^{(s)}, \underline{\sigma}^{(s)}) = 0 \end{array} \right. \qquad (23) $$

and

$$ \forall (s, s') \text{ such that } \Omega^{(s)} \text{ and } \Omega^{(s')} \text{ are adjacent} \quad \left\{ \begin{array}{l} \mathrm{tr}(\underline{u}^{(s)}) = \mathrm{tr}(\underline{u}^{(s')}) \text{ on } \Gamma^{(s,s')} \\ \underline{\sigma}^{(s)} \cdot \underline{n}^{(s)} + \underline{\sigma}^{(s')}, \underline{n}^{(s')} = \underline{0} \text{ on } \Gamma^{(s,s')} \end{array} \right. \qquad (24) $$

The set of fields $\underline{u}$ defined on $\Omega$ such that $\underline{u}_{|\Omega^{(s)}} \in \mathrm{KA}(\Omega^{(s)})$ without interface continuity is a broken space which we note $\mathrm{KA}(\bigcup \Omega^{(s)})$.

### 2.1.4 Finite element approximation for the substructured problem

We assume that the mesh of $\bar{\Omega}$ and the substructuring are conforming. This hypothesis implies that each element only belongs to one subdomain and nodes are matching on the interfaces. For each subdomain, let us denote by the subscript b the degrees of freedom on the boundary of the subdomain and by the subscript i the degrees of freedom inside the subdomain.

Let $\mathbf{t}^{(s)}$ be the discrete trace operator on the interface. $\mathbf{t}^{(s)}$ enables to cast degrees of freedom from a complete subdomain to its interface. Using an adapted ordering, we have

$$ \mathbf{t}^{(s)} = (\mathbf{I}_{bb} \quad \mathbf{0}_{bi}) \qquad (25) $$

Therefore $\mathbf{t}^{(s)T}$ is the extension by zero operator: it extends data supported by the boundary to the whole subdomain. Let us introduce the unknown nodal reaction on the interface $\lambda^{(s)}$, the equilibrium of each subdomain writes:

$$ \mathbf{K}^{(s)} \mathbf{u}^{(s)} = \mathbf{f}^{(s)} + \mathbf{t}^{(s)T} \lambda^{(s)} \qquad (26) $$

Let $(\mathbf{A}^{(s)})$ and $(\mathbf{B}^{(s)})$ be the primal and dual assembly operator. Those operators are signed boolean operators. The dual operator $(\mathbf{B}^{(s)})$ enables to express the continuity of displacements and the primal operator $(\mathbf{A}^{(s)})$ enables to express the mechanical equilibrium of interface. Their number of columns is equal to the number of boundary degrees of freedom. $\mathbf{A}^{(s)}$ injects the boundary degrees of freedom of $\Omega^{(s)}$ in the global interface. Thus the number of rows of $\mathbf{A}^{(s)}$ is equal to the number of degrees of freedom on the global interface. The number of rows of $\mathbf{B}^{(s)}$ is equal to the number of connections between pairs of neighboring degrees of freedom.

In the case of two subdomains, we have $\sum_s \mathbf{B}^{(s)} \mathbf{t}^{(s)} \mathbf{u}^{(s)} = \mathbf{t}^{(1)} \mathbf{u}^{(1)} - \mathbf{t}^{(2)} \mathbf{u}^{(2)}$ and $\sum_s \mathbf{A}^{(s)} \lambda^{(s)} = \lambda^{(1)} + \lambda^{(2)}$. For more details on the assembly operators, the reader can refer to [11].

The discrete counterpart of the interface admissibility equations is:

$$ \left\{ \begin{array}{l} \sum_s \mathbf{B}^{(s)} \mathbf{t}^{(s)} \mathbf{u}^{(s)} = 0 \\ \sum_s \mathbf{A}^{(s)} \lambda^{(s)} = 0 \end{array} \right. \qquad (27) $$

Equations (26) and (27) form the discrete substructured system, which is equivalent to the global problem (11).

### 2.1.5 Domain decomposition solvers and admissible fields

Domain decomposition solvers are well described in many papers (see for instance [11] and the associated bibliography). The principle is to condense the global problem on the interface to create a smaller problem. In classical algorithms such as BDD [16] and FETI [7], this new interface problem is solved iteratively thanks to a projected preconditioned conjugate gradient. Each iteration implies two parallel solvings on subdomains with two dual operators (one for the preconditioning step, one for the direct step) so that local problems with Neumann boundary conditions and local problems with Dirichlet boundary conditions are alternatively solved. In [22] it was proved that the following fields could be processed at no extra cost:

* $(\underline{u}_D^{(s)})_s \in \mathrm{KA}(\Omega)$: displacement field which results from a Dirichlet problem and which is thus globally admissible

* $(\lambda_N^{(s)})_s$: nodal reactions which are balanced at the interface.

* $(\underline{u}_N^{(s)})_s \in \mathrm{KA}(\bigcup \Omega^{(s)})$: displacement field which results from a Neumann problem and which is not globally admissible
• $\underline{\sigma}_N^{(s)}$: the stress field associated to $\underline{u}_N^{(s)}$ ($\underline{\sigma}_N^{(s)} = \mathbb{H} : \underline{\varepsilon}(\underline{u}_N^{(s)})$). It can be used (with additional input $\lambda_N^{(s)}$) to build in parallel stress fields $\hat{\underline{\sigma}}_N^{(s)}$ which are statically admissible $\hat{\underline{\sigma}}_N^{(s)} = (\hat{\underline{\sigma}}_N^{(s)})_s \in \text{SA}(\Omega)$ using dedicated methods such as [14, 20, 23, 26].

In [28] we proved the following result where $\alpha$ (denoted $\sqrt{\mathbf{r}^T\mathbf{z}}$ in [28]) is the preconditioner-norm of the residual, a quantity that is actually computed by the solver:

$$ \alpha := \|\underline{u}_N - \underline{u}_D\|_{\mathbb{H},\Omega} = \sqrt{\mathbf{r}^T\mathbf{z}} \quad (28) $$

## 2.2 A posteriori upper bound of the error in substructured context

In [22], a first parallel error estimator in substructured context was introduced. It is based on the error in constitutive relation and reads:

$$ \|\underline{u}_{ex} - \underline{u}_{D}\|_{\Omega} = \sqrt{\sum_{s} \|\underline{u}_{ex}^{(s)} - \underline{u}_{D}^{(s)}\|_{\Omega^{(s)}}^{2}} \le \sqrt{\sum_{s} e_{CR_{\Omega^{(s)}}}^{2}(\underline{u}_{D}^{(s)}, \hat{\underline{\sigma}}_{N}^{(s)})} \quad (29) $$

In [28], we showed that this estimator mixes two different sources of error: the discretization error which is inherent to the use of the finite element method and the algebraic error which is due to the use of an iterative solver which would not exist for a direct solver. The algebraic error monitors the convergence of the solver and can be made as small as wished. Therefore, a second parallel error estimator was proposed in [28]:

$$ \|\underline{u}_{ex} - \underline{u}_{N}\|_{\Omega} \le \alpha + \sqrt{\sum_{s} e_{CR_{\Omega^{(s)}}}^{2}(\underline{u}_{N}^{(s)}, \hat{\underline{\sigma}}_{N}^{(s)})} \quad (30) $$

$$ \|\underline{u}_{ex} - \underline{u}_{D}\|_{\Omega} \le \alpha + \sqrt{\sum_{s} e_{CR_{\Omega^{(s)}}}^{2}(\underline{u}_{N}^{(s)}, \hat{\underline{\sigma}}_{N}^{(s)})} \quad (31) $$

This estimator separates the two sources of error. When the solver has converged the two displacements fields $\underline{u}_N$ and $\underline{u}_D$ are identical and equal to $\underline{u}_H$ (the algebraic error $\alpha$ is very close to zero). $\sqrt{\sum_s e_{CR_{\Omega^{(s)}}}^2(\underline{u}_N^{(s)}, \hat{\underline{\sigma}}_N^{(s)})}$ is the estimation of the discretization error. As a consequence, at convergence, the estimators (29), (30) and (31) are identical.

# 3 Lower bound of the error in substructured context

In this section, we extend sequential results to demonstrate guaranteed lower bounds of the error in substructured context. Moreover, we prove a theorem that enables the separation of sources of error in the lower bound. We also develop the methodology to build a continuous error estimate from parallel error estimation procedure. Finally, we extend those results to goal-oriented error estimation.

## 3.1 A first lower bound of the error

**Theorem 1.** Let $\underline{u}_{ex} \in \text{KA}(\Omega)$ be the exact solution, $(\underline{u}_D^{(s)})_s \in \text{KA}(\Omega)$ the displacement field defined in 2.1.5 and $\underline{w} \in \text{KA}^0(\Omega) \setminus \{\underline{0}\}$ then

$$ \|\underline{u}_{ex} - \underline{u}_{D}\|_{\Omega} \ge \frac{|R_{D}(\underline{w})|}{\sqrt{\sum_{s} \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^{2}}} \quad (32) $$

with

$$ R_D(\underline{w}) = \sum_{s} R_D^{(s)}(\underline{w}^{(s)}) \quad (33) $$

$$ R_D^{(s)}(\underline{w}^{(s)}) := \int_{\Omega^{(s)}} \underline{f} \cdot \underline{w}^{(s)} d\Omega + \int_{\partial_g \Omega^{(s)}} \underline{g} \cdot \underline{w}^{(s)} dS - \int_{\Omega^{(s)}} \underline{\varepsilon}(\underline{u}_D^{(s)}) : \mathbb{H} : \underline{\varepsilon}(\underline{w}^{(s)}) d\Omega^{(s)} \quad (34) $$

*Proof.* This property is the direct application of (16) where we replace the displacement field $\underline{u}_H \in \text{KA}(\Omega)$ by $(\underline{u}_D^{(s)})_s \in \text{KA}(\Omega)$. The residual $R_D(\underline{w})$ can be rewritten:

$$ \begin{align} R_D(\underline{w}) &= \sum_{s} R_D^{(s)}(\underline{w}^{(s)}) \nonumber \\ &= \sum_{s} \left( \int_{\Omega^{(s)}} \underline{f} \cdot \underline{w}^{(s)} d\Omega + \int_{\partial_g \Omega^{(s)}} \underline{g} \cdot \underline{w}^{(s)} dS - \int_{\Omega^{(s)}} \underline{\varepsilon}(\underline{u}_D^{(s)}) : \mathbb{H} : \underline{\varepsilon}(\underline{w}^{(s)}) d\Omega^{(s)} \right) \tag{35} \\ &= L(\underline{w}) - a(\underline{u}_D, \underline{w}) \nonumber \end{align} $$
□

In practice, we choose $\underline{w}^{(s)} \in \text{KA}^{00}(\Omega^{(s)}) \setminus \{\underline{0}\} \subset \text{KA}^{0}(\Omega) \setminus \{\underline{0}\}$, which corresponds to imposing the nullity along the interface and which is inexpensive since it does not imply exchanges between subdomains. In subsection 3.3, we will give details about the computation of $\underline{w}$. Therefore the computation of a lower bound is as parallel as for the upper bound. In the assessments in section 4, we will verify that this lower bound is as accurate as the one obtained in the sequential context and that the quality neither depends on the approach (primal or dual) nor on the substructuring. Moreover, this lower bound is computable whatever the state of the iterative solver (converged or not).

## 3.2 Lower bound with separation of sources of error

Continuing the philosophy of separating the sources of error as detailed in [28], we propose a second lower bound :

**Theorem 2.** Let $\underline{u}_{ex} \in \text{KA}(\Omega)$ be the exact solution, $(\underline{u}_D^{(s)})_s \in \text{KA}(\Omega)$ and $(\underline{u}_N^{(s)})_s \in \text{KA}(\bigcup \Omega^{(s)})$ the displacement fields defined in 2.1.5 and $\underline{w} \in \text{KA}^0(\Omega) \setminus \{\underline{0}\}$, then

$$ \|\underline{u}_{ex} - \underline{u}_D\|_{\Omega} \ge \left| \frac{|R_N(\underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} - \frac{|a(\underline{u}_D - \underline{u}_N, \underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} \right| \quad (36) $$

which leads to the coarser bound

$$ \|\underline{u}_{ex} - \underline{u}_D\|_{\Omega} \ge \frac{|R_N(\underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} - \alpha \quad (37) $$

with

$$ R_N(\underline{w}) = \sum_s R_N^{(s)}(\underline{w}^{(s)}) \quad (38) $$

$$ R_N^{(s)}(\underline{w}^{(s)}) := \int_{\Omega^{(s)}} \underline{f} \cdot \underline{w}^{(s)} d\Omega + \int_{\partial g \cap \Omega^{(s)}} \underline{g} \cdot \underline{w}^{(s)} dS - \int_{\Omega^{(s)}} \underline{\varepsilon} \cdot (\underline{u}_N^{(s)}) : \mathbb{H} : \underline{\varepsilon} \cdot (\underline{w}^{(s)}) d\Omega^{(s)} \quad (39) $$

*Proof.* The proof of the first inequality is based on theorem 1 and on the triangle inequality :

$$ \begin{aligned} \|\underline{u}_{ex} - \underline{u}_D\|_{\Omega} &\ge \frac{|R_D(\underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} \\ \|\underline{u}_{ex} - \underline{u}_D\|_{\Omega} &\ge \frac{|L(\underline{w}) - a(\underline{u}_D, \underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} \\ &\ge \frac{|L(\underline{w}) - a(\underline{u}_N, \underline{w}) - a(\underline{u}_D - \underline{u}_N, \underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} \\ &\ge \frac{|R_N(\underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} - \frac{|a(\underline{u}_D - \underline{u}_N, \underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} \end{aligned} \quad (40) $$

which proves (36). (37) is simply based on the remark that

$$ \left| \frac{|R_N(\underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} - \frac{|a(\underline{u}_D - \underline{u}_N, \underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} \right| \ge \frac{|R_N(\underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} - \frac{|a(\underline{u}_D - \underline{u}_N, \underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} \quad (41) $$

and using twice the Cauchy-Schwarz inequality, we have :

$$ \begin{aligned} |a(\underline{u}_D - \underline{u}_N, \underline{w})| &= \left| \sum_s \int_{\Omega^{(s)}} \underline{\varepsilon} \cdot (\underline{u}_D^{(s)} - \underline{u}_N^{(s)}) : \mathbb{H} : \underline{\varepsilon} \cdot (\underline{w}^{(s)}) \right| \\ &\le \sum_s \|\underline{u}_N^{(s)} - \underline{u}_D^{(s)}\|_{\Omega^{(s)}} \|\underline{w}^{(s)}\|_{\Omega^{(s)}} \quad (42) \\ &\le \sqrt{\sum_s \|\underline{u}_N^{(s)} - \underline{u}_D^{(s)}\|_{\Omega^{(s)}}^2} \sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2} \end{aligned} $$

Finally, using the equality (28):

$$ \|\underline{u}_{ex} - \underline{u}_D\|_{\Omega} \ge \frac{|R_N(\underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}} - \alpha \quad (43) $$

□
As said earlier, the term $\alpha$ is a measure of the residual so it is purely algebraic whereas the first term of the inequality is mainly driven by the discretization error. During the first iterations, the second lower bound is not accurate because the algebraic error prevails so that $\frac{|L(w) - a(w_N, w)|}{\sqrt{\sum_s ||w^{(s)}||_{\Omega^{(s)}}^2}} - \alpha$ is negative and it is a trivial lower bound of the positive true error $||\underline{u}_{ex} - \underline{u}_D||_{\Omega}$. When the solver reaches convergence, the three lower bounds in theorems 1 and 2 are identical.

### 3.3 Reconstruction of admissible field $\underline{w}$

The upper bounds of the error (30) or (31) require the construction of a statically admissible field $(\hat{\underline{\sigma}}_N^{(s)})_s$. In case the flux-free technique is chosen, it is possible to construct a continuous field $\underline{w}^{(s)} \in \text{KA}^{00}(\Omega^{(s)})$ using the methodology developed in [20] for each subdomain in parallel. As a consequence,

$$ \left( \left( \prod_h \left( \sum_{i \in \mathcal{V}} \phi_H^i \underline{\varepsilon}^i \right) \right)^{(s)} \right)_s \in \text{KA}(\bigcup \Omega^{(s)}) \quad (44) $$

In order to have $\underline{w}^{(s)} \in \text{KA}^{00}(\Omega^{(s)})$, we choose not to sum the contributions from the nodes located on the interface of the subdomain $(s)$. They will be denoted as $\mathcal{V}_\Gamma$. $\underline{w}^{(s)}$ is defined by:

$$ \underline{w}^{(s)} = \left( \prod_h \left( \sum_{i \in \mathcal{V} \setminus \mathcal{V}_\Gamma} \phi_H^i \underline{\varepsilon}^i \right) \right)^{(s)} \quad (45) $$

Therefore $\underline{w}^{(s)} \in \text{KA}^{00}(\Omega^{(s)})$ and $\underline{w} = (\underline{w}^{(s)})_s \in \text{KA}^0(\Omega)$.

*Remark.* Using the subtle trick in [20], the computation of the discretized fields is eased. Indeed:

$$ \mathbf{w}^{(s)} = \sum_{i \in \mathcal{V} \setminus \mathcal{V}_\Gamma} \varphi_H^i \odot \mathbf{e}^{(s)} \quad (46) $$

and

$$ \mathbf{R}_D^{(s)}(\mathbf{w}^{(s)}) = \sum_{i \in \mathcal{V} - \mathcal{V}_\Gamma} \varphi_H^i \odot \mathbf{R}_D^{(s)}(\mathbf{e}^{(s)}) \quad (47) $$

where $\varphi_H^i$ gathers the nodal values of the shape function $\varphi_H^i$ projected on the richer space used to solve the star-patch problem whose $\mathbf{e}^{(s)}$ is the discretized solution and where $\odot$ represents the term by term multiplication.

*Remark.* If the method chosen to construct the admissible field is based on elements problems [4], it is always possible to construct a displacement field $\underline{w}^{(s)} \in \text{KA}^{00}(\Omega^{(s)})$ by computing the mean value along the edges inside the subdomains and imposing zero along the interfaces between subdomains.

### 3.4 Goal-oriented error estimation

Goal-oriented error estimation offers the possibility to have upper and lower bounds of the unknown exact value of a quantity of interest. Among various techniques, extractors (see [1] for instance) are the most common tools to define linear quantities of interest. They lead to the definition and the solving of an adjoint problem.

#### 3.4.1 Definition of the linear quantity of interest and of the adjoint problem

Let $\tilde{L}$ be the linear functional defining the quantity of interest $I$:

$$ I = \tilde{L}(\underline{u}) = \int_{\Omega} (\underline{\sigma}_{\Sigma} : \underline{\varepsilon}(\underline{u}) + f_{\Sigma} \underline{u}) d\Omega \quad (48) $$

where $\underline{\sigma}_{\Sigma}$ and $f_{\Sigma}$ are extractors.

We introduce the affine subspace of statically admissible fields (adjoint SA-fields) for the adjoint problem:

$$ \widetilde{\text{SA}}(\omega) = \left\{ \underline{\underline{\tau}} \in (L^2(\omega))_{\text{sym}}^{d \times d}; \forall \underline{u} \in \text{KA}^{00}(\omega), \int_{\omega} \underline{\underline{\tau}} : \underline{\underline{\varepsilon}}(\underline{u}) d\omega = \tilde{L}(\underline{u}) \right\} \quad (49) $$

The adjoint problem set on $\Omega$ can be formulated as:

$$ \text{Find } (\tilde{\underline{u}}_{ex}, \tilde{\underline{\sigma}}_{ex}) \in \text{KA}^0(\Omega) \times \widetilde{\text{SA}}(\Omega) \text{ such that } e_{CR_{\Omega}}(\tilde{\underline{u}}_{ex}, \tilde{\underline{\sigma}}_{ex}) = 0 \quad (50) $$

The solution to this problem, named exact solution, exists and is unique.
Remark. The formulation (50) is equivalent to the classical following formulation:

$$ \text{Find } \tilde{\underline{u}}_{ex} \in \text{KA}^0(\Omega) \text{ such that } \forall \underline{v} \in \text{KA}^0(\Omega), a(\tilde{\underline{u}}_{ex}, \underline{v}) = \tilde{L}(\underline{v}) \quad (51) $$

where a is the classical bilinear form (8).

The adjoint problem is usually solved with the finite element method. The mesh can differ from the one used for the forward problem. The approximated adjoint displacement is $\tilde{\underline{u}}_{\tilde{H}}$.

The discretization error for the adjoint problem is

$$ e_{discr} = \|\tilde{\underline{\epsilon}}_{discr}\|_{\Omega} = \|\tilde{\underline{u}}_{ex} - \tilde{\underline{u}}_{\tilde{H}}\|_{\Omega} \quad (52) $$

### 3.4.2 Error estimation of quantities of interest

As said earlier the adjoint problem is meant to extract one quantity of interest in the forward problem. Let $I_{ex} = \tilde{L}(\underline{u}_{ex})$ be the unknown exact value of the quantity of interest. $I_H = \tilde{L}(\underline{u}_D)$ is an approximation of this quantity of interest.

A bounding of the exact value of the quantity of interest $I_{ex}$ is [12, 13]:

$$ |I_{ex} - I_H - I_{HH2}| \le \frac{1}{2} e_{CR_{\Omega}}(\underline{u}_H, \underline{\hat{g}}_H) e_{CR_{\Omega}}(\tilde{\underline{u}}_{\tilde{H}}, \tilde{\underline{\hat{g}}}_{\tilde{H}}) \quad (53) $$

where

$$ I_{HH2} = \int_{\Omega} \frac{1}{2} (\underline{\hat{g}}_{\tilde{H}} + \mathbb{H} : \underline{\varepsilon}(\tilde{\underline{u}}_{\tilde{H}})) : \mathbb{H}^{-1} : (\underline{\hat{g}}_{\tilde{H}} - \mathbb{H} : \underline{\varepsilon}(\underline{u}_H)) d\Omega \quad (54) $$

and where $\underline{\hat{g}}_{\tilde{H}} \in \tilde{\mathbb{S}}\tilde{\mathbb{A}}_H(\Omega)$.

The error on the quantity of interest can also be estimated using the parallelogram identity [24]:

$$ \left\{ \begin{aligned} I_{ex} - I_H = \tilde{L}(\underline{\epsilon}_{discr}) &= a(\underline{\epsilon}_{discr}, \tilde{\underline{\epsilon}}_{discr}) \\ &= a(\kappa \underline{\epsilon}_{discr}, \frac{1}{\kappa} \tilde{\underline{\epsilon}}_{discr}) \\ &= \frac{1}{4} [\|\kappa \underline{\epsilon}_{discr} + \frac{1}{\kappa} \tilde{\underline{\epsilon}}_{discr}\|_{\Omega}^2 - \|\kappa \underline{\epsilon}_{discr} - \frac{1}{\kappa} \tilde{\underline{\epsilon}}_{discr}\|_{\Omega}^2] \end{aligned} \right. \quad (55) $$

where $\kappa$ is a scalar parameter whose optimal value is:

$$ \kappa = \frac{e_{CR_{\Omega}}(\tilde{\underline{u}}_H, \underline{\hat{g}}_H)}{e_{CR_{\Omega}}(\underline{u}_H, \underline{\hat{g}}_H)} \quad (56) $$

This optimal value minimizes the difference and thus improves the quality of the bounding. Introducing upper and lower bounds:

$$ \begin{aligned} \beta_{inf}^{+} &\le \|\kappa \underline{\epsilon}_{discr} + \frac{1}{\kappa} \tilde{\underline{\epsilon}}_{discr}\|_{\Omega}^2 \le \beta_{sup}^{+} \\ \beta_{inf}^{-} &\le \|\kappa \underline{\epsilon}_{discr} - \frac{1}{\kappa} \tilde{\underline{\epsilon}}_{discr}\|_{\Omega}^2 \le \beta_{sup}^{-} \end{aligned} \quad (57) $$

it is possible to obtain lower and upper bounds of the error on the quantity of interest:

$$ \frac{1}{4} \beta_{inf}^{+} - \frac{1}{4} \beta_{sup}^{-} \le I_{ex} - I_H \le \frac{1}{4} \beta_{sup}^{+} - \frac{1}{4} \beta_{inf}^{-} \quad (58) $$

### 3.4.3 Application to the substructured context

Let us suppose that the forward and adjoint problems are solved on the same mesh and on the same substructuring. Since the two problems share the same stiffness matrices on every subdomain, they can be solved simultaneously using a block algorithm. Following the same methodology described in section 2.1.5, one can compute the following admissible fields for the adjoint problem:

* $(\tilde{\underline{u}}_D^{(s)})_s \in \text{KA}^0(\Omega)$: displacement field which results from a Dirichlet problem and which is thus globally admissible

* $\tilde{\lambda}_N^{(s)}$: nodal reactions which are balanced at the interface.

* $(\tilde{\underline{u}}_N^{(s)})_s \in \text{KA}(\bigcup \Omega^{(s)})$: displacement field which results from a Neumann problem and which is not globally admissible

* $\tilde{\underline{\lambda}}_N^{(s)}$: stress field associated to $\tilde{\underline{u}}_N^{(s)}$ ($\tilde{\underline{\lambda}}_N^{(s)} = \mathbb{H} : \underline{\varepsilon}(\tilde{\underline{u}}_N^{(s)})$). It can be used (with additional input $\tilde{\lambda}_N^{(s)}$) to build in parallel stress fields $\tilde{\underline{\hat{g}}}_N^{(s)}$ which are statically admissible $\tilde{\underline{\hat{g}}}_N^{(s)} = (\tilde{\underline{\hat{g}}}_N^{(s)})_s \in \tilde{\mathbb{S}}\tilde{\mathbb{A}}(\Omega)$.
We also have the equality that expresses the distance between the Neumann and Dirichlet displacement fields in terms of the algebraic residual $\tilde{\alpha}$ of the adjoint problem :

$$ \| \tilde{u}_N - \tilde{u}_D \|_{H, \Omega} = \tilde{\alpha} \qquad (59) $$

For more details about the computation of those fields, the reader can refer to [27].

The upper bounds of the global error presented in section 2.2 and the lower bound in the theorem 1 can be applied on the adjoint problem.

**Upper and lower bounds for goal-oriented error estimation** We demonstrate two properties that give upper and lower bounds of the terms in the parallelogram identity. The properties are merely the application of the results on goal-oriented error estimation [12, 13, 20] into a substructured context.

**Corollary 1.** *Using notations of paragraph 2.1.5 and 3.4.3*

$$ \begin{aligned} \beta_{inf}^{+} &\le \|\kappa \underline{e}_{discr} + \frac{1}{\kappa} \tilde{\underline{e}}_{discr}\|_{\Omega}^{2} \le \beta_{sup}^{+} \\ \beta_{inf}^{-} &\le \|\kappa \underline{e}_{discr} - \frac{1}{\kappa} \tilde{\underline{e}}_{discr}\|_{\Omega}^{2} \le \beta_{sup}^{-} \end{aligned} \qquad (60) $$

with

$$ \left\{ \begin{aligned} \beta_{sup}^{+} &= 2\sqrt{\sum_{s} e_{CR_{\Omega(s)}}(\underline{u}_{D}^{(s)}, \underline{\hat{\sigma}}_{N}^{(s)})^{2}} \sqrt{\sum_{s} e_{CR_{\Omega(s)}}(\tilde{\underline{u}}_{D}^{(s)}, \tilde{\underline{\hat{\sigma}}}_{N}^{(s)})^{2}} \\ &\quad + 2\sum_{s} \int_{\Omega^{(s)}} (\underline{\hat{\sigma}}_{N}^{(s)} - \mathbb{H}\underline{\varepsilon}(\tilde{\underline{u}}_{D}^{(s)})) : \mathbb{H}^{-1} : (\underline{\hat{\sigma}}_{N}^{(s)} - \mathbb{H} : \underline{\varepsilon}(\underline{u}_{D}^{(s)})) d\Omega \\ \beta_{inf}^{+} &= \frac{(\kappa R_{D}(\underline{z}^{+}) + \frac{1}{\kappa} \tilde{R}_{D}(\underline{z}^{+}))^{2}}{\sum_{s} \|(\underline{z}^{+})^{(s)}\|_{\Omega^{(s)}}^{2}} \\ &\quad \text{with } (\underline{z}^{+})^{(s)} = \kappa \underline{w}^{(s)} + \frac{1}{\kappa} \tilde{\underline{w}}^{(s)} \end{aligned} \right. \qquad (61) $$

and

$$ \left\{ \begin{aligned} \beta_{sup}^{-} &= 2\sqrt{\sum_{s} e_{CR_{\Omega(s)}}(\underline{u}_{D}^{(s)}, \underline{\hat{\sigma}}_{N}^{(s)})^{2}} \sqrt{\sum_{s} e_{CR_{\Omega(s)}}(\tilde{\underline{u}}_{D}^{(s)}, \tilde{\underline{\hat{\sigma}}}_{N}^{(s)})^{2}} \\ &\quad - 2\sum_{s} \int_{\Omega^{(s)}} (\underline{\hat{\sigma}}_{N}^{(s)} - \mathbb{H}\underline{\varepsilon}(\tilde{\underline{u}}_{D}^{(s)})) : \mathbb{H}^{-1} : (\underline{\hat{\sigma}}_{N}^{(s)} - \mathbb{H} : \underline{\varepsilon}(\underline{u}_{D}^{(s)})) d\Omega \\ \beta_{inf}^{-} &= \frac{(\kappa R_{D}(\underline{z}^{-}) - \frac{1}{\kappa} \tilde{R}_{D}(\underline{z}^{-}))^{2}}{\sum_{s} \|(\underline{z}^{-})^{(s)}\|_{\Omega^{(s)}}^{2}} \\ &\quad \text{with } (\underline{z}^{-})^{(s)} = \kappa \underline{w}^{-}(s) - \frac{1}{\kappa} \tilde{\underline{w}}^{-}(s) \end{aligned} \right. \qquad (62) $$

One has to pay attention to the computation of the parameter $\kappa$ which implies exchanges between subdomains. Indeed, this coefficient is defined by:

$$ \kappa = \frac{\sqrt{\sum_{s} e_{CR_{\Omega(s)}}(\tilde{\underline{u}}_{D}^{(s)}, \tilde{\underline{\hat{\sigma}}}_{N}^{(s)})}}{\sqrt{\sum_{s} e_{CR_{\Omega(s)}}(\underline{u}_{D}^{(s)}, \underline{\hat{\sigma}}_{N}^{(s)})}} \qquad (63) $$

Anyhow this exchange between subdomains is already done at the end of the parallel error estimation to obtain global measures.

Despite the possibility to separate contributions in terms $\sqrt{\sum_{s} e_{CR_{\Omega(s)}}(\tilde{\underline{u}}_{D}^{(s)}, \tilde{\underline{\hat{\sigma}}}_{N}^{(s)})}$ and $\sqrt{\sum_{s} e_{CR_{\Omega(s)}}(\underline{u}_{D}^{(s)}, \underline{\hat{\sigma}}_{N}^{(s)})}$, the full separation in the lower bounds $\beta_{inf}^{-}$ and $\beta_{inf}^{-}$ is a complex task since the parameter $\kappa$ is the ratio of errors mixing algebraic and discretization sources.

However, the separation of sources for both global errors enables steering the iterative solver by an objective of precision of the quantity of interest (see [27]). The computation of $\beta_{inf}^{-}$ and $\beta_{inf}^{-}$ after convergence improves the bounding.

# 4 Numerical assessment

For all numerical examples, the behavior is linear, isotropic and elastic. The Young modulus is 1 Pa and the Poisson coefficient is 0.3.
## 4.1 Structure with exact solution

Let us consider a square linear elastic structure $\Omega = [-3l; 3l] \times [-3l; 3l]$ with homogeneous Dirichlet boundary conditions and plane strain hypothesis. The domain is subjected to a polynomial body force such that the exact solution is known:

$$ \underline{u}_{ex} = (x + 3l)(x - 3l)(y + 3l)(y - 3l) ((y - 3l)^2 \underline{e}_x + (y + 3l)\underline{e}_y) $$

The mesh is made out of first order Lagrange triangles. As shown in Figure 1 the structure is decomposed into 9 regular subdomains.

<figure><img src="image_2.png" /><figcaption>Figure 1: Substructuring</figcaption></figure>

We solve the problem with a BDD solver (primal approach). We use the Flux-free technique [20] to build statically admissible stress fields; each star-patch problem is solved by subdividing each element into 12 elements (h-refinement technique). For the sake of simplicity, we note:

* $\rho = \frac{|L(\underline{w}) - a(\underline{u}_D, \underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}}$ the lower bound of the error

* $\rho_{discr} = \frac{|L(\underline{w}) - a(\underline{u}_N, \underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}}$ the discretization part of the lower bound

* $\rho_{alg} = \frac{|a(\underline{u}_D - \underline{u}_N, \underline{w})|}{\sqrt{\sum_s \|\underline{w}^{(s)}\|_{\Omega^{(s)}}^2}}$ the algebraic part of the lower bound

* $\rho_{bis} = \rho_{discr} - \alpha$ the lower bound with separation of sources of error

* $\theta = \sqrt{\sum_s e_{CR_{\Omega^{(s)}}}^2 (\underline{u}_D^{(s)}, \underline{\hat{u}}_N^{(s)})}$ the upper bound of the error

* $\theta_{discr} = \sqrt{\sum_s e_{CR_{\Omega^{(s)}}}^2 (\underline{u}_N^{(s)}, \underline{\hat{u}}_N^{(s)})}$ the discretization part of the upper bound

The quantities with the superscript $seq$ are computed with a sequential simulation (no substructuring and use of a direct solver).

### 4.1.1 Effects of the substructuring on the computation of the lower bound at convergence

In this subsection, we compare the lower bound obtained in a sequential simulation with the lower bound obtained in the substructured context when the solver has converged. Since the study is done at convergence, the primal and dual approaches are equivalent.

On figure 2, we observe that the bounds are the same for sequential and substructured computations. We also verify that the exact error is between upper and lower bounds. The convergence is the one expected for such a regular problem (h-slope).

Then, we compare the bounds for several substructuring as illustrated in figure 3. Table 1 gathers the lower bounds for sequential, primal and dual approaches computed at convergence normalized by the bounds for sequential computation. We observe that the substructuring has quasi no influence on the accuracy of the lower bound.

### 4.1.2 Separation of sources of error

In this subsection, we illustrate the separation of sources of error in the lower bound.

On the first graph in figure 4, we give the evolution of the upper and lower bounds and of the true error until the fifth iteration. We observe the fast convergence of those bounds. On the same graph, we also visualize the
<figure><img src="image_3.png" /><figcaption>Figure 2: Evolution of the upper and lower bounds of the error in function of the mesh size h</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Different substructures</figcaption></figure>

discretization parts $\theta_{discr}$ and $\rho_{discr}$. The bounding is more precise and enables to define after a few iterations the interval in which the true error is located at convergence. The second graph in figure 4 gives the evolution of the terms in theorem 2. The quantity $\alpha$ strictly decreases along the iterations, the evolution of the quantity $\rho_{alg}$ is not as smooth.

Figure 5 represents the two lower bounds of theorem 2. As expected, the second lower bound is not precise at the beginning since it gives a negative value. However, the zero of this bound enables to tell the moment when the algebraic error becomes smaller than the discretization error.

## 4.2 Pre-cracked structure

We now consider a pre-cracked structure decomposed into 16 subdomains as illustrated in Figure 6. The displacements at the base of the structure and on the larger hole are imposed to be zero. The upper-left part and the second hole are subjected to a constant unit pressure. We made the hypothesis of plane stress. The quantity of interest is the mean of the stress component $\sigma_{xx}$ on a region $\omega$ close to the crack. In Figure 6, the loading of the reference problem is in blue and the loading of adjoint problem is in orange. We used the FETI algorithm (dual approach) to solve the interface problem and the statically admissible stress fields are built using the Flux-free technique [3, 20] with h-refinement technique for the solving of local problems on star-patches (each element is divided into 16 elements).

### 4.2.1 Separation of sources of error in upper and lower bounds

For this paragraph, we consider only the forward problem. The mesh used for this computation is composed of 4370 degrees of freedom. We compute the global upper and lower bounds during the iterations and the discretization upper and lower bounds along the iterations on the same graph in figure 7.
<table><caption>Table 1: Relative upper and lower bounds of the error for various substructuring</caption><thead><tr><th></th><th>seq</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th></tr></thead><tbody><tr><td>θ<sub>seq</sub></td><td>1</td><td>1.0002</td><td>0.9993</td><td>0.9989</td><td>1.0004</td><td>0.9993</td></tr><tr><td>ρ<sub>seq</sub></td><td>1</td><td>0.9981</td><td>0.9964</td><td>0.9926</td><td>0.9958</td><td>0.9886</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: Separation of sources of error in upper and lower bounds</figcaption></figure>

As for the previous example, we observe the fast convergence of the upper and lower bounds. This graph also illustrates the separation of sources of error. Finally, we observe that the residual $\alpha$ decreases with the iterations and that, again, the evolution of $\rho_{alg}$ is not as regular.

In figure 8, we give the discretization upper and lower bounds, which define a discretization envelope, for all iterations on the first graph and until the seventh iteration on the second graph.

The discretization envelope can be used to define the stopping criterion of the solver. In [28], the proposed criterion was to stop when the algebraic error was ten times smaller than the discretization error. On this example, the solver would stop at the sixth iteration. A new stopping criterion could be to stop when the algebraic error is smaller than the discretization part of the lower bound, which would lead to stop to at the fourth iteration.

## 4.3 Goal-oriented error estimation

In this subsection, we propose an auto-adaptive strategy to steer the iterative solver by an objective of precision on a quantity of interest. In this example, the objective of precision will be five percent. We used the FETI algorithm (dual approach) to solve the interface problem and the statically admissible stress fields are built using the Flux-free technique [20] with h-refinement technique for the solving of local problems on star-patches (each element is divided into 4 elements). Since the forward and adjoint problems are auto-adjoint (due to the symmetry of the bilinear form), we solve the two problems simultaneously using a block conjugate gradient (see [27] for more details). The separation of sources of error in the lower bounds of global error on forward and adjoint problem enables the definition of the following criterion: **STOP when $\alpha < \rho_{discr}$ and $\tilde{\alpha} < \tilde{\rho}_{discr}$** which expresses the fact that the residual is out of the discretization envelope so that the algebraic error is negligible in comparison with the discretization error. Using equation (58), we have the following upper and lower bounds on the unknown exact value of the quantity of interest $I_{ex}$:

$$I_{ex}^{-} = I_{H} + \frac{1}{4}\beta_{inf}^{+} - \frac{1}{4}\beta_{sup}^{-} \le I_{ex} \le I_{H} + \frac{1}{4}\beta_{sup}^{+} - \frac{1}{4}\beta_{inf}^{-} = I_{ex}^{+} \quad (64)$$

We will compare the bounds $I_{ex}^{-}$ and $I_{ex}^{+}$ in case the quantities $\beta_{inf}^{+}$ and $\beta_{inf}^{-}$ are computed using expressions in Corollary 1 and in case they are chosen equal to zero ($\beta_{inf}^{+} = \beta_{inf}^{-} = 0$), which is equivalent to the bounding in equation (53) with $u_H = u_D$, $\tilde{u}_H = \tilde{u}_D$, $\hat{u}_H = \hat{u}_D$ and $\tilde{\hat{u}}_H = \tilde{\hat{u}}_D$.

We start with a first mesh which is a little bit refined near the quantity of interest in order to have several elements in the region $\omega$. The discretization error is computed at iteration 1 (which is more relevant than the initialization Iteration 0). The criterion is defined and the solver iterates until the criterion is reached. The error is estimated once again to verify that the discretization error has not changed too much. We give in table 2 the evolution of the residual for forward and adjoint problems and the bounds on the global errors for the two problems.

Regarding the quantity of interest, at the sixth iteration, we obtain the values presented in table 3.

In table 4, we give the upper and lower bounds of the unknown exact value $I_{ex}$ of the quantity of interest with and without the use of the lower bounds. We observe that the use of the lower bounds enables to reduce the width by 44 %. At the end of the first solving, the error on the quantity of interest is 22.224 %.
<figure><img src="image_6.png" /><figcaption>Figure 5: Evaluation of the two lower bounds during iterations</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Loading of forward (blue) and adjoint problems (orange), domain decomposition</figcaption></figure>

We also give, for both problems, the contribution from each subdomain thanks to the following quantities :

$$ \eta^{(s)} = \frac{e_{CR_{\Omega^{(s)}}}(\underline{u}_N^{(s)}, \underline{\hat{g}}_N^{(s)})}{\sqrt{\sum_s e_{CR_{\Omega^{(s)}}}^2(\underline{u}_N^{(s)}, \underline{\hat{g}}_N^{(s)})}} \quad \text{and} \quad \tilde{\eta}^{(s)} = \frac{e_{CR_{\Omega^{(s)}}}(\tilde{\underline{u}}_N^{(s)}, \tilde{\underline{\hat{g}}}_N^{(s)})}{\sqrt{\sum_s e_{CR_{\Omega^{(s)}}}^2(\tilde{\underline{u}}_N^{(s)}, \tilde{\underline{\hat{g}}}_N^{(s)})}} \quad (65) $$

and plot the contributions on the Figure 9.

As expected, for the adjoint problem, the error is mainly located in the sixth subdomain, which is the subdomain with the load. For the forward problem, the error is more diffuse.

In order to reach the objective of precision, we decide to first improve the quality of the solution of the forward problem. To do so, since the error is not located in few subdomains, we decide to refine the mesh on the whole structure. Of course, we could have used the error map provided by the error estimator and a remeshing criterion to process adaptive remeshing (see for instance [33, 5, 17, 2, 6]). For sake of simplicity, a refinement by splitting is performed. We also reuse the search directions computed during the first solving in order to speed up the next one. The 12 interface vectors (2 interface vectors -one for the forward problem, one for the adjoint problem- computed at each iteration of the first solving that converged in 6 iterations) corresponding to the search directions are projected on the new mesh and used as additional constraints thanks to augmented-Krylov methods [31]. Since the dual approach was chosen, we use a two-level FETI algorithm [8] to take into account the additional constraints. It does not modify the methodology to construct admissible fields nor the error estimator. Recycling search directions leading to a better initialization, we compute the error estimation at iteration 0 and define the stopping criterion. Once the criterion is reached, the error is estimated once again and if the criterion is checked, the solver is stopped. We give in the table 5 the evolution of the residuals and of the bounds on the global errors for the forward and adjoint problems. We can observe that the first residual is comparable to the last residual of the first solving.

Regarding the quantity of interest, we obtain the values presented in table 6.
<figure><img src="image_8.png" /><figcaption>Figure 7: Pre-cracked structure : Evolution of the upper and lower bounds (global, discretization part, algebraic part) and of the residual during the iterations</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Pre-cracked structure : Discretization envelope and residual against iterations</figcaption></figure>

In table 7, we give the upper and lower bounds of the unknown exact value $I_{ex}$ of the quantity of interest with and without the use of the lower bounds. We observe that the use of the lower bounds improves the bounding. At the end of the second solving, the uncertainty on the quantity of interest is 6.7893 %.

In order to reach the objective of precision, we decide to refine the mesh only for the sixth subdomain to improve the quality of the adjoint solution. Refining this subdomain's discretization introduces incompatibilities at the interface that can be easily managed thanks to transfer matrix as explained in [27]. This incompatibility does not affect the error estimator since the quantity of interest is located far from the subdomain's boundary. Once again, we reuse the search directions of the first two solutions to speed up the third solving. As the discretization of the interface is not modified (see [27]), the interface vectors can be directly used as additional constraints.

Regarding the quantity of interest, we obtain the values presented in table 9.

In table 10, we give the upper and lower bounds of the unknown exact value $I_{ex}$ of the quantity of interest with and without the use of the lower bounds. We observe that the use of the lower bounds enables to reduce the width by 44 %. At the end of the third solving, the objective of precision is reached and the error on the exact value of the quantity of interest is smaller than 4 %.

Finally, we give the global errors and residuals against cumulative iteration on Figure 10 and the evolution of the approximated value of the quantity of interest $I_H$ and the upper and lower bounds for $I_{ex}$ in Figure 11.

# 5 Conclusion

In this paper, we proposed a strict lower bound of the error in a substructured context which can be computed in parallel using the admissible fields built for the computation of the upper bound. Moreover, a theorem gives a second lower bound that separates the algebraic error from the discretization error. As illustrated on mechanical examples, the lower bounds are quasi independent from the substructuring and are as accurate as the sequential lower bound. The examples also show the separation of sources of error for the lower bound. Finally, we proposed an auto-adaptive strategy to steer the iterative solver by an objective of precision on a quantity of interest. Benefiting from the separation of sources of error to avoid oversolving and the recycling of Krylov subspaces, the strategy automatically defines a sequence of optimized solvings.
<table><caption>Table 2: Pre-cracked structure: First mesh</caption><thead><tr><th>Iteration</th><th>θ<sub>discr</sub></th><th>ρ<sub>discr</sub></th><th>θ̂<sub>discr</sub></th><th>ρ̃<sub>discr</sub></th><th>α</th><th>α̃</th></tr></thead><tbody><tr><td>0</td><td></td><td></td><td></td><td></td><td>260.07</td><td>0.26041</td></tr><tr><td>1</td><td>11.45</td><td>7.0008</td><td>0.12867</td><td>1.5457 10<sup>-3</sup></td><td>48.339</td><td>0.14931</td></tr><tr><td>2</td><td></td><td></td><td></td><td></td><td>35.991</td><td>8.5664 10<sup>-2</sup></td></tr><tr><td>3</td><td></td><td></td><td></td><td></td><td>16.649</td><td>4.8205 10<sup>-2</sup></td></tr><tr><td>4</td><td></td><td></td><td></td><td></td><td>5.5966</td><td>1.2462 10<sup>-2</sup></td></tr><tr><td>5</td><td></td><td></td><td></td><td></td><td>3.5765</td><td>8.4996 10<sup>-3</sup></td></tr><tr><td>6</td><td>9.9004</td><td>6.7105</td><td>0.12682</td><td>1.5893 10<sup>-3</sup></td><td>0.94805</td><td>1.5156 10<sup>-3</sup></td></tr></tbody></table>

<table><caption>Table 3: Pre-cracked structure: First mesh : error estimation</caption><thead><tr><th>I<sub>H</sub></th><th>β<sup>-</sup><sub>inf</sub></th><th>β<sup>+</sup><sub>inf</sub></th><th>1/4β<sup>+</sup><sub>sup</sub></th><th>1/4β<sup>-</sup><sub>sup</sub></th></tr></thead><tbody><tr><td>3.1505</td><td>1.0087</td><td>1.2129</td><td>0.68424</td><td>0.57132</td></tr></tbody></table>

## Acknoledgment

The author would like to thank Professor Pedro Diez (Universitat Politècnica de Catalunya) for the helpful and fruitful discussions.

## References

[1] R. Becker and R. Rannacher. A feed-back approach to error control in finite element methods: Basic analysis and examples. *Journal of Numerical Mathematics*, 4:237–264, 1996.

[2] Emmanuel Bellenger and Patrice Coorevits. Adaptive mesh refinement for the control of cost and quality in finite element analysis. *Finite Elements in Analysis and Design*, 41(15):1413–1440, 2005.

[3] R. Cottereau, P. Díez, and A. Huerta. Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method. *Computational Mechanics*, 44(4):533–547, 2009.

[4] P. Díez, N. Parés, and A. Huerta. Recovering lower bounds of the error by postprocessing implicit residual a posteriori error estimates. *International Journal for Numerical Methods in Engineering*, 56(10):1465–1488, 2003.

[5] Pedro Díez and Antonio Huerta. A unified approach to remeshing strategies for finite element h-adaptivity. *Computer Methods in Applied Mechanics and Engineering*, 176(1-4):215–229, 1999.

[6] Pedro Díez, Juan José Ródenas, and Olgierd C. Zienkiewicz. Equilibrated patch recovery error estimates: simple and accurate upper bounds of the error. *International Journal for Numerical Methods in Engineering*, 69(10):2075–2098, 2007.

[7] C. Farhat and F. X. Roux. The dual schur complement method with well-posed local neumann problems. *Contemporary Mathematics*, 157:193–201, 1994.

[8] Charbel Farhat and Jan Mandel. The two-level FETI method for static and dynamic plate problems part I: An optimal iterative solver for biharmonic systems. *Computer Methods in Applied Mechanics and Engineering*, 155(12):129–151, 1998.

[9] L. Gallimard. A constitutive relation error estimator based on traction-free recovery of the equilibrated stress. *International Journal for Numerical Methods in Engineering*, 78(4):460–482, 2009.

[10] P. Gosselet, C. Rey, and J. Pebrel. Total and selective reuse of krylov subspaces for the resolution of sequences of nonlinear structural problems. *International Journal for Numerical Methods in Engineering*, 94(1):60–83, 2013.

[11] Pierre Gosselet and Christian Rey. Non-overlapping domain decomposition methods in structural mechanics. *Archives of Computational Methods in Engineering*, 13(4):515–572, 2006.

[12] P. Ladevèze. Upper error bounds on calculated outputs of interest for linear and nonlinear structural problems. *Comptes Rendus Académie des Sciences - Mécanique, Paris*, 334(7):399–407, 2006.

[13] P. Ladevèze. Strict upper error bounds on computed outputs of interest in computational structural mechanics. *Computational Mechanics*, 42(2):271–286, 2008.

[14] P. Ladevèze and D. Leguillon. Error estimate procedure in the finite element method and application. *SIAM Journal of Numerical Analysis*, 20(3):485–509, 1983.
<table><caption>Table 4: Pre-cracked structure: First mesh : bounds for the exact quantity of interest</caption><thead><tr><th></th><th>I<sub>ex</sub><sup>+</sup></th><th>I<sub>ex</sub><sup>-</sup></th><th>width</th><th>precision</th></tr></thead><tbody><tr><td>With β<sub>inf</sub><sup>+</sup> = β<sub>inf</sub><sup>-</sup> = 0</td><td>3.8347</td><td>2.5791</td><td>1.2555</td><td>39.852 %</td></tr><tr><td>With β<sub>inf</sub><sup>+</sup> and β<sub>inf</sub><sup>-</sup> from table 3</td><td>3.5825</td><td>2.8824</td><td>0.7001</td><td>22.224 %</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>Figure 9: First mesh : Distribution of the error estimator within subdomains</figcaption></figure>

[15] P. Ladevèze and J. P. Pelle. *Mastering Calculations in Linear and Nonlinear Mechanics*. Springer, New York, 2004.

[16] Patrick Le Tallec. Domain decomposition methods in computational mechanics. *Comput. Mech. Adv.*, 1(2):121–220, 1994.

[17] J. T. Oden and S. Prudhomme. Goal-oriented error estimation and adaptivity for the finite element method. *Computers & Mathematics with Applications*, 41(5-6):735–756, 2001.

[18] S. Ohnimus, E. Stein, and E. Walhorn. Local error estimates of fem for displacements and stresses in linear elasticity by solving local neumann problems. *International Journal for Numerical Methods in Engineering*, 52(7):727–746, 2001.

[19] M. Paraschivoiu, J. Peraire, and A. T. Patera. A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations. *Computer Methods in Applied Mechanics and Engineering*, 150(1-4):289–312, 1997.

[20] N. Parés, P. Díez, and A. Huerta. Subdomain-based flux-free a posteriori error estimators. *Computer Methods in Applied Mechanics and Engineering*, 195(4-6):297–323, 2006.

[21] N. Parés, H. Santos, and P. Díez. Guaranteed energy error bounds for the poisson equation using a flux-free approach: Solving the local problems in subdomains. *International Journal for Numerical Methods in Engineering*, 79(10):1203–1244, 2009.

[22] A. Parret-Fréaud, C. Rey, P. Gosselet, and F. Feyel. Fast estimation of discretization error for fe problems solved by domain decomposition. *Computer Methods in Applied Mechanics and Engineering*, 199(49-52):3315–3323, 2010.

[23] F. Pled, L. Chamoin, and P. Ladevèze. On the techniques for constructing admissible stress fields in model verification: Performances on engineering examples. *International Journal for Numerical Methods in Engineering*, 88(5):409–441, 2011.

[24] S. Prudhomme and J. T. Oden. On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors. *Computer Methods in Applied Mechanics and Engineering*, 176(1-4):313–331, 1999.
<table><caption>Table 5: Pre-cracked structure: Second mesh</caption><thead><tr><th>Iteration</th><th>θ<sub>discr</sub></th><th>ρ<sub>discr</sub></th><th>θ̂<sub>discr</sub></th><th>ρ̃<sub>discr</sub></th><th>α</th><th>α̃</th></tr></thead><tbody><tr><td>0</td><td>5.7611</td><td>3.8882</td><td>7.5633 10<sup>-2</sup></td><td>7.551 10<sup>-4</sup></td><td>5.7017</td><td>7.4954 10<sup>-3</sup></td></tr><tr><td>1</td><td></td><td></td><td></td><td></td><td>3.3265</td><td>4.4442 10<sup>-3</sup></td></tr><tr><td>2</td><td></td><td></td><td></td><td></td><td>1.2869</td><td>8.4996 10<sup>-4</sup></td></tr><tr><td>3</td><td>5.4304</td><td>3.7882</td><td>7.5564 10<sup>-2</sup></td><td>7.6332 10<sup>-4</sup></td><td>0.31197</td><td>2.4632 10<sup>-4</sup></td></tr></tbody></table>

<table><caption>Table 6: Pre-cracked structure: Second mesh : error estimation</caption><thead><tr><th>I<sub>H</sub></th><th>β<sup>-</sup><sub>inf</sub></th><th>β<sup>+</sup><sub>inf</sub></th><th>1/4 β<sup>+</sup><sub>sup</sub></th><th>1/4 β<sup>-</sup><sub>sup</sub></th></tr></thead><tbody><tr><td>3.2266</td><td>0.32673</td><td>0.43838</td><td>0.23661</td><td>0.17372</td></tr></tbody></table>

[25] Christian Rey and Franck Risler. A rayleigh-ritz preconditioner for the iterative solution to large scale nonlinear problems. *Numerical Algorithms*, 17(3/4):279-311, 1998.

[26] V. Rey, P. Gosselet, and C. Rey. Study of the strong prolongation equation for the construction of statically admissible stress fields: Implementation and optimization. *Computer Methods in Applied Mechanics and Engineering*, 268(0):82-104, 2014.

[27] Valentine Rey, Pierre Gosselet, and Christian Rey. Strict bounding of quantities of interest in computations based on domain decomposition. *Computer Methods in Applied Mechanics and Engineering*, 287(0):212-228, 2015.

[28] Valentine Rey, Christian Rey, and Pierre Gosselet. A strict error bound with separated contributions of the discretization and of the iterative solver in non-overlapping domain decomposition methods. *Computer Methods in Applied Mechanics and Engineering*, 270(0):293-303, 2014.

[29] F. Risler and C. Rey. Iterative accelerating algorithms with krylov subspaces for the solution to large-scale non-linear problems. *Numerical algorithms*, 23:1-30, 2000.

[30] Y. Saad. *Iterative methods for sparse linear systems*. Society for Industrial and Applied Mathematics, Philadelphia, 3rd edition edition, 2003.

[31] Yousef Saad. Analysis of augmented krylov subspace methods. *SIAM Journal on Matrix Analysis and Applications*, 18(2):435-449, 1997.

[32] T. Strouboulis, I. Babuška, D. K. Datta, K. Copps, and S. K. Gangaraj. A posteriori estimation and adaptive control of the error in the quantity of interest. part I: A posteriori estimation of the error in the von mises stress and the stress intensity factor. *Computer Methods in Applied Mechanics and Engineering*, 181(1-3):261-294, 2000.

[33] R. Verfürth. *A Review of A Posteriori Error Estimation and Adaptive Mesh-refinement Techniques*. Wiley-Teubner, Stuttgart, 1996.
<table><caption>Table 7: Pre-cracked structure: Second mesh : bounds for the exact quantity of interest</caption><thead><tr><th></th><th><i>I</i><sub>ex</sub><sup>+</sup></th><th><i>I</i><sub>ex</sub><sup>-</sup></th><th>width</th><th>precision</th></tr></thead><tbody><tr><td>With <i>β</i><sub>inf</sub><sup>+</sup> = <i>β</i><sub>inf</sub><sup>-</sup> = 0</td><td>3.4632</td><td>3.0528</td><td>0.41034</td><td>12.717 %</td></tr><tr><td>With <i>β</i><sub>inf</sub><sup>+</sup> and <i>β</i><sub>inf</sub><sup>-</sup> from table 6</td><td>3.3815</td><td>3.1624</td><td>0.21906</td><td>6.7893 %</td></tr></tbody></table>

<table><caption>Table 8: Pre-cracked structure: Third mesh</caption><thead><tr><th>Iteration</th><th>θ<sub>discr</sub></th><th>ρ<sub>discr</sub></th><th>θ̃<sub>discr</sub></th><th>ρ̃<sub>discr</sub></th><th>α</th><th>α̃</th></tr></thead><tbody><tr><td>0</td><td>5.1501</td><td>3.5841</td><td>4.6197 10<sup>-2</sup></td><td>2.9485 10<sup>-4</sup></td><td>0.32764</td><td>6.2878 10<sup>-4</sup></td></tr><tr><td>1</td><td>5.1267</td><td>3.5824</td><td>4.619 10<sup>-2</sup></td><td>2.9482 10<sup>-4</sup></td><td>0.10619</td><td>1.6588 10<sup>-4</sup></td></tr></tbody></table>

<table><caption>Table 9: Pre-cracked structure: Third mesh : error estimation</caption><thead><tr><th><i>I</i><sub>H</sub></th><th><i>β</i><sub>inf</sub><sup>-</sup></th><th><i>β</i><sub>inf</sub><sup>+</sup></th><th><sup>1</sup>&frasl;<sub>4</sub> <i>β</i><sub>sup</sub><sup>+</sup></th><th><sup>1</sup>&frasl;<sub>4</sub> <i>β</i><sub>sup</sub><sup>-</sup></th></tr></thead><tbody><tr><td>3.2625</td><td>0.18986</td><td>0.25873</td><td>0.13618</td><td>0.10058</td></tr></tbody></table>

<table><caption>Table 10: Pre-cracked structure: Third mesh : bounds for the exact quantity of interest</caption><thead><tr><th></th><th><i>I</i><sub>ex</sub><sup>+</sup></th><th><i>I</i><sub>ex</sub><sup>-</sup></th><th>width</th><th>precision</th></tr></thead><tbody><tr><td>With <i>β</i><sub>inf</sub><sup>+</sup> = <i>β</i><sub>inf</sub><sup>-</sup> = 0</td><td>3.3986</td><td>3.1619</td><td>0.2367</td><td>7.2574 %</td></tr><tr><td>With <i>β</i><sub>inf</sub><sup>+</sup> and <i>β</i><sub>inf</sub><sup>-</sup> from table 9</td><td>3.3512</td><td>3.2265</td><td>0.1247</td><td>3.8199%</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Figure 10: Pre-cracked structure : Discretization envelope and residual against cumulative iterations</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Pre-cracked structure : Evolution of the approximated quantity of interest $I_H$ and of the bounds of the exact value $I_{ex}$</figcaption></figure>