O A González-Estrada

E Nadal

J J Ródenas
email: jjrodena@mcm.upv.es

P Kerfriden
email: kerfridenp@cardiff.ac.uk

S P A Bordas

F J Fuenmayor
email: ffuenmay@mcm.upv.es

Mesh adaptivity driven by goal-oriented locally equilibrated superconvergent patch recovery

Keywords: goal-oriented, error estimation, recovery, quantities of interest, error control, mesh adaptivity

Goal-oriented error estimates (GOEE) have become popular tools to quantify and control the local error in quantities of interest (QoI), which are often more pertinent than local errors in energy for design purposes (e.g. the mean stress or mean displacement in a particular area, the stress intensity factor for fracture problems). These GOEE are one of the key unsolved problems of advanced engineering applications in, for example, the aerospace industry. This work presents a simple recovery-based error estimation technique for QoIs whose main characteristic is the use of an enhanced version of the Superconvergent Patch Recovery (SPR) technique previously used for error estimation in the energy norm. This enhanced SPR technique is used to recover both the primal and dual solutions. It provides a nearly statically admissible stress field that results in accurate estimations of the local contributions to the discretisation error in the QoI and, therefore, in an accurate estimation of this magnitude. This approach leads to a technique with a reasonable computational cost that could easily be implemented into already available finite element codes, or as an independent postprocessing tool.

Introduction

Assessing the quality of numerical simulations has been a critical area of research for many years. In the context of the finite element method, verification tools have been extensively developed since the late 70s. Historically, these tools have aimed at quantifying a global distance between the finite element solution and the (unknown) exact solution. In the context of linear elasticity, this distance is naturally measured in the energy norm, as the optimal properties underlying the finite element method are established in this sense. Adaptivity schemes have been developed to control this global error measure by local mesh refinement. However, the use of energy estimates for verification is rather limited. Indeed, this particular distance does not necessarily reflect the precision of the actual output of the calculation, for instance, a local stress measure or the displacement of part of the structure. In the last 15 years, a revolution has happened in this area of research: the development of the so-called error estimates in quantities of interest. The idea is to quantify the effect of local errors on the precision of the output of the calculation, which leads to adaptivity schemes and quality assessment that are more relevant to engineering applications. Technically, errors in quantities of interest rely on the evaluation of the quantity of interest by means of solving an auxiliary "mechanical" problem. Once the quantity of interest is globalised, a combination of energy estimates for the initial problem and auxiliary problem can be deployed to estimate the error in the quantity of interest.

The energy estimates developed during the last 40 years can be classified into two categories: explicit estimates that try to link a certain measure of the residual to the error in the displacement field, and estimates that post-process solution fields that are somehow more accurate than the finite element solution.

The first category of estimates mimics the optimal bounds used to prove the convergence of finite element discretisation schemes. Very easy to implement, the explicit estimates include interpolation constants that are problem dependent and difficult to evaluate. For this reason, they are usually not favoured by engineers. Energy estimates based on post-processing reconstruct an approximation of the exact solution by performing local computations. A first family of such estimates, usually called implicit residual-based approaches, rely on multilevel finite element approaches. A finite element problem with a very fine discretisation is solved over each of the local subdomains (either individual elements [START_REF] Ainsworth | A posteriori Error Estimation in Finite Element Analysis[END_REF], patches of elements [START_REF] Babuška | A-posteriori error estimates for the finite element method[END_REF][START_REF] Cottereau | Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method[END_REF] or subdomains consisting of an element and its immediate neighbourhood [START_REF] Larsson | Strategies for computing goal-oriented a posteriori error measures in non-linear elasticity[END_REF]). Depending on how the local problem is linked to the global finite element solution whose accuracy is to be evaluated, different properties can be obtained for the estimate. For instance, the equilibrated element residual approach, the fluxfree approach and the constitutive relation error, yield estimates that bound the error from above. Error estimates based on local elements with Dirichlet boundary conditions bound the error from below (see for instance [START_REF] Díez | Recovering lower bounds of the error by postprocessing implicit residual a posteriori error estimates[END_REF]). Other implicit residualbased approaches tailor the link between the local and global problems so that the efficiency of the estimate is maximised, without focussing on the bounding properties [START_REF] Larsson | Strategies for computing goal-oriented a posteriori error measures in non-linear elasticity[END_REF]. More detailed discussion about this class of estimates can be found in [START_REF] Verfürth | A review of a posteriori error estimation and adaptive meshrefinement techniques[END_REF][START_REF] Stein | Error-controlled adaptive finite elements in solid mechanics[END_REF][START_REF] Díez | A posteriori error estimation for standard finite element analysis[END_REF].

The second post-processing approach consists of deriving simple smoothing techniques that yield a solution field that converges faster than the finite element solution. In the context of the finite element method, a very popular prototype for such approaches is the Zienkiewicz-Zhu estimate (ZZ) [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF], associated, for example, with nodal averaging or with the superconvergent patch recovery (SPR) [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique[END_REF]. In the latter case, the idea is to fit a higher-order polynomial approximation of the stress field over each patch of elements. The distance in the energy norm between this smoothed stress field and the finite element stress field is defined as an error estimate. The success of this approach in the engineering community relies on an intuitive mechanical definition and a certain ease of implementation compared to the first class of estimates, without sacrificing the numerical effectivity. These smoothing techniques were extended to enriched approximations in [START_REF] Bordas | Derivative recovery and a posteriori error estimate for extended finite elements[END_REF][START_REF] Bordas | A simple error estimator for extended finite elements[END_REF][START_REF] Ródenas | A recoverytype error estimator for the extended finite element method based on singu-lar+smooth stress field splitting[END_REF][START_REF] Duflot | A posteriori error estimation for extended finite elements by an extended global recovery[END_REF][START_REF] Prange | Error estimation for crack simulations using the XFEM[END_REF][START_REF] Hild | A residual error estimator for the XFEM approximation of the elasticity problem[END_REF] and to smoothed finite elements (SFEM) in [START_REF] González-Estrada | Efficient recovery-based error estimation for the smoothed finite element method for smooth and singular linear elasticity[END_REF]. The role of enrichment and statical admissibility in such procedures was discussed in [START_REF] González-Estrada | On the role of enrichment and statical admissibility of recovered fields in aposteriori error estimation for enriched finite element methods[END_REF]. Some contributions have suggested bridging these two approaches, mostly in an attempt to obtain the guaranteed upper bound that some of the residual-based techniques propose, while retaining the ease of implementation of the ZZ framework. The basic identity is that when the recovered stress field is exactly statically admissible, then the ZZ estimate coincides with the constitutive relation error, and bounds the energy error from above. Following this line of thought, Díez et al. [START_REF] Díez | Equilibrated patch recovery error estimates: simple and accurate upper bounds of the error[END_REF] and Ródenas et al. [START_REF] Ródenas | Accurate recoverybased upper error bounds for the extended finite element framework[END_REF] presented a methodology to obtain practical upper bounds of the error in the energy norm using an SPR-based approach where equilibrium was locally imposed on each patch. Other approaches propose a partial fulfilment of the equilibrium to improve the robustness of the smoothing procedure, which can be seen as a trade-off between exactly equilibrated approaches and more simple empirical averaging. This contribution considers the ZZ-type estimate and the use of one of these smoothing techniques, namely SPR-CX (where C stands for the enforcement of equilibrium constraints and X for the smooth+singular decomposition of the recovered stress field). In this approach, the internal equilibrium and boundary equilibrium are constrained to be locally satisfied in the smoothing-based recovery operation. Such an approach has been shown to yield very good practical global and local effectivities, in particular on the boundaries where basic averaging techniques fail to obtain accurate results. Note that the method is conceptually and practically different from what is proposed by authors working on equilibrated residual (ER) approaches [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF][START_REF] Stein | Error-controlled adaptive finite elements in solid mechanics[END_REF][START_REF] Cottereau | Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method[END_REF] or constitutive relation error methods [START_REF] Ladevèze | Upper error bounds on calculated outputs of interest for linear and nonlinear structural problems[END_REF]. Conceptually, methods using the ER concept focus on upper bounding properties, whilst trying to achieve good effectivities. This is achieved by recovering a globally equilibrated stress field. Here, our focus is effectivity only. The patch-equilibration is performed to ensure that the recovered stress field converges faster than the finite element solution. It is not required to satisfy global equilibrium. Practically, most ER techniques solve for local equilibria using polynomial finite elements, yielding asymptotic error bounds. The use of local divergence-conforming approximation spaces associated with energy minimisation principles [START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF][START_REF] Cottereau | Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method[END_REF] permits to achieve strict error bounding. In our case, we aim at achieving superconvergence by fitting our polynomial approximation to the finite element solution at particular points of the patch. The equilibrium is introduced as a constraint to this minimisation problem.

As mentioned previously, the error in quantities of interest can be related to the energy estimates by means of auxiliary problems. This was done successfully in a number of cases, including implicit residual approaches [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF][START_REF] Rüter | Goal-oriented a posteriori error estimates in linear elastic fracture mechanics[END_REF][START_REF] Larsson | Strategies for computing goal-oriented a posteriori error measures in non-linear elasticity[END_REF][START_REF] Ladevèze | Local error estimators for finite element linear analysis[END_REF], dual analysis [START_REF] De Almeida | Upper bounds of the error in local quantities using equilibrated and compatible finite element solutions for linear elastic problems[END_REF] and the ZZ smoothing approaches [START_REF] Cirak | A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem[END_REF][START_REF] Rüter | Goal-oriented a posteriori error estimates in linear elastic fracture mechanics[END_REF] with basic nodal averaging or SPR technique. Note that if an energy estimate with bounding properties is used, then a similar bounding of the error in the quantity of interest can be obtained [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF][START_REF] Ladevèze | Local error estimators for finite element linear analysis[END_REF]. However, it is usually difficult to obtain guaranteed error bounds in local quantities of interest while maintaining the effectivity of the estimate. The necessity of obtaining such bounds in an engineering context is also arguable, as the reliability of an a posteriori error estimate, which is quantified by its local effectivity, can be verified beforehand on a number of practical cases.

This contribution proposes to use the SPR-CX approach to derive an efficient and simple goal-oriented adaptivity procedure in the context of elasticity for linear quantities of interest. In this context, a simple stress smoothing procedure with partial fulfilment of the equilibrium is applied to both the initial problem and the auxiliary problem associated with the quantity of interest. As the enhanced smoothing technique relies on the explicit definition of the load, we emphasise that the auxiliary or dual problem can often be seen as a mechanical problem with body and boundary loads, and prescribed initial stresses and strains that depend on the particular quantity of interest. We derive these mechanical loads for a number of practical quantities of interest, in relation with the SPR-CX smoothing technique. We also extend this approach to handle singular elasticity problems, where the quantity of interest is a generalised stress intensity factor. The methodology is verified in a number of numerical examples. We also implement the hierarchy of practical upper bounds proposed in [START_REF] Rüter | Goal-oriented a posteriori error estimates in linear elastic fracture mechanics[END_REF] with the ZZ/SPR-CX technique as the energy estimate. We show that the respective property that we observe for each of the estimates follows the trend observed by the authors of that work.

The paper is organised as follows. In Section 2 we define the model problem and the finite element approximation used. Section 3 is related to the representation of the error in the energy norm and it is used to introduce the nearly equilibrated recovery technique for obtaining enhanced stress fields. In Section 4 we provide the basics behind error estimation in quantities of interest. We define the dual problem and present the smoothing-based error estimate, in particular, the expressions for the dual loads required for the stress recovery procedure. We present the smoothingbased error estimate and describe the expressions for the dual loads, which are required for the stress recovery process. In Section 5 we present some numerical examples and, finally, conclusions are drawn in Section 6.

2 Linear elasticity problem solved by the finite element method

Problem statement

In this section, we introduce the 2D linear elasticity problem to which the proposed methodology is dedicated. Using Voigt notation, we denote by σ = {σ x , σ y , σ xy } T the Cauchy stress, by u the displacement, and by ε the strain, all these fields being defined over domain Ω ⊂ R 2 , of boundary denoted by ∂Ω. Prescribed tractions, denoted by t, are imposed over part Γ N of the boundary, while displacements denoted by ū are prescribed over the complementary part Γ D of the boundary. b denotes the body load. The elasticity problem takes the following form. We look for (σ, u) satisfying:

• statical admissibility

L T σ + b = 0 in Ω (1) Gσ = t on Γ N ( 2
)
where L is the differential operator for linear elasticity, and G is the projection operator that projects the stress field into tractions over the boundary considering the unit normal n to Γ N such that

L T = ∂/∂x 0 ∂/∂y 0 ∂/∂y ∂/∂x , G = n x 0 n y 0 n y n x • kinematic admissibility u = ū on Γ D , (3)
• constitutive relation

ε(u) = Lu in Ω (4) σ = D(ε(u) -ε 0 ) + σ 0 in Ω (5)
where operator D contains the elasticity coefficients of the usual linear isotropic constitutive law relating stress and strain, and σ 0 and ε 0 are initial stresses and strains respectively.

Using the notations introduced in [START_REF] Verfürth | A review of a posteriori error estimation techniques for elasticity problems[END_REF] the problem above takes the primal variational form:

Find u ∈ (V + {w}) such that ∀v ∈ V : Ω ε(u) T Dε(v)dΩ = Ω v T bdΩ + Γ N v T tdΓ + Ω ε(v) T Dε 0 dΩ - Ω ε(v) T σ 0 dΩ (6)
where

V = {v | v ∈ [H 1 (Ω)] 2 , v| Γ D = 0}
and w is a particular displacement field satisfying the Dirichlet boundary conditions.

Finite element discretisation

Let us introduce a classical finite element discretisation scheme for the elasticity problem. The approximate displacement field u h is searched for in a space of finite dimension (V h + {w}) ⊂ (V + {w}) such that V h is spanned by locally supported finite element shape functions.

Using the Galerkin framework, the primal variational formulation ( 6) is recast in the form:

Find u h ∈ V h + {w} such that ∀v ∈ V h : Ω ε(u h ) T Dε(v)dΩ = Ω v T bdΩ + Γ N v T tdΓ + Ω ε(v) T Dε 0 dΩ - Ω ε(v) T σ 0 dΩ (7)
which can be solved using classical finite element technology [START_REF] Zienkiewicz | The Finite Element Method: The Basis[END_REF].

3 Smoothing-based error estimates in the energy norm

Zienkiewicz-Zhu estimate

In the absence of other types of errors, the finite element discretisation error is defined by e = uu h . To quantify the quality of u h one of the standard approaches is to evaluate the energy norm of e, which is defined by:

e 2 = Ω ε(e) T Dε(e)dΩ. (8)
If we introduce the error in the stress field e σ = σ-σ h , where σ h = D ε(u h ) -ε 0 + σ 0 is the finite element stress field, and make use of the constitutive relation, the previous expression can be written as

e 2 = Ω e T σ D -1 e σ dΩ (9)
Following Zienkiewicz-Zhu [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF], an estimate E zz of the exact error measure e can be obtained by introducing the approximation

e 2 ≈ (E zz ) 2 = Ω (e * σ ) T D -1 (e * σ ) dΩ, ( 10
)
where the approximate stress error e * σ is defined by e * σ = σ * -σ h and σ * is the recovered stress field.

If σ * converges to the exact solution at a higher rate (superconvergent) than the FE stress solution, then, the ZZ estimate is asymptotically exact, meaning that the approximate error tends towards the exact error with mesh refinement [START_REF] Zienkiewicz | The Finite Element Method: The Basis[END_REF]. This estimate does not have guaranteed bounding properties unless the recovered field is statically admissible, in which case it coincides (formally) with the Constitutive Relation Error [START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF][START_REF] Pled | An enhanced method with local energy minimization for the robust a posteriori construction of equilibrated stress fields in finite element analyses[END_REF] and the Equilibrated Residual Approach [START_REF] Ainsworth | A posteriori Error Estimation in Finite Element Analysis[END_REF], which are considered difficult to implement. The most appealing feature of the ZZ approach to error estimation is that very simple recovery techniques based on stress smoothing permit to obtain good effectivities, which explains the popularity of this method in the engineering community. However, the basic smoothing-based recovery techniques, such as the original SPR, which only use the FE results for the recovery, provide good global effectivities but suffer from a local lack of accuracy. For example, the SPR does not use the information along the boundary, where the imposed tractions are known and the patches can have a reduced number of elements, thus, leading to a well-known loss in accuracy of the recovered stresses along the Neumann boundary. This local lack of accuracy is a major issue in engineering applications, where boundary stresses are typically of interest, and makes such basic estimates difficult to use to drive adaptivity. To tackle this problem, trade-offs between exactly equilibrated approaches and smoothing-techniques were developed. The basics of one such method, the SPR-CX technique, are recalled in the next section as they constitute the corner stone to derive the efficient error estimates in quantities of interest presented in this contribution.

Enhanced SPR-based stress recovery

Stress smoothing

As noted in previous sections, a widely used technique to control the error in the energy norm in the finite element discretisation is the Zienkiewicz-Zhu error estimator shown in [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique[END_REF], which is based on the recovery (often called smoothing) of an enhanced stress field σ * .

To obtain the recovered field σ * , first, we define the field σ -such that we subtract the initial stress and strain from the field σ:

σ -= σ -σ 0 + Dε 0 , (11)
and perform the smoothing on σ -. Then, the recovered field is

σ * = (σ -) * + σ 0 -Dε 0 , (12)
where (σ -) * is the smoothed field that corresponds to σ -. In the SPR-CX technique, a patch P (J) is defined as the set of elements connected to a vertex node J. On each patch, a polynomial expansion for each one of the components of the recovered stress field is expressed in the form

σ * k (x) = p(x)a k k = xx, yy, xy (13)
where p represents a polynomial basis and a k are unknown coefficients to be evaluated. Usually, the polynomial basis is chosen equal to the finite element basis for the displacements. The coefficients a k are evaluated using a least squares approximation to the values of FE stresses evaluated at the integration points of the elements within the patch, x G ∈ P (J) . The recovered stress field coupling the three stress components for the 2D case reads:

σ * (x) =    σ * xx (x) σ * yy (x) σ * xy (x)    = P(x)A =   p(x) 0 0 0 p(x) 0 0 0 p(x)      a xx a yy a xy    (14)
To obtain the coefficients A we solve a linear system of equations resulting from the minimisation of the functional

F (J) (A) = P (J) (PA -σ -h ) 2 dΩ ( 15
)
where σ -h = Dε(u h ).

To obtain a continuous field, a partition of unity procedure [START_REF] Blacker | Superconvergent patch recovery with equilibrium and conjoint interpolant enhancements[END_REF] properly weighting the stress interpolation polynomials, obtained from the different patches formed at the vertex nodes of the element containing point x, is used. The field σ * is interpolated using linear shape functions N (J) associated with the n v vertex nodes such that

σ * (x) = nv J=1 N (J) (x) σ * (J) (x) -Dε 0 (x) + σ 0 (x). ( 16
)
Note that in [START_REF] Hild | A residual error estimator for the XFEM approximation of the elasticity problem[END_REF] we add back the contribution of the initial stresses and strains subtracted in [START_REF] Bordas | Derivative recovery and a posteriori error estimate for extended finite elements[END_REF].

Equilibrium conditions

The accuracy of such estimates depends on the quality of the recovered field. To improve the quality of the recovered fields, numerical results indicate that when recovering singular fields, statical admissibility and suitably chosen enrichment functions improve the accuracy [START_REF] González-Estrada | On the role of enrichment and statical admissibility of recovered fields in aposteriori error estimation for enriched finite element methods[END_REF][START_REF] González-Estrada | Error estimation in quantities of interest for XFEM using recovery techniques[END_REF] 1 . In this work we consider the SPR-CX recovery technique, which is an enhancement of the error estimator introduced in [START_REF] Díez | Equilibrated patch recovery error estimates: simple and accurate upper bounds of the error[END_REF], to recover the solutions for both the primal and dual problems ( see Section 4). The technique incorporates the ideas in [START_REF] Ródenas | Improvement of the superconvergent patch recovery technique by the use of constraint equations: the SPR-C technique[END_REF] to guarantee locally on patches the exact satisfaction of the equilibrium equations, and the extension in [START_REF] Ródenas | A recoverytype error estimator for the extended finite element method based on singu-lar+smooth stress field splitting[END_REF] to singular problems. Constraint equations are introduced via Lagrange multipliers into the linear system used to solve for the coefficients A on each patch in order to enforce the satisfaction of the:

• Internal equilibrium equations: We define the constraint equation for the internal equilibrium in the patch as:

∀x j ∈ P (J) L T σ * (J) (x j ) + L T (σ 0 (x j ) -Dε 0 (x j )) + b(x j ) := c int (x j ) = 0 (17)
We consider b(x) a polynomial least squares fit of degree p -1, p being the degree of the recovered stress field σ * (J) , to the actual body forces b(x). We enforce c int (x j ) at a sufficient number of j non-aligned points (nie) to guarantee the exact representation of b(x) and solve the linear system. For example, in the case of b(x) being represented by a plane, we choose a set of three nonaligned points in P (J) . In 2D, this constraint will add two equations (one per spatial direction) per point to the linear system assembled for the patch, .i.e. 2nie equations. 1 The use of enrichment to improve recovery based error estimates for enriched approximations was discussed in detail in some of the first papers discussing derivative recovery techniques for enriched finite element approximations, see References [START_REF] Ródenas | Accurate recoverybased upper error bounds for the extended finite element framework[END_REF][START_REF] Bordas | Derivative recovery and a posteriori error estimate for extended finite elements[END_REF][START_REF] Bordas | A simple error estimator for extended finite elements[END_REF][START_REF] Ródenas | A recoverytype error estimator for the extended finite element method based on singu-lar+smooth stress field splitting[END_REF][START_REF] Duflot | A posteriori error estimation for extended finite elements by an extended global recovery[END_REF][START_REF] González-Estrada | On the role of enrichment and statical admissibility of recovered fields in aposteriori error estimation for enriched finite element methods[END_REF]. A very detailed and clear discussion of a wide variety of error estimators and adaptive procedures for discontinuous failure is provided in [START_REF] Pannachet | Error estimation and adaptivity for discontinuous failure[END_REF].

• Boundary equilibrium equations: A point collocation approach is used to impose the satisfaction of a polynomial approximation to the tractions along the Neumann boundary intersecting the patch. The constraint equation reads

∀x j ∈ Γ N ∩P (J) G σ * (J) (x j )+GL T (σ 0 (x j )-Dε 0 (x j ))-t(x j ) := c ext (x j ) = 0 (18)
We enforce c ext (x j ) in nbe = p + 1 points along the part of the boundary crossing the patch. In the case that more than one boundary intersects the patch, only one curve is considered in order to avoid over-constraining. The boundary constraint adds 2nbe equations to the linear system of equations to solve.

• Compatibility equations: c cmp (x j ) is only imposed in the case that p ≥ 2 in a sufficient number of non-aligned points nc. For example, for p = 2 we have nc = 1. σ * directly satisfies c cmp for p = 1. The 2D compatibility equation expressed in terms of stresses is:

∀x j ∈ P (J) ∂ 2 ∂y 2 (kσ xx -νqσ yy ) + ∂ 2 ∂x 2 (kσ yy -νqσ xx ) -2(1 + ν) ∂ 2 σxy ∂x∂y = 0 (19) where k, q are functions of the Poisson's coefficient ν k = 1, q = 1 for plane stress k = (1 -ν) 2 , q = (1 + ν) for plane strain
Enforcement of compatibility adds nc equations to the linear system for the patch.

Thus, after adding the 2nie + nc equations and the additional 2nbe equations in patches intersected by the Neumann boundary, the Lagrange functional enforcing the constraint equations for a patch P (J) can be written as

L (J) (A, λ) = F (J) (A) + nie i=1 λ int i c int (x i ) + nbe j=1 λ ext j c ext (x j ) + nc k=1 λ cmp k (c cmp (x k )) (20)
Problems involving an internal interface Γ I along which continuity of the field is not maintained (e.g. bimaterial problems), we use a different polynomial expansion on each side of the boundary and enforce statical admissibility (traction continuity along Γ I ). Suppose that we have a patch intersected by Γ I such that Ω e = Ω 1,e ∪ Ω 2,e for intersected elements, as indicated in Figure 1. To enforce equilibrium conditions along Γ I we define the stresses σ *

Ω 1 | Γ I , σ * Ω 2 | Γ I on
each side of the internal boundary. Then, the boundary equilibrium along Γ I given the prescribed tractions t Γ I = [t x t y ] T is:

G( σ * Ω 1 | Γ I -σ * Ω 2 | Γ I ) = t Γ I (21)
Once we have the equilibrated recovered fields on each of the patches σ * (J) , ( 16) yields a continuous recovered stress field σ * . Note that this postprocess of the local solutions introduces in the global solution σ * a small lack of equilibrium s = nv J=1 ∇N (J) σ * (J) when evaluating the divergence of the internal equilibrium equation, as explained in [START_REF] Díez | Equilibrated patch recovery error estimates: simple and accurate upper bounds of the error[END_REF][START_REF] Ródenas | Accurate recoverybased upper error bounds for the extended finite element framework[END_REF].

Singular fields

Different techniques were proposed to account for the eventual non smooth, e.g. singular, part during the recovery process [START_REF] Ródenas | A recoverytype error estimator for the extended finite element method based on singu-lar+smooth stress field splitting[END_REF][START_REF] Bordas | Derivative recovery and a posteriori error estimate for extended finite elements[END_REF][START_REF] Duflot | A posteriori error estimation for extended finite elements by an extended global recovery[END_REF]. Here, following the ideas in [START_REF] Ródenas | A recoverytype error estimator for the extended finite element method based on singu-lar+smooth stress field splitting[END_REF] for singular problems, the exact stress field σ is decomposed into two stress fields, a smooth field σ smo and a singular field σ sing :

σ = σ smo + σ sing . (22)
Then, on each patch, the recovered field σ * required to compute the error estimate given in [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique[END_REF] can be expressed as the contribution of two recovered stress fields, one smooth σ * smo and one singular σ * sing :

σ * = σ * smo + σ * sing . ( 23
)
For the recovery of the singular part, the expression which describes the asymptotic fields near the crack tip with respect to a coordinate system (r, φ) at the tip is used [START_REF] Szabó | Finite Element Analysis[END_REF]:

σ sing (r, φ) = K I λ I r λ I -1 Φ I (λ I , φ) + K II λ II r λ II -1 Φ II (λ II , φ) (24)
where λ m (with m = I, II) are the eigenvalues that determine the order of the singularity, Φ m are a set of trigonometric functions that depend on the angular position φ and K m are the so-called generalised stress intensity factors (GSIFs) [START_REF] Szabó | Finite Element Analysis[END_REF].

To evaluate σ * sing we first obtain estimated values of the generalised stress intensity factors K I and K II using a domain integral method based on extraction functions [START_REF] Szabó | Finite Element Analysis[END_REF][START_REF] Giner | A domain integral for the calculation of generalized stress intensity factors in sliding complete contacts[END_REF]. The estimated values of K I and K II are evaluated once for the singularity and then used to define σ * on the patches. Notice that the recovered part σ * sing is an equilibrated field as it satisfies the equilibrium equations.

Once the field σ * sing is evaluated, an FE-type (discontinuous) approximation to the smooth part σh smo can be obtained subtracting σ * sing from the raw FE field:

σh smo = σh -σ * sing . (25)
Then, the field σ * smo is evaluated applying the enhancements of the SPR technique previously described, i.e. satisfaction of equilibrium and compatibility equations on each patch. Note that as both σ * smo and σ * sing satisfy the equilibrium equations, σ * also satisfies equilibrium on each patch. The splitting procedure is only applied on patches within a prescribed distance to the singularity as, away from the singular point, the solution is smooth enough and can be accurately recovered by the smoothing process.

4 Error estimation in quantities of interest by the SPR-CX recovery technique

Quantity of interest

In this section we show how the ZZ estimate may be extended to evaluate the error in a particular quantity of interest [START_REF] Ainsworth | A posteriori Error Estimation in Finite Element Analysis[END_REF]. For simplicity, in this section we assume that the Dirichlet boundary conditions are homogeneous. We will consider a quantity of interest Q : V → R, defined as a linear functional of the displacement field2 . For example, this quantity can be an average of the displacement over a part of the domain. The aim is to estimate the error in functional Q, which is expressed by

Q(u) -Q(u h ) = Q(u -u h ) = Q(e). ( 26
)

Auxiliary problem and exact error representation

Standard procedures to evaluate Q(e) introduce the dual problem (this terminology comes from the optimal control community) or auxiliary problem

Find ũ ∈ V such that ∀v ∈ V, Ω ε(v) T Dε(ũ)dΩ = Q(v). ( 27
)
Problem ( 27) can be seen as the variational form of an auxiliary mechanical problem. Dual field ũ ∈ V is a displacement that vanishes over Γ D . Test function v is a virtual displacement. Field σ = D(ε(ũ) -ε0 ) + σ0 , where σ0 and ε0 are known quantities (initial stress and strain) that will be detailed later on, can be interpreted as a mechanical stress field. The left-hand side of ( 27) is the work of internal forces of the auxiliary mechanical problem. As detailed later on, Q(v) is the work of an abstract external load for the auxiliary mechanical problem. The primal and dual problem are solved by the finite element method. Here, we will make use of the same finite element space for both problems. Therefore, we will look for an approximation of ũ ∈ V using the Galerkin approach:

Find ũh ∈ V h such that ∀v ∈ V h , Ω ε(v) T Dε(ũ h )dΩ = Q(v). ( 28
)
By choosing v = e in [START_REF] Cirak | A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem[END_REF], and making use of the Galerkin orthogonality for the primal problem (see [START_REF] Ainsworth | A posteriori Error Estimation in Finite Element Analysis[END_REF] for more details), we can recast expression [START_REF] Cirak | A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem[END_REF] as follows:

Q(e) = Ω ε(e) T Dε(ẽ)dΩ (29)
where ẽ = ũũh is the discretisation error of the dual problem [START_REF] Cirak | A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem[END_REF]. In the context of the ZZ framework, we make use of the constitutive relation in order to obtain an expression in terms of mechanical stresses:

Q(e) = Ω e T σ D -1 ẽσ dΩ ( 30
)
where ẽσ = σ-σh is the stress error of the dual problem and σh = D(ε(ũ h )-ε 0 )+ σ0 the finite element stress field.

Smoothing-based error estimate

In order to obtain a computable estimate of the error in the quantity of interest using the ZZ methodology, the stress error fields in expression [START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF] are replaced by their smoothing-based approximations:

Q(e) ≈ E 1 = Ω e * σ T D -1 ẽ * σ dΩ (31)
where the approximate dual error is ẽ * σ = σ * -σh and σ * is the recovered auxiliary stress field.

The recovered auxiliary stress field can be computed in different ways, for instance by using the standard SPR approach, as proposed in [START_REF] Cirak | A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem[END_REF][START_REF] Rüter | Goal-oriented a posteriori error estimates in linear elastic fracture mechanics[END_REF]. Here, we propose to recover both the primal and dual stress field by the SPR-CX technique described in Section 3.2. Two remarks are worth being made here. First, the analytical expressions that define the tractions and body forces for the dual problem are obtained from the interpretation of the functional Q in terms of tractions, body loads, initial stresses and strains, as seen in Section 4.4. Second, to enforce equilibrium conditions during the recovery process along the boundary of the domain of interest (DoI), we consider it as an internal interface. We use different polynomial expansions on each side of the boundary and enforce statical admissibility of the normal and tangential stresses as previously explained in Section 3.2.

Analytical definitions in the dual problem for enforcing the partial equilibrium

The SPR-CX procedure relies on the enforcement of equilibrium for the recovered stress fields evaluated on each patch. For this technique to be applied to recover the dual stress field, the corresponding mechanical equilibrium must be made explicit.

In order to do so, the right-hand side of ( 27) is interpreted as the work of mechanical external forces, and the analytical expression of these forces is derived, depending on the quantity of interest:

Find ũ ∈ V such that ∀v ∈ V : Ω ε(v) T Dε(ũ)dΩ = Q(v) = Ω v T bdΩ + Γ N v T tdΓ + Ω ε(v) T Dε 0 dΩ - Ω ε(v) T σ0 dΩ (32)
The problem in ( 32) is solved using a FE approximation with test and trial functions in V h . The finite element solution is denoted by ũh ∈ V h . Such derivations were presented in [START_REF] Ródenas | Goal Oriented Adaptivity: Una introducción a través del problema elástico lineal[END_REF][START_REF] González-Estrada | Equilibrated patch recovery for accurate evaluation of upper error bounds in quantities of interest[END_REF][START_REF] Verdugo | Natural quantities of interest in linear elastodynamics for goal oriented error estimation and adaptivity[END_REF]. Here, we only recall some of the results presented in these papers. Additionally, we provide the analytical expression of the dual load when the quantity of interest is the generalised stress intensity factor (GSIF).

Mean displacement in Ω I

Let us assume that the objective is to evaluate a weighted average of the displacement in a sub-domain of interest Ω I ⊂ Ω. In this case, the quantity of interest is

Q(u) = 1 |Ω I | Ω I u T c u dΩ (33)
where |Ω I | is the measure of Ω I and c u is an extraction vector used to select the combination of the components of the displacement field that are of interest. For example, c u = {1, 0} T to extract the first component of u. Now, the right-hand side of expression [START_REF] Cirak | A posteriori error estimation and adaptivity for linear elasticity using the reciprocal theorem[END_REF] becomes

Q(v) = Ω I v T c u |Ω I | dΩ (34)
By comparing this expression to the variational form (32) of a general mechanical problem, it is clear that quantity b defined by b = c u /|Ω I | corresponds to a field of body forces for the auxiliary mechanical problem.

Mean displacement along Γ I

The quantity of interest is now the mean value of the displacement along a subset Γ I of the Neumann boundary Γ N :

Q(u) = 1 |Γ I | Γ I u T c u dΓ ( 35
)
where |Γ I | is the length of Γ I and c u an extractor acting on u. Note that the quantity t = c u /|Γ I | can be interpreted as a vector of tractions applied along the boundary in the problem defined in [START_REF] Blacker | Superconvergent patch recovery with equilibrium and conjoint interpolant enhancements[END_REF]. Thus, t is a vector of tractions applied on Γ I which can be used in the dual problem to extract the mean displacements along Γ I .

Mean strain in Ω I

In this case we are interested in some combination of the components of the strain over a subdomain Ω I such that the QoI is given by:

Q(u) = 1 |Ω I | Ω I c T ε ε(u)dΩ = Ω I c T ε |Ω I | ε(u)dΩ (36)
where c ε is the extraction operator used to define the combination of strains under consideration. Thus, the term σ0 = c T ε /|Ω I | represents the vector of initial stresses that are used to define the auxiliary problem for this particular QoI.

Mean stress value in Ω I

Let us consider now as QoI the mean stress value given by a combination c σ of the stress components σ = D(ε(u) -ε 0 ) + σ 0 in a domain of interest which reads:

Q(u) = 1 |Ω I | Ω I c T σ (D(ε(u) -ε 0 ) + σ 0 )dΩ. ( 37
)
Q is an affine functional. Let us define

Q(v) = Ω c T σ D(ε(v))dΩ (38)
for v an arbitrary vector of H 1 (Ω). Q is such that Q(e) = Q(e), so that by solving the dual problem

Ω ε(v) T Dε(ũ)dΩ = Q(v) (39)
for ũ, the derivations of Section 4.2 apply.

Similarly to the previous quantity, the right-hand side of the auxiliary problem is defined by the term ε0 = c T σ /|Ω I |, which represents in this case a vector of initial strains.

Mean tractions along Γ I included in Γ D

Let us assume that we want to evaluate, for example, the mean value of a combination of the tractions T R on a part Γ I of the Dirichlet boundary Γ D . The application of the principle of virtual work with test functions v ∈ H 1 that do not necessarily vanish over Γ D gives:

Γ D v T T R dΓ = Ω ε(u) T Dε(v)dΩ - Γ N v T tdΓ - Ω v T bdΩ + Ω ε(v) T σ 0 dΩ - Ω ε(v) T Dε 0 dΩ (40)
Extracting the quantity c T R T R dΓ , where c R is an extractor defined over Γ I , is done by defining the prolongation δ ∈ H 1 (Ω) of extractor c R such that δ| Γ I = c R /|Γ I |. For instance, δ is the finite element field that is null at every node that does not belong to Γ I . Substituting δ for v in [START_REF] Verdugo | Natural quantities of interest in linear elastodynamics for goal oriented error estimation and adaptivity[END_REF] yields

Q(u) = Ω ε(u) T Dε(δ)dΩ - Γ N δ T tdΓ - Ω δ T bdΩ + Ω ε(δ) T σ 0 dΩ - Ω ε(δ) T Dε 0 dΩ (41)
Q is an affine functional. We define

Q(v) = Ω ε(v) T Dε(δ)dΩ (42)
Q is such that Q(e) = Q(e), so that by solving the dual problem for ũ

Ω ε(v) T Dε(ũ)dΩ = Q(v) (43)
the derivations of Section 4.2 applies. The dual load is an initial strain ε(δ). By recalling the variational form in [START_REF] Verfürth | A review of a posteriori error estimation and adaptive meshrefinement techniques[END_REF] we see that the dual problem is a mechanical problem in ū = ũδ, where ū|

Γ D = -δ| Γ D = -c R /|Γ I |,

Generalised stress intensity factor

Consider the problem of evaluating the generalised stress intensity factor (GSIF), that characterises the singular solution in problems with reentrant corners and cracks, as the quantity of interest.

From [START_REF] Szabó | Finite Element Analysis[END_REF][START_REF] Ródenas | A recovery error estimator for singular problems using singular+smooth field splitting[END_REF] we take the expression to evaluate the GSIF, which reads

Q(u) = K = - 1 C Ω σ jk u aux k -σ aux jk u k ∂q ∂x j dΩ ( 44
)
where u aux , σ aux are the auxiliary fields used to extract the GSIFs in mode I or mode II and C is a constant that is dependent on the geometry and the loading mode. q is an arbitrary C 0 function used to define the extraction zone Ω I which must take the value of 1 at the singular point and 0 at the boundaries and x j refers to the local coordinate system defined at the singularity.

Rearranging terms in the integral in [START_REF] Moumnassi | Finite element analysis on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces[END_REF] we obtain:

Q(u) = K = Ω I (σ) T - 1 C   u aux 1 q ,1 u aux 2 q ,2 u aux 2 q ,1 + u aux 1 q ,2   - (u) T - 1 C σ aux 11 q ,1 + σ aux 21 q ,2 σ aux 12 q ,1 + σ aux 22 q ,2 dΩ (45)
which can be rewritten as a function of initial strains ε0 and body loads b:

Q(u) = K = Ω I σ(u) T ε0 + (u) T bdΩ. ( 46
)
Thus, if we replace u with the vector of arbitrary displacements v, the quantity of interest can be evaluated from

Q(v) = Ω I σ(v) T ε0 dΩ + Ω I v T bdΩ. (47)
Hence, the initial strains and the body loads per unit volume that must be applied in the dual problem to extract the GSIF are defined as

ε0 = - 1 C   u aux 1 q ,1 u aux 2 q ,2 u aux 2 q ,1 + u aux 1 q ,2   , b = 1 C σ aux 11 q ,1 + σ aux 21 q ,2 σ aux 12 q ,1 + σ aux 22 q ,2 (48)

Local contributions

The ZZ error estimate in quantities of interest given in equation ( 31) can be written in terms of local contributions, which proves useful for adaptivity purposes. For a discretisation with n e elements, element e occupying domain Ω e such that Ω e : Ω =

ne Ω e , we can write:

E 1 = ne Ωe e * σ T D -1 ẽ * σ dΩ (49)

Related estimates

In [START_REF] Rüter | Goal-oriented a posteriori error estimates in linear elastic fracture mechanics[END_REF], the authors introduce a hierarchy of practical bounds, by using the triangle and Cauchy-Schwarz inequalities on the exact error measure [START_REF] Zienkiewicz | The Finite Element Method: The Basis[END_REF]. We recall these estimates as they will be used for benchmarking and verification purposes:

E 2 = ne Ωe e * σ T D -1 ẽ * σ dΩ (50)
E 3 = ne Ωe e * σ T D -1 e * σ dΩ Ωe ẽ * T σ D -1 ẽ * σ dΩ (51)
E 4 = Ω e * σ T D -1 e * σ dΩ Ω ẽ * T σ D -1 ẽ * σ dΩ. ( 52
)
The properties of these error estimates are discussed in [START_REF] Rüter | Goal-oriented a posteriori error estimates in linear elastic fracture mechanics[END_REF]. In particular, E 4 is an upper bound if the recovered primal and dual stresses are exactly equilibrated. Such equilibrated recovered stresses were proposed in the context of ZZ estimates in [START_REF] Díez | Equilibrated patch recovery error estimates: simple and accurate upper bounds of the error[END_REF]. Generally speaking, the authors of [START_REF] Rüter | Goal-oriented a posteriori error estimates in linear elastic fracture mechanics[END_REF] observe that the more accurate the estimate, the less likely the practical upper bounding property. In [START_REF] Rüter | Goal-oriented a posteriori error estimates in linear elastic fracture mechanics[END_REF] the authors used the original SPR technique [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique[END_REF]. The following section will show how the use of the SPR-CX improves the accuracy of the error estimator in the quantity of interest both at the local and global level.

Numerical results

In this section 2D benchmark problems with exact solutions are used to investigate the quality of the proposed technique. The first and second problems have a smooth solution whilst the third is a singular problem. For all models we assume plane strain conditions. The h-adaptive FE code used in the numerical examples is based on Cartesian meshes independent of the problem geometry [START_REF] Nadal | On the advantages of the use of cartesan grid solvers in shape optimization problems[END_REF][START_REF] Nadal | Efficient Finite Element methodology based on Cartesian Grids. Application to structural shape optimization[END_REF], with bilinear quadrilateral (Q4) elements3 . To represent the domain geometry accurately, the integrals in elements cut by the boundary are restricted to the part of the element within the domain as in [START_REF] Belytschko | Structured extended finite element methods for solids defined by implicit surfaces[END_REF][START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF]. The integration procedure in these elements is based on the definition of triangular integration subdomains within Ω e and aligning with the geometry of the domain. In elements cut by curved boundaries, these triangular subdomains can have curved boundaries. In that case, the actual geometry is reproduced using a transfinite mapping [START_REF] Gordon | Construction of curvilinear co-ordinate systems and applications to mesh generation[END_REF][START_REF] Sevilla | NURBS-enhanced finite element method ( NE-FEM )[END_REF]. Dirichlet boundary conditions are imposed using Lagrange multipliers following a procedure similar to that described in [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] and, more recently, [START_REF] Moumnassi | Finite element analysis on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces[END_REF] 4 .

To assess the performance of the proposed technique we consider the effectivity index of the error estimator θ defined as the quotient of the estimated error E in the quantity of interest over the exact error Q(e):

θ = E Q(e) . ( 53
)
We can also represent the effectivity in the QoI, θ QoI , defined as the corrected value of the QoI Q(u h ) using the error estimate, divided by the exact value Q(u):

θ QoI = Q(u h ) + E Q(u) . ( 54
)
The relative error in the QoI for the exact and estimated error are

η Q (e) = |Q(e)| |Q(u)| , η Q (e es ) = |E| |Q(u h ) + E| . ( 55
)
The distribution of the local effectivity index D is analysed at the element level, following the definitions described in [START_REF] Ródenas | Improvement of the superconvergent patch recovery technique by the use of constraint equations: the SPR-C technique[END_REF] for the error in the energy norm, adapted here to the case of the error in QoI:

D = θ e -1 if θ e ≥ 1 D = 1 - 1 θ e if θ e < 1 with θ e = E e Q(e e ) , (56)
where superscript e denotes evaluation at the element level. To evaluate Q(e e ) and E e we use [START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF][START_REF] Pled | An enhanced method with local energy minimization for the robust a posteriori construction of equilibrated stress fields in finite element analyses[END_REF] locally on each element. Nonetheless, we should remark that this is only possible for some problems with analytical solutions as the exact value of σ is unknown in the vast majority of cases, especially for the dual problem.

Once the error in the QoI is estimated, the local error estimates in each element can be used to perform h-adaptive analyses using similar techniques to those already available for the error in the energy norm. The refinement of the mesh using the error estimate as the guiding parameter considers a stopping criterion that checks the value of the global estimated error against a prescribed or desired error. If the estimated error is higher than the desired error then the mesh is refined. Several procedures to perform the refinement are available in the literature. To define the size of the elements in the new mesh we follow the adaptive process described in [START_REF] Ladevèze | Accuracy and optimal meshes in finite element computation for nearly incompressible materials[END_REF][START_REF] Coorevits | An automatic procedure with a control of accuracy for finite element analysis in 2D elasticity[END_REF][START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF] which minimises the number of elements in the new mesh for a given accuracy level. This criterion is equivalent to the traditional approach of equally distributing the error in each element of the new mesh as shown in [START_REF] Li | Notes on mesh optimal criteria in adaptive finite element computations[END_REF][START_REF] Fuenmayor | Criteria to achieve nearly optimal meshes in the hadaptive finite element method[END_REF]. These criteria provide the size of the elements in the new mesh as a function of (i ) the ratio of the estimated error in the energy norm in the current mesh to the desired error in the new mesh, and (ii ) the estimated error in the energy norm on each element, which always takes non negative values. They cannot be directly used in goal-oriented adaptivity because the local contributions to the global error in the QoI, evaluated in each element using [START_REF] Pled | An enhanced method with local energy minimization for the robust a posteriori construction of equilibrated stress fields in finite element analyses[END_REF], can take negative values. Thus, for our implementation using h-adaptive routines developed for the error in the energy norm we prepare as input the square root of the absolute value of the error in the QoI -following the relation from the expressions in [START_REF] Díez | A posteriori error estimation for standard finite element analysis[END_REF][START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF] -and the ratio of the estimated error in the QoI in the current mesh to the desired error in the new mesh.

The refinement technique provides a distribution of the required new element sizes. These sizes are specified for each element of the current mesh, which will be recursively split into 4 new elements until the sizes of the elements are smaller than the required size. A maximum difference of only one refinement level will be allowed between adjacent elements. In these cases C 0 continuity will be enforced by means of multipoint constraints [START_REF] Abel | An algorithm for multipoint constraints in finite element analysis[END_REF][START_REF] Farhat | Incorporation of linear multipoint constraints in substructure based iterative solvers. Part 1: a numerically scalable algorithm[END_REF].

Problem 1: Thick-wall cylinder subjected to internal

pressure.

The geometrical model for this problem and the initial mesh are represented in The exact solution for the radial displacement assuming plane strain conditions is given by

u r (r) = P (1 + ν) E(c 2 -1) r (1 -2ν) + b 2 r ( 57
)
where c = b/a, r = x 2 + y 2 and φ = arctan(y/x). Stresses in cylindrical coordinates are

σ r (r) = P c 2 -1 1 - b 2 r 2 σ φ (r) = P c 2 -1 1 + b 2 r 2 σ z (r, φ) = 2ν P c 2 -1 ( 58
)
Figure 2: Thick-wall cylinder subjected to an internal pressure. Model and analytical solution, the domain of interest Ω I is indicated in dark.

Several linear quantities of interest were considered for this problem: the mean radial displacements along Γ o , the mean displacements u x in the domain of interest Ω I and the mean stresses σ x in Ω I .

Problem 1.a.: Mean radial displacements along Γ o

Let Q be the functional that evaluates the mean normal displacement ūn along Γ o such that:

Q(u) = ūn = 1 |Γ o | Γo (Ru) T c u dΓ ( 59
)
where R is the standard rotation matrix for the displacements that aligns the coordinate system with the boundary Γ o and c u = {1, 0} T is the extraction vector that selects the normal component. The exact value of the QoI given by (57) for r = b is ūn = 2.42 6 • 10 -3 .

To characterise the error before using it in the adaptive process we use the set of uniformly refined meshes shown in Figure 3. Figure 4 shows the set of meshes resulting from the h-adaptive process guided by the error estimate in this QoI.

In Tables 1 and2 we show the results for the error estimate E SPR-CX 1 = E 1 from (31), evaluated using the SPR-CX recovery technique in both the primal and the dual problem, and the exact error Q(e) for uniform and h-adapted meshes, respectively. The recovery technique accurately captures the exact error and provides good effectivities, with values of θ = 0.9473 for the coarsest h-adapted mesh with 180 degrees of freedom (DOF). The evolution of the effectivity index for the proposed technique and the standard SPR is represented in Figure 5. In this case, SPR-CX converges faster to the optimal value θ = 1 and shows a high level of accuracy (θ = 0.9928 for 3294 DOF in h-adapted meshes).

Note that in this case, the dual problem corresponds to a traction t = R T c u /|Γ o | that represents a constant traction normal to the external boundary. Therefore, the solution to the dual problem can be evaluated exactly for this quantity of interest using the analytical solution for a cylinder under external pressure:

σ r (r) = - P o b 2 b 2 -a 2 1 - a 2 r 2 σ φ (r) = - P o b 2 b 2 -a 2 1 + a 2 r 2 (60)
where P o represents the applied external pressure. As the exact solution for the dual problem is available, it is possible to evaluate the local effectivity D in ( 56 and effectivities considering the uniform meshes in Figure 3. , for uniform (left) and h-adapted (right) meshes. standard deviation σ(D) of the local effectivity for h-adapted meshes. The ideal scenario is that both parameters are as small as possible and go to zero as we increase the number of DOFs. In the figure we see that the SPR-CX gives a better local estimation, with values of σ(D) and m(|D|) that are smaller than those for the SPR -for the mesh with approx. 3,300 DOF, σ(D) = 0.11 for the SPR-CX, compared to σ(D) = 0.36 for the SPR. As the h-adaptive procedures use local information to refine the mesh, providing an accurate local error estimate to the adaptive algorithm is of critical importance. For this reason, the local performance of the proposed technique indicates that the error estimator based on equilibrated recovered fields for the primal and dual problems is superior to the standard recovery techniques to guide the h-adaptive process, even in cases in which the global effectivity is similar to the effectivity of the standard SPR.

Problem 1.b.: Mean displacements ūx in Ω I

Let us consider the mean displacement ūx in Ω I as the quantity of interest. The objective is to evaluate the error when evaluating ūx defined by the functional:

Q(u) = ūx = 1 |Ω I | Ω I u x dΩ (61)
The exact value of the QoI can be computed for this problem and is ūx = 0.002238239291713. Figure 7 shows the first four meshes used in the h-adaptive refinement process guided by the error estimate for the QoI.

Figure 8 shows the relative error (in percentage) for the error estimates in (31,50-52) evaluated using the proposed recovery technique. The relative error for Q(e) is also shown for comparison. For all the curves the relative error decreases monotonically when increasing the number of DOF, indicating that the h-adaptive process has a stable convergence. The most accurate estimation is obtained for the estimate E 1 , which practically coincides with the exact relative error. The other estimates tend to overestimate the exact error, although strictly speaking they do not have Table 3 shows the estimated error in the QoI, E 1 , the exact error, Q(e), the effectivity in the quantity of interest, θ QoI and the effectivity of the error estimator, gives a good approximation of the exact error. The effectivity of the error estimator θ converges and is very close to the optimal value θ = 1 (with θ = 1.0123 for 75,323 DOF). As expected from these results, the effectivity θ QoI is very accurate as well, with θ QoI = 1.00000012 for 75,323 DOF.

10 0 DOF η Q (%) E 1 E 2 E 3 E 4 Q(e)
If in [START_REF] Pled | An enhanced method with local energy minimization for the robust a posteriori construction of equilibrated stress fields in finite element analyses[END_REF] we consider the case of the non-equilibrated superconvergent patch recovery procedure, resembling the averaging error estimators presented in [START_REF] Rüter | Goal-oriented a posteriori error estimates in linear elastic fracture mechanics[END_REF], we obtain the results shown in Figure 10. This figure shows that, in this case, the effectivity of the error estimator provided by the standard SPR technique is similar to the effectivity obtained with the SPR-CX technique here proposed, although the latter results in more accurate values for the coarsest mesh (θ = 1.1659 for the SPR and θ = 1.0023 for the SPR-CX, considering 528 DOF). In any case, we should recall that the local behaviour with the SPR-CX is generally better than with the SPR.

Problem 1.c.: Mean stress σx in Ω I

Consider as QoI the mean stress value σx given in [START_REF] Giner | A domain integral for the calculation of generalized stress intensity factors in sliding complete contacts[END_REF]. Figure 11 shows the first four meshes of bilinear elements used in the refinement process guided by the error The evolution of the relative error for the estimates presented in (31, 50-52) using the proposed recovery technique and the exact error is shown in Figure 12. Similar to the observations done for the previous examples, the most accurate results are obtained when considering the estimate E 1 . In this case, the other estimates considerably overestimate the true error. Figure 13 shows the effectivity index for E SPR-CX 1 , which uses the locally equilibrated SPR-CX recovery technique, together with the effectivity obtained with the non-equilibrated SPR technique (E SPR 1 curve). This graph clearly shows the improvement obtained with the SPR-CX recovery technique, with effectivities very close to 1 (θ = 0.9797 for the SPR-CX whilst θ = 1.6209 for the SPR, considering 19,573 DOF), in contrast with the oscillatory effectivities provided by the SPR technique.

DOF

η Q (%) Table 4 shows the estimate E SPR-CX 1 , the exact error Q(e), the effectivity index for the QoI θ QoI and the effectivity of the error estimator θ. For this problem the exact value of the QoI is σx = 0.0 6. Table 4 indicates that the equilibrated recovery procedure (SPR-CX) provides very accurate estimations of the error in the QoI and in the value of the QoI itself, with θ = 1.0259 and θ QoI = 1.00000063 for 79,442 DOF. In this test case, in the primal problem the displacements are imposed along the inner boundary of the cylinder (see Figure 14) and the quantity of interest is the mean normal traction along the inner boundary. This is opposite to Problem 1, where tractions were imposed and we were interested in the mean value of radial displacements. In Problem 2, the radial displacements u r (a) to impose are evaluated from [START_REF] Fuenmayor | Criteria to achieve nearly optimal meshes in the hadaptive finite element method[END_REF] such that the exact value for the QoI is tn = 1.

E 1 E 2 E 3 E 4 Q(e)
Table 5 shows the results for the error estimate E SPR-CX 1 using the proposed recovery technique and the exact error Q(e). The procedure accurately captures the error and yields good effectivities for this QoI, with values of θ = 0.9109 for the last mesh with 2,218 DOF. Figure 15 shows the evolution of the effectivity with the number of DOF for the SPR-CX and the SPR. Results show a better performance of the proposed technique when compared with the SPR which presents an oscillatory behaviour (θ = 0.9466 for mesh 2 with 260 DOF and θ = 0.5681 for mesh 3 with 732 DOF). Figure 16 represents the evolution of the mean absolute value m(|D|) and standard deviation σ(D) of the local effectivity. Again, for this example, the SPR-CX gives lower values of these parameters than the SPR. The improved local performance of the SPR-CX is particularly useful for the adaptive algorithm. Notice that for the second mesh, although the global effectivity is close to unity for the SPR, m(|D|) and σ(D) indicate that the SPR-CX is a superior choice when estimating the true error at the element level. Thus, the apparent satisfactory behaviour of the SPR in mesh 2 is due to error compensations of large values of local effectivities.

Problem 3: L-Shape plate

Let us consider the singular problem of a finite portion of an infinite domain with a reentrant corner. The model is loaded on the boundary with the tractions corresponding to the first terms of the asymptotic expansion that describes the exact solution under mixed mode loading conditions around the singular vertex, see Figure 17. The exact values of boundary tractions on the boundaries represented by discontinuous thick lines were imposed in the FE analyses.

The exact displacement and stress fields for this singular elasticity problem can be found in [START_REF] Szabó | Finite Element Analysis[END_REF]. Exact values of the generalised stress intensity factors (GSIF) [START_REF] Szabó | Finite Element Analysis[END_REF] under mixed mode were taken as K I = 1 and K II = 1. The material parameters are Young's modulus E = 1000, and Poisson's ratio ν = 0.3. As the analytical solution of this problem is singular at the reentrant corner of the plate, for the recovery of the dual and primal fields we apply the singular +smooth decomposition of the stresses as explained in Section 3.2. We use a domain integral method based on extraction functions to obtain an approximation of the recovered singular part as explained in Section 3.2.

In this example, we consider the GSIFs K I and K II as the quantities of interest. Figure 18 shows the Cartesian meshes used to solve the primal and dual problems when the mesh is h-adapted for the evaluation of K I . For the dual problem, we use the same Dirichlet conditions as shown in Figure 17 and the set of nodal forces used to extract the QoI in the annular domain Ω I , defined by a plateau function q, shown Table 6: Problem 3. Stress intensity factor K I as QoI using the initial strains and body loads in [START_REF] Sevilla | NURBS-enhanced finite element method ( NE-FEM )[END_REF]. 18. Function q is defined such that q = 1 for r ≤ r 1 = 0.6, q = 0 for r ≥ r 2 = 0.8 and has a smooth transition for r 1 < r < r 2 given by a quartic spline q(s) = 1 -6s 2 + 8s 3 -3s 4 with s = (r -r 1 /r 2 )/(1 -r 1 /r 2 ). In order to impose equilibrium conditions during the recovery of the stresses we use the following approach. For the primal solution, on each patch, we enforce internal equilibrium and compatibility in Ω, and boundary equilibrium all along the Neumann boundary. For the dual problem, we enforce internal equilibrium using the initial strains and body loads given by [START_REF] Sevilla | NURBS-enhanced finite element method ( NE-FEM )[END_REF] and homogeneous Neumann boundary conditions. Compatibility is enforced within all the domain.

Table 6 shows the results for the stress intensity factor K I . Similarly to the results for other QoIs, we observe that the proposed technique provides an accurate estimate E of the exact error Q(e). The effectivity index θ is always close to the optimal value θ = 1 and for the last mesh with 39,193 DOF we obtain θ = 0.9940. As a result, and in agreement with the previous cases, the effectivity in the QoI is highly accurate, with θ QoI = 0.9999 for the same mesh. Table 7 shows the same results for the stress intensity factor K II . Again, we observe a satisfactory behaviour of the error indicator and very accurate effectivities, both for the error estimate and for the QoI itself (θ = 1.0222 and θ QoI = 1.00000053 for the last mesh).

Figure 19 shows the evolution of the relative error with respect to the number of DOF for the two QoIs. Figures 20 and21 show the evolution of the effectivity index of the error estimators obtained with the locally equilibrated SPR-CX technique, and with the non-equilibrated recovery as we increase the number of DOF. The results indicate that the proposed methodology accurately evaluates the error in the QoI, giving values of θ close to 1 and considerably improving the results obtained with Table 7: Problem 3. Stress intensity factor K II as QoI using the initial strains and body loads in [START_REF] Sevilla | NURBS-enhanced finite element method ( NE-FEM )[END_REF]. In this paper, we presented an a posteriori recovery-based strategy that aims to control the error in quantities of interest. The proposed technique is the first attempt to use a recovery procedure that constructs locally equilibrated stress fields for both the primal and the dual problem.

To recover the solution for the primal problem the formulation enforces equilibrium and compatibility and, for singular problems, relies on a singular +smooth stress splitting [START_REF] Ródenas | A recoverytype error estimator for the extended finite element method based on singu-lar+smooth stress field splitting[END_REF]. In order to enforce equilibrium conditions when evaluating the recovered stress fields for the dual problem, we used an approach that allows us to express the functional which defines a given QoIs in terms of body loads, boundary tractions and initial strains and stresses applied to the dual problem. The proposed technique was tested on different quantities of interest: mean displacements and stresses on a domain of interest, mean displacements and tractions along a boundary and the generalised stress intensity factor for a particular singular problem.

The methodology we proposed provides accurate global and local evaluations of the error in the different quantities of interest analysed, improving the results obtained with the original SPR technique, which in some cases can produce satisfactory estimations of the global error due to locally compensating errors. In particular, our approach proves superior to the standard SPR in terms of local error. We also showed that these accurate, local, estimations of the approximation error can be used to drive adaptive procedures.

An extension of the work presented here for extended finite element approximations is now under development. The final goal is to be able to guide the adaptive process in 3D XFEM problems, modelling fatigue crack growth using the stress intensity factors as design parameters for industrial applications of the XFEM, similar to those discussed in the early work of References [START_REF] Bordas | Enriched finite elements and level sets for damage tolerance assessment of complex structures[END_REF][START_REF] Bordas | A simulation-based design paradigm for complex cast components[END_REF][START_REF] Bordas | An extended finite element library[END_REF] and [START_REF] Wyart | A substructured FE-shell/XFE-3D method for crack analysis in thin-walled structures[END_REF][START_REF] Wyart | Substructuring FEXFE approaches applied to three-dimensional crack propagation[END_REF].

Figure 1 :

 1 Figure 1: Equilibrium conditions along internal boundaries.

  hence, a boundary value problem with Dirichlet boundary conditions ũ = -c R /|Γ I | on Γ D .

Figure 2 .

 2 Due to symmetry, only 1/4 of the section is modelled. The domain of interest (DoI) Ω I is denoted by a dark square whereas the contours of interest Γ i and Γ o are the internal and external surfaces of radius a and b. Young's modulus is E = 1000, Poisson's ratio is ν = 0.3, a = 5, b = 20 and the load P = 1.

4 Figure 3 :

 43 Figure 3: Problem 1.a. Sequence of meshes with uniform refinement for the cylinder under internal pressure.

4 4 Figure 4 : 1 Figure 5 : 1 ,

 444151 Figure 4: Problem 1.a. Sequence of h-adapted meshes for the mean radial displacement along Γ o .

1 Figure 6 : 1 ,

 161 Figure 6: Problem 1.a. Evolution of the mean absolute value m(|D|) and standard deviation σ(D) of the local effectivity considering locally equilibrated , E SPR-CX 1 , and non-equilibrated recovery, E SPR 1 .

Figure 7 :

 7 Figure 7: Problem 1.b. Sequence of h-adaptive refined meshes.

Figure 8 :

 8 Figure 8: Problem 1.b. Evolution of the relative error η Q obtained with the SPR-CX technique, considering the error estimates in (31, 50-52) and the exact error Q(e).

4 Figure 9 :

 49 Figure 9: Problem 1.b. Evolution of the effectivity index θ for the error estimates in (31, 50-52) obtained with the SPR-CX technique.

1 Figure 10 : 1 .Figure 11 :

 110111 Figure 10: Problem 1.b. Evolution of the effectivity index θ considering equilibrated E SPR-CX 1

Figure 12 : 1 Figure 13 : 1 ,

 121131 Figure 12: Problem 1.c. Evolution of the relative error η Q considering the error estimates in (31, 50-52) and the exact error Q(e), obtained with the SPR-CX technique.

Figure 14 :

 14 Figure 14: Thick-wall cylinder subjected to an internal pressure.Dirichlet boundary conditions. To formulate the dual problem we consider displacement constraints along Γ I such that u r (r = a) = 1 |Γ I | .(62)

1 Figure 16 :Figure 17 :

 11617 Figure 16: Problem 2. Evolution of the mean absolute value m(|D|) and standard deviation σ(D) of the local effectivity considering locally equilibrated, E SPR-CX 1 , and non-equilibrated recovery, E SPR 1 .

1 ,

 1 101 7.302834 • 10 -3 9.566751 • 10 -3 0.76335570 0.99773608 2,927 2.065838 • 10 -3 2.206343 • 10 -3 0.93631774 0.99985950 10,399 5.194502 • 10 -4 5.280162 • 10 -4 0.98377707 0.99999143 39,193 1.300133 • 10 -4 1.307862 • 10 -4 0.99409058 0.99999923 in Figure

Figure 18 :

 18 Figure 18: Problem 3. Cartesian meshes with h-adaptive refinement.

1 ,

 1 101 2.031449 • 10 -3 1.797716 • 10 -3 1.13001659 1.00023373 2,927 5.273469 • 10 -4 4.342458 • 10 -4 1.21439713 1.00009310 10,399 1.013432 • 10 -4 9.947112 • 10 -5 1.01881993 1.00000187 39,193 2.440128 • 10 -5 2.387141 • 10 -5 1.02219675 1.00000053 the original SPR technique (for the finest mesh considering K II we have θ = 0.8610 for the standard SPR and θ = 1.0222 for the SPR-CX).

KFigure 19 : 1 Figure 20 : 1 Figure 21 :

 19120121 Figure 19: Problem 3. Evolution of the relative error η Q for QoIs K I and K II , for the exact and estimated error as defined in (55), obtained with the SPR-CX technique.

Table 1 :

 1 Problem 1.a. Values for the error estimate E SPR-CX 1

Table 3 :

 3 Problem 1.b. Error estimate E SPR-CX 10 -5 1.666535 • 10 -51.00233710 1.00001740 1,126 1.849599 • 10 -6 1.868635 • 10 -6 0.98981317 0.99999150 4,348 4.240859 • 10 -7 4.046075 • 10 -7 1.04814162 1.00000870 18,076 1.017070 • 10 -7 9.922496 • 10 -8 1.02501376 1.00000111 75,328 2.128196 • 10 -8 2.102331 • 10 -8 1.01230285 1.00000012 θ, for the sharp error estimate E SPR-CX 1 = E 1 . Comparing Q(e) and E 1 we can notice that both values decrease as we refine and that the estimate E SPR-CX 1

	1	and its effectivities.

Table 4 :

 4 Problem 1.c. Values for the error estimate E SPR-CX 10 -5 7.687113 • 10 -5 1.07004501 1.00008077 5,131 2.678276 • 10 -5 3.016986 • 10 -5 0.88773219 0.99994919 19,573 1.108891 • 10 -5 1.131773 • 10 -5 0.97978268 0.99999657 79,442 1.655151 • 10 -6 1.613254 • 10 -6 1.02597049 1.00000063 5.2 Problem 2: Mean normal traction tn along Γ I

	1	and effectivities.

Table 5 :

 5 Problem 2. Values for the error estimate E SPR-CX

	1	and effectivities.

1 Figure 15: Problem 2. Evolution of the effectivity index θ considering locally equilibrated , E SPR-CX 1 , and non-equilibrated recovery, E SPR 1 .

The explanations are restricted to linear quantities of interest. In the developments, affine quantities of the displacement will also be considered, but we will show that this particular case can be recast into the linear case.

The interested reader is referred to the recent paper by Moumnassi and colleagues[START_REF] Moumnassi | Finite element analysis on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces[END_REF] which discusses recent advances in "ambient space finite elements" and proposes hybrid level set/FEMs able to handle sharp corners and edges.

The interested reader may want to refer to[50,

[START_REF] Menk | A robust preconditioning technique for the extended finite element method[END_REF][START_REF] Hiriyur | A quasi-algebraic multigrid approach to fracture problems based on extended finite elements[END_REF][START_REF] Gerstenberger | An algebraic multigrid approach to solve xfem based fracture problems[END_REF] regarding preconditioning techniques for systems sharing similar features.

7 Acknowledgements This work was supported by the EPSRC grant EP/G042705/1 "Increased Reliability for Industrially Relevant Automatic Crack Growth Simulation with the eXtended Finite Element Method". Stéphane Bordas also thanks partial funding for his time provided by the European Research Council Starting Independent Research Grant (ERC Stg grant agreement No. 279578) "RealTCut Towards real time multiscale simulation of cutting in non-linear materials with applications to surgical simulation and computer guided surgery". This work has received partial support from the research project DPI2010-20542 of the Ministerio de Economía y Competitividad (Spain). The financial support of the FPU program (AP2008-01086), the funding from Universitat Politècnica de València and Generalitat Valenciana (PROMETEO/2012/023) are also acknowledged. All authors also thank the partial support of the Framework Programme 7 Initial Training Network Funding under grant number 289361 "Integrating Numerical Simulation and Geometric Design Technology."


Gilles Audemard
email: audemard@cril.fr

Jean-Marie Lagniez
email: lagniez@cril.fr

Bertrand Mazure
email: mazure@cril.fr

Lakhdar Saïs
email: sais@cril.fr

Sathys

SATHYS: Sat Hybrid Solver

Keywords:

à 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.

1 Overview SATHYS (Sat Hybrid Solver) is a new hybridization of local search and modern SAT solver. In our approach, both components heavily cooperate through relevant information gathered during search. More precisely, our hybrid solver alternatively performs the search process using local search and CDCL based SAT solvers. On the one hand, at each node of the search tree, the local search component is used to extend to the model, the current (consistent) partial assignment built by the CDCL based component. On the other hand, the CDCL part is conditionally invoked by the local search component when a local minima is encountered. Each solver benefits from the other in several ways. First, each time a local minima is reached, the local search technique update the activity of the boundary variables [START_REF] Goldberg | Boundary points and resolution[END_REF]. The idea is to direct the CDCL search towards boundary points proven to be important by Goldberg in [START_REF] Goldberg | Boundary points and resolution[END_REF]. Secondly, the polarities of the literals involved in the best complete assignment found during local search are exploited by the CDCL component. From the other side, the CDCL solver shares with local search the current partial assignment together with the learnt clauses. The originality of our proposed hybrid SAT solver arises in alternating search of both components while exchanging relevant information.

The CDCL part includes all the modern enhancements of the DPLL procedure as they can be found in solvers such as RSAT [START_REF] Pipatsrisawat | A lightweight component caching scheme for satisfiability solvers[END_REF] and MINISAT [START_REF] Een | An extensible SAT-solver[END_REF]. These enhancements include watched literal to the unit propagation, first-UIP learning scheme, frequent restarts (Luby strategy [START_REF] Luby | Optimal speedup of las vegas algorithms[END_REF]), activity-based decision heuristics (VSIDS), and the phase learning policy is used [START_REF] Pipatsrisawat | A lightweight component caching scheme for satisfiability solvers[END_REF]. The main improvement consists in storing the binary clauses in an adjacency list, and in reducing the learnt clauses . These learnt clauses are reduced by resolution with the binary clauses. For reducing learnt clauses database, clauses are sorted by using the phase. This step allows to associate a weight with each clause. This weight is the number of falsified literals. Once the clauses are sorted, half of them are kept.

The local search part includes the watched literal speedup [START_REF] Fukunaga | Efficient implementations of sat local search[END_REF] and the metaheuristic used is Rnovelty [START_REF] Mcallester | Evidence for invariants in local search[END_REF].

For not too large instances, we use SatElite as a pre-processor [START_REF] Eén | Effective preprocessing in SAT through variable and clause elimination[END_REF].

The system is written in C and has about 4000 lines of code. It was submitted to the race as a 64 bit binary. It is written from scratch.


Gilles Audemard
email: audemard@cril.fr

Jean-Marie Lagniez
email: lagniez@cril.fr

Bertrand Mazure
email: mazure@cril.fr

Lakhdar Saïs
email: sais@cril.fr

RCL: Reduce learnt clauses

come

Overview

RCL is a SAT solver which includes all the modern enhancements of the DPLL procedure as they can be found in solvers such as RSAT [START_REF] Pipatsrisawat | A lightweight component caching scheme for satisfiability solvers[END_REF] and MINISAT [START_REF] Een | An extensible SAT-solver[END_REF]. These enhancements include watched literal to the unit propagation, first-UIP learning scheme, frequent restarts (Luby strategy [START_REF] Luby | Optimal speedup of las vegas algorithms[END_REF]), activity-based decision heuristics (VSIDS), and the phase learning policy is used [START_REF] Pipatsrisawat | A lightweight component caching scheme for satisfiability solvers[END_REF]. The main improvement consists in storing the binary clauses in an adjacency list, and in reducing the learnt clauses . These learnt clauses are reduced by resolution with the binary clauses. For reducing learnt clauses database, clauses are sorted by using the phase. This step allows to associate a weight with each clause. This weight is the number of falsified literals. Once the clauses are sorted, half of them are kept. For not too large instances, we use SatElite as a pre-processor [START_REF] Eén | Effective preprocessing in SAT through variable and clause elimination[END_REF].

Code

The system is written in C and has about 3000 lines of code. It was submitted to the race as a 64 bit binary. It is written from scratch.