Pierre-Eric Allier

Ludovic Chamoin
email: chamoin@ens-paris-saclay.fr

LMT Pierre Ladevèze

Lmt

Towards simplified and optimized a posteriori error estimation using PGD reduced models

Keywords: constitutive relation error, error estimation, finite element methods, model reduction, proper generalized decomposition

The paper deals with the use of model order reduction within a posteriori error estimation procedures in the context of the finite element method. More specifically, it focuses on the constitutive relation error concept, which has been widely used over the last 40 years for FEM verification of computational mechanics models. A technical key-point when using constitutive relation error is the construction of admissible fields, and we propose here to use the proper generalized decomposition to facilitate this task. In addition to making the implementation into commercial FE software easier, it is shown that the use of proper generalized decomposition enables to optimize the verification procedure and to get both accurate and reasonably expensive upper bounds on the discretization error. Numerical illustrations are presented to assess the performance of the proposed approach.

INTRODUCTION

Numerical models, such as these obtained using the finite element method (FEM), are fundamental in science and engineering activities as they constitute the basic ingredient of simulations that enable to predict the behavior of physical phenomena. Consequently, a permanent issue is the verification of these models, which nowadays can attain very high levels of complexity, to certify the quality of numerical simulations. Model verification deals with the assessment of the numerical (FE) model with respect to an initial mathematical model and implies the estimation of discretization error to control the quality of the approximate numerical solution. It delivers error estimates that can further be used in mesh adaptivity processes. In this context, a large set of a posteriori error estimates has appeared over the last 30 years (see previous studies [START_REF] Verfürth | A Review of a Posteriori Error Estimation and Adaptive Mesh Refinement Techniques[END_REF][START_REF] Babus ˘ka | The Finite Element Method and Its Reliability[END_REF][START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF][START_REF]Verifying Calculations, Forty Years On: An Overview of Classical Verification Techniques for FEM Simulations[END_REF] for an overview). Methods can be roughly classified in 3 categories: (1) those based on field smoothing (initiated in Zienkiewicz and Zhu 5 ); (2) those based on residuals of equilibrium (initiated in Babus ˘ka and Rheinboldt 6 ); and (3) those based on residuals of the constitutive relation (initiated in Ladevèze and Leguillon [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF] ). In this work, we focus on the latter methods that involve the constitutive relation error (CRE) concept.

In the context of Computational Mechanics models in which the constitutive relation is a major component, the CRE concept is a convenient and powerful tool. The idea of CRE is rather simple: so-called admissible fields verifying all equations of the model except the constitutive relation are constructed, then the residual associated with the constitutive relation is measured. The CRE concept was first introduced as a robust a posteriori error estimator in FE computations, [START_REF] Comparaison | De Modèles De Milieux Continus[END_REF] enabling to compute both strict and effective discretization error bounds for linear and more generally convex Structural Mechanics problems and to lead mesh adaptivity processes. It was primarily used for linear thermal and elasticity problems [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF][START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF] before being extended to nonlinear time-dependent problems [START_REF] Ladevèze | A new a posteriori error estimation for nonlinear time-dependent finite element analysis[END_REF][START_REF] Ladevèze | Constitutive relation error estimations for finite element analyses considering (visco)-plasticity and damage[END_REF] and to goal-oriented error estimation. [START_REF] Ladevèze | Strict upper error bounds for calculated outputs of interest in computational structural mechanics[END_REF][START_REF] Ladevèze | Calculation of strict error bounds for finite element approximations of nonlinear pointwise quantities of interest[END_REF][START_REF] Ladevèze | New bounding techniques for goal-oriented error estimation applied to linear problems[END_REF][START_REF] Wang | Computable upper and lower bounds on eigenfrequencies[END_REF] The use of CRE for model verification, for which a general overview can be found in other works, [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF][START_REF] Ladevèze | The Relation Error Method: General Verification Tool. Verifying Calculations, Forty Years on: an Overview of Classical Verification Techniques for FEM Simulations[END_REF] requires in particular the computation of admissible dual fields that are fully equilibrated. This requirement, which is the main practical issue both in terms of computational cost and implementation technicality, was addressed by means of several techniques [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF][START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF][START_REF] Ladevèze | A general method for recovering equilibrating element tractions[END_REF][START_REF] Carstensen | Fully reliable localized error control in the FEM[END_REF][START_REF] Machiels | A flux-free nodal Neumann subproblem approach to output bounds for partial differential equations[END_REF][START_REF] Cottereau | Strict error bounds for linear solid mechanics problems using a subdomain based flux-free method[END_REF][START_REF] Gallimard | A constitutive relation error estimator based on traction-free recovery of the equilibrated stress[END_REF][START_REF] De Almeida | Recovery of equilibrium on star patches using a partition of unity technique[END_REF][START_REF] Parès | Guaranteed energy error bounds for the Poisson equation using a flux-free approach: solving the local problems in subdomains[END_REF][START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF] Ern | A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF][START_REF] Ladevèze | A new non-intrusive technique for the construction of admissible stress fields in model verification[END_REF][START_REF] Pled | On the techniques for constructing admissible stress fields in model verification: performances on engineering examples[END_REF][START_REF] Rey | Study of the strong prolongation equation for the construction of statically admissible stress fields: implementation and optimization[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF] ; we concentrate here on the hybrid-flux approach, referred as Element Equilibration Technique (EET) in the recent literature, that post-processes the primal FE solution at hand and requires local fine-scale computations.

The objective of the paper is to present new numerical tools, on the basis of model reduction techniques and an offline/online strategy that can be coupled to the CRE concept to make this latter fully implementable and exploitable for practical industrial applications. They particularly aim at decreasing the computational cost and technicality level, which are required when computing admissible fields, leading to fast and inexpensive V&V procedures. For that purpose, we decide to refer to the proper generalized decomposition (PGD), which is an a priori model reduction technique that has been extensively used over the last decade to solve multiparametric problems (see previous works [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF][START_REF] Ladevèze | On reduced models in nonlinear solid mechanics[END_REF] ).

First, starting from the classical hybrid-flux approach, we introduce a reduced model at the element level issued from a parameterization of local problems. This approach, later referred as EET-PGD, was primarily introduced in Chamoin et al [START_REF] Chamoin | Synergies between the constitutive relation error concept and PGD model reduction for simplified V&V procedures[END_REF] for a simple case (it is generalized here). It leads to a precomputation of local equilibrium from given data (tractions, element geometry, and material), in an offline step, using PGD. The PGD solution is then used and post-processed in the online step. Second, from the PGD representation mentioned previously, we address optimization of tractions (or hybrid fluxes) by introducing dual unknowns over each element of the mesh. A global optimization problem is thus defined for which all expensive computations (such as numerical integrations) can be performed in the offline step. The associated technique is later referred as HFT-PGD where HFT stands for Hybrid Flux Technique, as equilibration at the element level is indirectly performed here.

In the context of model verification, the CRE concept was already used to control PGD approximations (see a posteriori error estimates developed in other studies [START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF][START_REF] Ladevèze | Toward Guaranteed PGD-reduced Models[END_REF] ) or to directly drive the PGD process with CRE minimization. [START_REF] Allier | Proper generalized decomposition computational methods on a benchmark problem: introducing a new strategy based on constitutive relation error minimization[END_REF] Nevertheless, the use of PGD in CRE implementation has never been investigated and we wish to show here that there are major advantages to do so, in particular for the construction of equilibrated fields. The goal is to give a general framework on the effective use of PGD for model verification with CRE. For the sake of simplicity and clarity, we consider scalar linear elliptic (stationary thermal) problems to introduce technicalities of the proposed approach, but extensions to elasticity, time-dependent, or nonlinear problems are also given with regard to existing literature [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF][START_REF] Ladevèze | A new a posteriori error estimation for nonlinear time-dependent finite element analysis[END_REF] ).

The paper outline is as follows: We review in Section 2 the CRE concept and classical construction of equilibrated fields for a general nonlinear mechanical problem; the use of PGD in addition to CRE for the construction of admissible fields (EET-PGD technique) is shown in detail in Section 3; optimization of tractions using PGD (HFT-PGD technique) is addressed in Section 4; illustrative numerical results are reported in Section 5; eventually, conclusions are drawn in Section 6.

MODEL PROBLEM AND A POSTERIORI ERROR ESTIMATION WITH CRE

Model problem

We consider a structure that occupies a continuum open bounded domain Ω ⊂ R d with boundary 𝜕Ω.Weconsidersmall displacements, quasi-static loadings, and isothermal conditions. The time interval is denoted by [0, T]. The structure, for all t ∈[ 0, T], is placed in an environment characterized by a given displacement u d on 𝜕 1 Ω ⊂𝜕Ω, given tractions F d on 𝜕 2 Ω, and a given body force f d in Ω. 𝜕 1 Ω and 𝜕 2 Ω are complementary parts of 𝜕Ω. The problem that describes the evolution of the structure on [0, T] is to find the displacement-stress pair (u,𝜎 𝜎) that verifies the following:

• kinematic constraints u ∈  [0,T] ; u |𝜕 1 Ω = u d ∀t ∈]0, T[; 𝜖(u)= 1 2 (∇u +∇ T u), (1)
• equilibrium equations (principle of virtual works)

𝜎 𝜎 ∈  [0,T] ; ∫ Ω 𝜎 𝜎 ∶ 𝜖(u * )= ∫ Ω f d u * + ∫ 𝜕 2 Ω F d u * ∀u * ∈  0 , ∀t ∈]0, T[, (2)
• constitutive relations

𝜎𝜎 t = (𝜖( u) |𝜏 , 0 ≤ 𝜏 ≤ t)∀ ( x, t) ∈ Ω×]0, T[. (3)
 [0,T] = L 2 ([0, T]; [H 1 (Ω)] d ) is the space containing the displacement field u defined on Ω×]0, T[.  [0,T] = L 2 ([0, T]; H(div, Ω))
is the space of symmetric tensors containing the stress field 𝜎 𝜎 defined on Ω×]0, T[.  0 is the vector space of specified virtual displacements. The operator  characterizes the mechanical behavior of the material; it may be nonlinear and dependent on history.

Let us denote by  [0,T] ad ⊂  [0,T] the space of displacements verifying the kinematic constraints (1), and  [0,T] ad ⊂  [0,T] the space of stress fields, which are solutions to equilibrium equation 2.

The CRE concept 2.2.1 General idea

The CRE concept is built from a dual approach and measures the residual on the constitutive relation; it is a robust error estimate in the context of model verification.

The notion of error in the constitutive law (CRE) was introduced in Ladevèze and Leguillon. [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF] The general concept is to split the equations of the problem into two groups with different mechanical contents. The first group involves "safe" equations, ie, kinematic and static constraints whereas the second group involves less reliable equations, ie, constitutive relation equations. A solution satisfying the first group is called admissible. An admissible solution is the exact solution to the problem if and only if it satisfies the equations of the second group. Thus, the quality of the admissible solution may be evaluated through the residuals associated to the equations of the second group. Since the second group is related to the constitutive relations, the terminology of error in the constitutive law was naturally introduced. The splitting of the equations into the two groups is not always obvious. For instance, when the material behavior is described using internal variables, the state laws may be seen as "internal" equilibrium relations and thus included in the first group, or seen as constitutive law like the evolution laws and thus included in the second group. The first choice was investigated in Ladevèze and Moës. [START_REF] Ladevèze | A new a posteriori error estimation for nonlinear time-dependent finite element analysis[END_REF] This choice leads to the terminology of dissipation error since the error is measured through the residuals of the sole evolution laws, the state laws being a priori satisfied. An alternative is to consider the state laws as part of the second group; in this case, the error estimate involves the residuals associated to both the state and evolution laws. We shall see that this second choice allows to fit elasticity, plasticity, and viscoplasticity with or without softening in the same framework. The dissipation error is recovered as a particular case by enforcing strongly the state laws.

The CRE definition is based on duality pairing. Constitutive relation error gives guaranteed bounds on a large set of problems. In the following, we give examples of CRE functionals for classical mechanical models.

Linear elasticity

Let us start with the simplest family of problems, ie, elasticity problems. One looks for the final state of the structure at t = T, and the problem does not depend on time. The constitutive relation reads

𝜎 𝜎 = 𝜖(u), ( 4
)
where  is the Hooke tensor, which is symmetric and positive definite. The admissibility spaces are denoted  ad and  ad . An FE solution (u h ,𝜎 𝜎 h ) is computed, where 𝜎𝜎 h = 𝜖(u h ) is equilibrated in the finite element sense. The principle is to associate a new displacement-stress pair denoted ( û,σ 𝜎), which is admissible, ie, in  ad ×  ad .The pair verifies the most reliable equations of the problem; its quality and that of the FE solution is measured through the residual related to the constitutive relation. The crucial point is to build a stress σ 𝜎 that verifies equilibrium equations. Once ( û,σ 𝜎) is obtained, the constitutive relation error reads

E 2 CRE ( û,σ 𝜎)= ∫ Ω [𝜓(𝜖( û)) + 𝜓 * ( σ 𝜎)-σ 𝜎 ∶ 𝜖( û)] = 1 2 || σ 𝜎 -𝜖( û)|| 2  , (5)
where

|| • ||  = √ ∫ Ω •  -1 • is the energy norm. 𝜓(𝜖( û)) = 1 2 𝜖( û)∶𝜖( û) and 𝜓 * ( σ 𝜎)= 1 2  -1 σ 𝜎 ∶ σ
𝜎 are dual (in the Legendre Fenchel sense) convex potentials related to free energy. This error estimate is related to the error defined from the exact solution through the Prager-Synge theorem:

|| û -u|| 2  + || σ 𝜎 -𝜎 𝜎|| 2  = 2E 2 CRE ( û,σ 𝜎). (6)

Viscoplasticity

We consider materials with standard thermodynamics formulation involving internal variables and small displacements without softening such as viscoplasticity. [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF][START_REF] Ladevèze | A new a posteriori error estimation for nonlinear time-dependent finite element analysis[END_REF][START_REF] Ladevèze | Strict upper error bounds for calculated outputs of interest in computational structural mechanics[END_REF] The reference problem is reformulated introducing internal variables (rather than formulating the constitutive relation in functional form). Internal variables other than the plastic strain tensor 𝜖 p (such as the accumulated plastic strain p for isotropic hardening) are gathered in a vector X. ), where Ψ and Ψ * are free energy potentials, and evolution laws ̇ep = B(s).

Λ is assumed to be linear, symmetric, and positive; this is the case of most viscoplastic materials, after a change of variables (normal description). B can be nonlinear and multivalued, as in plasticity; it is assumed to be expressed with two pseudo-potentials of dissipation 𝜑 and 𝜑 * that are dual convex functions (family of standard materials): ̇ep ∈ 𝜕 s 𝜑 * (s) with 𝜕 s denoting subdifferentials; these correspond to the indicatrix function of the elasticity domain for plasticity problems.

The equations of the problem are then divided in two groups:

• The first one related to the free energy defines the admissibility and combines equilibrium equations (s is S-admissible), kinematic constraints ( ̇e = ̇ee + ̇ep is K-admissible), initial conditions (e p = 0att = 0), and state equations; and • The second group, related to dissipation, only includes evolution laws.

The quality of an admissible solution is estimated through the quality of satisfaction of the evolution laws. This defines the CRE error of dissipation type, for which there is a link (similar to the Prager-Synge theorem) with the solution error:

E 2 CRE ( ̇ê p , ŝ)= ∫ T 0 ∫ Ω e 2 CRE ( ̇ê p , ŝ); e 2 CRE ( ̇ê p , ŝ)=𝜑 * (ŝ)+𝜑( ̇ê p )-ŝ • ̇ê p . ( 7
)
Remark 1. A more general framework for the dissipation error (7) was proposed in Ladevèze [START_REF] Ladevèze | Strict upper error bounds for calculated outputs of interest in computational structural mechanics[END_REF] when evolution laws are not given by potentials (but the operator B remains monotonous).

Remark 2. An alternative definition of the CRE functional for nonlinear dissipative models is as follows:

E 2 CRE (ê, ŝ)= ∫ Ω [ Ψ * (ŝ)+Ψ(ê)-ŝ • ê] |T + ∫ T 0 ∫ Ω [ 𝜑 * (ŝ)+𝜑( ̇ê p )-ŝ • ̇ê p ] . ( 8
)
It is based on a definition of admissibility similar to elasticity and involves residuals on both state equations and evolution laws. It is a suitable definition that preserves convexity properties when addressing damage models (softening materials). [START_REF] Ladevèze | Constitutive relation error estimations for finite element analyses considering (visco)-plasticity and damage[END_REF][START_REF] Ladevèze | Constitutive relation error estimators for (visco)plastic finite element analysis with softening[END_REF]

Recovering admissible fields 2.3.1 Elasticity case

We usually choose û = u h . The construction of σ 𝜎 is more technical, and several methods have been proposed. The optimal method would be a full dual analysis (minimizing the complementary energy using equilibrium elements, see previous studies [START_REF] De Veubeke | Dual analysis for heat conduction problems by finite elements[END_REF][START_REF] Oden | On dual complementary variational principles in mathematical physics[END_REF][START_REF] De Veubeke | Displacement and equilibrium models in the finite element method[END_REF][START_REF] De Almeida | Equilibrium Finite Element Formulations[END_REF] ). However, this is hardly applicable in practice, as it requires the solution of an additional global problem, with substantial computational efforts and nonconventional FE spaces that are usually not available in commercial software.

We focus on the hybrid-flux method (or EET) that enables to compute an admissible field σ 𝜎 h ∈  ad from a post-processing of the FE field 𝜎 𝜎 h at hand, with local independent computations. [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF][START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF][START_REF] Ladevèze | A general method for recovering equilibrating element tractions[END_REF] T h ek e yi st h ep r o l o n g a t i o n condition:

∫ K ( σ 𝜎 h -𝜎 𝜎 h ) 𝛁 𝜙 i = 0, (9)
which is enforced for all elements K and all nodes i connected to K. 𝜙 i is the scalar shape function associated with node i. The recovery method can be split in two steps:

Step 1: recovery of equilibrated tractions FK|Γ defined on edges (or faces in 3D) Γ of each element K,with FK|Γ = F d if Γ ⊂𝜕 2 Ω, such that equilibrium at the element level is verified:

∫ K f d • u * R + ∫ 𝜕K FK • u * R = 0 ∀u * R ∈  R (K), (10)
where  R (K) denotes the space of rigid body motions on K. The prolongation condition (9) then yields (using

σ 𝜎 h|𝜕K • n K = FK and 𝛁 • σ 𝜎 + f d = 0): ∫ 𝜕K FK 𝜙 i = ∫ K (𝜎 𝜎 h 𝛁 𝜙 i -f d 𝜙 i ), (11)
and leads to small well-posed systems of equations (over patches of elements connected to each node i), which can be solved explicitly. In practice, tractions are defined as FK|Γ = 𝜂 Γ K FΓ ,with𝜂 Γ K = ±1, to ensure continuity across element boundaries and are searched as a linear combination of FE shape functions:

FK|Γ (x)= ∑ j∈ Γ F j K|Γ 𝜙 j (x)
.  Γ denotes the set of nodes connected to the edge/face Γ. All technical details on the construction of equilibrated tractions can be found in other works [START_REF] Ladevèze | Mastering Calculations in Linear and Nonlinear Mechanics[END_REF][START_REF] Pled | On the techniques for constructing admissible stress fields in model verification: performances on engineering examples[END_REF] ;

Step 2: local recovery of an equilibrated stress field σ 𝜎 h ,forgiventractions FK , over each element K solving the Neumann problem:

∫ K σ 𝜎 h ∶ 𝜖(u * )= ∫ K f d • u * + ∫ 𝜕K FK • u * ∀u * ∈  (K), (12)
with 12) may be obtained analytically, using polynomial functions with sufficiently high degree, provided the source term f d is polynomial as well. [START_REF] Ladevèze | New advances on a posteriori error on constitutive relation in finite element analysis[END_REF] In practice, an alternative approach with numerical solution and higher-order elements is preferred. Indeed, the optimal stress field verifying (12) is the one that minimizes over K the local error estimate || σ 𝜎 -𝜎 𝜎 h || ,K (or equivalently || σ 𝜎|| ,K ). Duality arguments show that this is equivalent to taking σ 𝜎 h|K = 𝜖(𝝆),with𝝆 ∈  (K) verifying:

 (K)=[ H 1 (K)] d . A solution to (
∫ K 𝜖(𝝆)∶𝜖(u * )= ∫ K f d • u * + ∫ 𝜕K FK • u * ∀u * ∈  (K). (13)
A numerical approximation of the solution of (13) (defined up to rigid body motions in  R (K))canbeobtained using the FEM with a single finite element of degree p+k,wherep denotes the polynomial degree used to compute u h and k denotes the extra degree. Numerical studies performed in Babus ˘ka et al [START_REF] Babus ˘ka | Validation of a posteriori error estimators by numerical approach[END_REF] showed that analytical and numerical approaches give similar CRE error estimates choosing k ≥ 3, even though the flux field is not rigorously equilibrated in each element K with the latter approach. We consider the numerical approach in the remainder of the paper.

Solving ( 13) is in practice the most costly part in the hybrid-flux method (in particular for 3D applications), as it involves high-order elements and has to be performed for each element K. We wish to use model order reduction to circumvent this issue.

Time-dependent nonlinear case

The FE solution is computed at time points t n , considering the classical incremental FEM. Assuming that data are piecewise linear on [0, T], we complete the FE solution on [0, T] and get (u h ,𝜎 𝜎 h ). The stress field 𝜎 𝜎 h is FE equilibrated on [0, T]. We usually choose ûh = u h (except when incompressibility needs to be enforced [START_REF] Ladevèze | Accuracy and optimal meshes in finite element computation for nearly incompressible materials[END_REF] ).

The crucial point is to build an admissible stress field. For that purpose, the classical recovery used for elasticity is still used for any time point t n , then a piecewise linear interpolation is performed. Other internal variables, which are not involved in equilibrium equations, are obtained solving local and explicit minimization problems.

ACCELERATION OF LOCAL RECOVERY USING PGD

Scalar problem

For the sake of simplicity, we first consider a scalar problem, ie, a steady-state thermal problem that consists in finding the temperature/flux pair (u, q) such that

u ∈  ; u = u d on 𝜕 1 Ω (kinematic constraints) q ∈  ; ∫ Ω q • 𝛁 u * = ∫ Ω f d u * + ∫ 𝜕 2 Ω F d u * ∀u * ∈  0 (balance equations) q =  𝛁 u (constitutive relation) (14)
with

 = H 1 (Ω) and  = H(div, Ω) = {𝝅 ∈[ L 2 (Ω)] d , 𝛁 •𝝅 ∈ L 2 (Ω)}.
 is a symmetric positive definite operator (second-order tensor). We introduce the energy semi-norm

|| • ||  = √ ∫ Ω 𝛁(•) •  𝛁(•).
For an admissible pair ( û, q)∈  ad ×  ad , the CRE functional is defined as follows:

E 2 CRE ( û, q)= 1 2 ∫ Ω ( q - 𝛁 û)• -1 ( q - 𝛁 û)= 1 2 || q - 𝛁 û|| 2  . ( 15
)
It measures the nonverification of the constitutive relation for ( û, q) and can be used as an error estimate. Indeed, after computing an FE solution (u h , q h ), choosing û = u h and constructing qh ∈  ad lead to the Prager-Synge theorem ( 6) and to the following:

|| u h -u||  ≤ √ 2E CRE (u h , qh ). (16)
The quality of the estimate depends on that of qh . Classically, qh is constructed using the hybrid-flux approach as described in Section 2.3, with the following two steps and specifying that tractions refer here to normal thermal fluxes on element boundaries:

Step 1: construction of equilibrated tractions FK|Γ = 𝜂 Γ K FΓ on edges/faces Γ of each element K, such that equilibration at the element level is verified:

∫ K f d + ∫ 𝜕K FK = 0. ( 17
)
This is performed using the prolongation condition, and tractions are searched as linear combinations of FE shape functions: 1). Step 2: local construction, for given tractions FK and over each element K,of qh|K =  𝛁 𝜌 solving the following local Neumann problem with higher-order elements:

FK|Γ (x)= ∑ j∈ Γ F j K|Γ 𝜙 j (x) (see Figure
∫ K  𝛁 𝜌 • 𝛁 u * = ∫ K f d u * + ∫ 𝜕K FK u * ∀u * ∈  (K). (18)
Solving ( 18) is in practice the most costly part in the hybrid-flux method (in particular for 3D applications), as it involves high-order elements and has to be performed for each element K. We wish to use the PGD technique to find, in an offline phase, a parameterized solution to (18), valid for any configuration of the geometry and the loading. From equilibrated tractions computed in step 1 and with respect to problem data (material parameters, mesh geometry, … ), this PGD solution would then be directly used in the online error estimation phase, and in a multiquery context, for each element K of the mesh.

Parametrization of the local problems

In the following, we consider that the material behavior is isotropic and that material parameters are constant over each element K, so that their values have no influence on qh|K ;wethusset = I when solving (18) and define qh|K = 𝛁 𝜌.In cases where  is not constant over each element or is not isotropic, its structure could be parameterized and additional material parameters could be introduced in the PGD decomposition.

On each edge/face Γ of any element K, tractions are linear combinations of FE shape functions and thus read as follows:

FK|Γ (x)= ∑ j∈ Γ F j K|Γ 𝜙 j (x); F j K|Γ = 𝜂 Γ K F j Γ ∈ R. (19)
Consequently, the solution 𝜌 to (18) can be written as a linear combination of elementary solutions:

𝜌(x)= ∑ Γ⊂𝜕K ∑ j∈ Γ F j K|Γ 𝜌 j Γ (x), (20)
where 𝜌 j Γ is the solution (up to a constant) to the elementary problem:

∫ K 𝛁 𝜌 j Γ • 𝛁 u * = ∫ Γ 𝜙 j u * -∫ K ( 1 |K| ∫ Γ 𝜙 j ) u * ∀u * ∈  (K). (21)
The solution 𝜌 j Γ to each problem ( 21) can be computed with the PGD technique, for any element K, parameterizing the geometry of K with a set of parameters p geo ∈ . Following the approach described in other studies, [START_REF] Ammar | Parametric solutions involving geometry: a step towards efficient shape optimization[END_REF][START_REF] Zlotnik | Proper generalized decomposition of a geometrically parametrized heat problem with geophysical applications[END_REF][START_REF] Courard | Integration of PGD-virtual charts into an engineering design process[END_REF] we reformulate the weak problem ( 21) by introducing a parameter-dependent mapping (p geo )∶K ref → K(p geo ) from a reference fixed element K ref to the geometrically parameterized element K(p geo ). Such a geometrical transformation then allows defining the weak problem in a tensor product space and applying the PGD method, to compute generic parameterized solutions 𝜌 j Γ (p geo ) that can be used for any element geometry. Remark 3. In the presence of geometrical variabilities, an alternative approach described in other works [START_REF] Canuto | A fictitious domain approach to the numerical solution of PDEs in stochastic domains[END_REF][START_REF] Nouy | Fictitious domain method and separated representations for the solution of boundary value problems on uncertain parameterized domains[END_REF] could also be used. It consists in embedding the parameterized domain into a fixed fictitious domain.

In the present case, the mapping is defined from two transformations:

• a first scaling mapping  1 ∶ K → K maps a homothetic element K with diameter 1 to the actual element K with diameter 𝛼. This mapping involves a transformation matrix T 1 = 𝛼I d such that x = T 1 x; • a second mapping  2 ∶ K ref → K maps a reference element K ref to element K. This mapping uses an isoparametric formulation and involves a transformation matrix T 2 such that x = T 2 x ref .

The global mapping is therefore (p geo )= 1 •  2 with transformation matrix T(p geo )=𝛼T 2 , Jacobian matrix J = T, and Jacobian J(p geo )=det(J).

As an illustration, we consider two element types. We first consider 3-node triangle elements (Figure 2). The mappings read

( x ȳ ) = ( ∑ 3 i=1 xi 𝜙 i (𝜂, 𝜉) ∑ 3 i=1 ȳi 𝜙 i (𝜂, 𝜉) ) = T 2 ( 𝜂 𝜉 ) ; T 2 = [ 1 x3 0 ȳ3 ] , (22)
where (x 3 , ȳ3 ) are local coordinates of node 3 in the coordinate system associated with element K and (𝜂, 𝜉) are local coordinates in the coordinate system associated with element K ref .

The global mapping  thus involves 3 parameters, ie, p geo =(𝛼, x3 , ȳ3 ) and J = 𝛼 2 ȳ3 . We now consider 4-node tetrahedron elements (Figure 3). The mappings read

( x ȳ z ) = ⎛ ⎜ ⎜ ⎝ ∑ 4 i=1 xi 𝜙 i (𝜂, 𝜉, 𝜏) ∑ 4 i=1 ȳi 𝜙 i (𝜂, 𝜉, 𝜏) ∑ 4 i=1 zi 𝜙 j (𝜂, 𝜉, 𝜏) ⎞ ⎟ ⎟ ⎠ = T 2 ( 𝜂 𝜉 𝜏 ) ; T 2 = [ 1 x3 x4 0 ȳ3 ȳ4 00z 4 ] , (23)
where

(
∫ K ref J.J -T 𝛁 𝜌 j Γ • J -T 𝛁 u * = ∫ Γ ref J s 𝜙 j u * -∫ K ref ( 1 |K| ∫ Γ ref J s 𝜙 j ) u * ∀u * ∈  (K ref ), (24) 3 3 3
K K ref 1 2 1 1 2 2 1 1 1 y 1 2

FIGURE 2

Mapping between reference and parameterized elements for a 3-node triangle element with J s the Jacobian associated with the transformation of element edge/face Γ. For 3-node triangle elements (Figure 2), it reads

J s = 𝛼 ⎧ ⎪ ⎨ ⎪ ⎩ 1o n Γ 12 √ (1 -x3 ) 2 + ȳ2 3 on Γ 23 √ x2 3 + ȳ2 3 on Γ 13 , ( 25
)
whereas for 4-node tetrahedron elements (Figure 3), it reads

J s = 𝛼 2 ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ȳ3 on Γ 123 √ ȳ2 4 + z2 4 on Γ 124 √ z2 4 (x 2 3 + ȳ2 3 )+(x 4 ȳ3 -x3 ȳ4 ) 2 on Γ 134 √ z2 4 ((x 3 -1) 2 + ȳ2 3 )+((x 4 -1)ȳ 3 -(x 3 -1)ȳ 4 ) 2 on Γ 234 . ( 26
)
Remark 4. Equivalently, ( 24) can be recast as follows:

∫ K ref 𝛁 𝜌 j Γ • K 𝛁 u * = ∫ Γ ref Fd u * + ∫ K ref fd u * ∀u * ∈  (K ref ), (27)
with K = J.J -1 J -T , Fd = J s 𝜙 j ,and fd =-J.

(

|K| ∫ Γ ref J s 𝜙 j 1
) .

Remark 5. The number of elementary problems (21) (or (24)) and the number of parameters involved in the mapping  (size of p geo ) depend on the type of FE element. For instance,

• 3-node triangle elements involve 6 elementary problems (2 for each of the 3 edges) and 3 geometrical parameters (6 degrees of freedom with 3 rigid body motions); • 6-node triangle elements involve 9 elementary problems (3 for each of the 3 edges) and 9 geometrical parameters (12 degrees of freedom with 3 rigid body motions); • 4-node tetrahedron elements involve 12 elementary problems (3 for each of the 4 faces) and 6 geometrical parameters (12 degrees of freedom with 6 rigid body motions).

Solution using PGD

We now introduce a reduced model, based on PGD, at the element level to facilitate step 2 of the hybrid-flux approach.

For that purpose, approximations of solutions 𝜌 j Γ (x ref , p geo ) are computed offline and once for all using the PGD technique with separated-variable modal decomposition.

Considering a general linear D-dimensional problem defined on the tensor space

 =  1 ⊗  2 ⊗ •••⊗  D ,
PGD is a low-rank tensor method that consists in searching an approximation u m of u in a low-dimensional tensor subspace of  m ⊂  made of canonical format tensors of rank m:

u m = m ∑ i=1 w 1 i ⊗ w 2 i •••⊗ w D i , w 𝜇 i ∈  𝜇 . ( 28
)
Among the various strategies to construct u m , [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time dependent partial differential equations[END_REF] we focus on the one called progressive Galerkin. Introducing a global weak formulation of the problem,

find u ∈  such that B(u, v)=F(v)∀ v ∈ , (29)
and assuming that the rank m -1 decomposition u m-1 is known, the rank

m decomposition u m = u m-1 + w 1 ⊗ w 2 •••⊗ w D is searched such that B(u m ,𝛿v)=F(𝛿v)∀ 𝛿v = 𝛿w 1 ⊗ w 2 •••⊗ w D + w 1 ⊗𝛿w 2 •••⊗ w D +•••+w 1 ⊗ w 2 •••⊗𝛿w D , ( 30
)
with 𝛿w 𝜇 ∈  𝜇 . This formulation naturally leads to a nonlinear problem where a set of coupled low-dimensional problems have to be solved:

B(w 1 ⊗ w 2 •••⊗ w D ,𝛿w 1 ⊗ w 2 •••⊗ w D )=R m-1 (𝛿w 1 ⊗ w 2 •••⊗ w D )∀ 𝛿w 1 ∈  1 , B(w 1 ⊗ w 2 •••⊗ w D , w 1 ⊗𝛿w 2 •••⊗ w D )=R m-1 (w 1 ⊗𝛿w 2 •••⊗ w D )∀ 𝛿w 2 ∈  2 , ⋮ =⋮ B(w 1 ⊗ w 2 •••⊗ w D , w 1 ⊗ w 2 •••⊗𝛿w D )=R m-1 (w 1 ⊗ w 2 •••⊗𝛿w D )∀ 𝛿w D ∈  D , (31) with R m-1 (v)=F(v)-B(u m-1 , v).
This problem is in practice solved using alternated minimization with an iterative (fixed point) strategy.

To solve (24) with the progressive Galerkin approach of PGD, we introduce the following bilinear form B and linear form F constructed from the parameterized separated variable Jacobian transformation (all technical details can be found in other works [START_REF] Ammar | Parametric solutions involving geometry: a step towards efficient shape optimization[END_REF][START_REF] Zlotnik | Proper generalized decomposition of a geometrically parametrized heat problem with geophysical applications[END_REF] ). Introducing the interval I p k geo in which parameter p k geo evolves (1 ≤ k ≤ P = size(p geo )), these forms read

B(𝜌 j Γ , v)= ∫ I 𝛼 ∫ I p geo,1 … ∫ I p geo,P ∫ K ref J.J -T 𝛁 𝜌 j Γ • J -T 𝛁 v, F(v)= ∫ I 𝛼 ∫ I p geo,1 … ∫ I p geo,P [ ∫ Γ ref J s 𝜙 j v -∫ K ref J.
(

|K| ∫ Γ ref J s 𝜙 j ) v ] . ( 1
) 32
The parameter 𝛼 is included only for completeness of the description; it acts as a multiplicative constant in the solution 𝜌 and disappears when computing qh|K =  𝛁 𝜌. PGD solutions thus read

𝜌 j Γ,m (x ref , p geo )=𝛼 m ∑ i=1 ( 𝜓 i (x ref ) ∏ 2≤k≤P 𝛿 p geo,k i (p geo,k ) ) . ( 33
)
Space functions 𝜓 i (x ref ) in ( 33) are computed using the FEM with a single element of degree p + k. Other functions are discretized using a fine grid over individual spaces I p geo,k . Additionally, square roots involved in J s are approximated using SVD and polynomial interpolation with high degree (degree 8 in practice). Remark 6. The number m of PGD modes that is required to get accurate solutions 𝜌 j Γ,m can be rigorously defined using classical a posteriori error estimation tools devoted to PGD. [START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF][START_REF] Ladevèze | Toward Guaranteed PGD-reduced Models[END_REF][START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] De Almeida | A basis for bounding the errors of proper generalised decomposition solutions in solid mechanics[END_REF] For extremal cases (elements with bad shape), the number of PGD modes should be chosen carefully; however, too much distorted elements are usually avoided in practice (regular meshes) and a few number of PGD modes leads to a PGD error, which is negligible compared to the discretization error that we wish to assess. A numerical evaluation of the value m that yields sufficient accuracy is provided in the numerical results.

Combining (20) and (33), a PGD representation of 𝜌(x) solution to (18) reads

𝜌 m (x ref , p geo , { F j K|Γ }) = ∑ Γ⊂𝜕K ∑ j∈ Γ F j K|Γ 𝜌 j Γ,m (x ref , p geo ). (34)
It provides for a parameterized equilibrated flux field at the element level:

qh,m|K (x ref , p geo , { F j K|Γ }) = J -T (p geo ) 𝛁 ( ∑ Γ⊂𝜕K ∑ j∈ Γ F j K|Γ 𝜌 j Γ,m (x ref , p geo ) ) , ( 35
)
which can be directly used online in the a posteriori error estimation procedure.

A representation of the use of PGD in step 2 of the hybrid-flux approach (so-called EET-PGD technique) is given in Figure 4.

Remark 7. In practice, the numerical implementation of the EET-PGD technique is organized so that the loop over all elements of the mesh is avoided. It is rather chosen to store all PGD space modes in a single vector and to perform global matrix multiplications which are optimized in classical software.

Extension to elasticity problems

Considering elasticity problems, the local Neumann problems to solve in step 2 read

∫ K 𝜖(𝝆)∶𝜖(u * )= ∫ K f d • u * + ∫ 𝜕K FK • u * ∀u * ∈  (K), ( 36
)
where  is the fourth-order symmetric elasticity tensor. On each edge/face Γ of any element K, tractions are linear combinations of FE shape functions and thus read (in the 3D case)

FK|Γ (x)= ∑ j∈ Γ [ F jx K|Γ 𝝓 x j (x)+ F jy K|Γ 𝝓 y j (x)+ F jz K|Γ 𝝓 z j (x)] ; ( F jx K|Γ , F jy K|Γ , F jz K|Γ ) T = 𝜂 Γ K ( F jx Γ , F jy Γ , F jz Γ ) T ∈ R 3 , ( 37
)
with 𝝓 x j (x)=( 𝜙 j (x), 0, 0) T , 𝝓 y j (x)=( 0,𝜙 j (x), 0) T ,and𝝓 z j (x)=( 0, 0,𝜙 j (x)) T . Consequently, the solution 𝝆 to (36) can be written as follows:

𝝆(x)= ∑ Γ⊂𝜕K ∑ j∈ Γ [ F jx K|Γ 𝝆 jx Γ (x)+ F jy K|Γ 𝝆 jy Γ (x)+ F jz K|Γ 𝝆 jz Γ (x)], (38)
where

𝝆 jx∕y∕z Γ
is the solution, up to a rigid body motion, to the elementary problem (generalization of ( 21)):

∫ K 𝜖(𝝆 jx∕y∕z Γ )∶𝜖(u * )= ∫ Γ 𝝓 x∕y∕z j • u * -∫ K (a 1 ∧ X + a 2 )•u * ∀u * ∈  (K) (39)
X are barycentric coordinates in element K,anda 1 and a 2 are defined as follows:

a 1 = ∫ Γ X ∧ 𝝓 x∕y∕z j ∫ K X • X ; a 2 = 1 |K| ∫ Γ 𝝓 x∕y∕z j . ( 40
)
As for the scalar case, a PGD representation of 𝝆 is then constructed as follows:

𝝆 m (x ref , p geo , { F jx∕y∕z K|Γ }) = ∑ Γ⊂𝜕K ∑ j∈ Γ [ F jx K|Γ 𝝆 jx Γ,m (x ref , p geo )+ F jy K|Γ 𝝆 jy Γ,m (x ref , p geo )+ F jz K|Γ 𝝆 jz Γ,m (x ref , p geo )] (41)
and a parameterized equilibrated stress field at the element level reads

σ 𝜎 h,m|K (x ref , p geo , { F jx∕y∕z K|Γ }) = (J -T (p geo )∇) s 𝝆 m (x ref , p geo , { F jx∕y∕z K|Γ }). (42)

OPTIMIZATION OF TRACTIONS

The quality of the CRE estimate is related to that of constructed equilibrated tractions, so that optimizing tractions FK|Γ is a relevant task. However, this leads to a global optimization problem that is complex to solve without sophisticated numerical methods. In previous studies, [START_REF] Ladevèze | New advances on a posteriori error on constitutive relation in finite element analysis[END_REF][START_REF] Florentin | Evaluation of the local quality of stresses in 3d finite element analysis[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] a partial optimization was investigated. The basic idea was to give greater freedom in the construction of equilibrated tractions, applying the prolongation condition on nonvertex nodes alone, to improve the quality of the recovered admissible flux fields. To keep a reasonable computational cost, such an enhanced construction of equilibrated tractions may in practice be applied only in zones of high element aspect ratios or sharp gradients, in which the quality of admissible flux fields may be affected with usual procedures.

In this section, we wish to take advantage of the PGD representation (35) (which provides for a reduced model at the element level) to address a full optimization of tractions FK|Γ by minimizing the CRE functional. This procedure, which would be out of reach with a classical approach, is highly facilitated by the explicit dependency of the admissible field qh,m on { F j K|Γ } provided by the PGD. It is referred as the HFT-PGD technique.

Optimization problem

We define the (infinite dimensional) space of equilibrated tractions:

 = { { FΓ (x)}; ∫ K f d + ∑ Γ⊂𝜕K ∫ Γ 𝜂 Γ K FΓ = 0f o r a l lK,𝜂 Γ K FΓ = F d if Γ ⊂𝜕 2 Ω } . ( 43
)
In this space, we consider the (finite dimensional) subspace of tractions that have a given form on each edge/face (linear combination of shape functions):

 𝓁 = { { FΓ (x)} ∈  ; FΓ (x)= ∑ j∈ Γ F j Γ 𝜙 j (x) } ⊂  . ( 44
)
The dimension of  𝓁 depends on the polynomial degree of FE shape functions and the number of edges/faces. Consequently, the unknowns of the problem are interface data { F j Γ }, which are gathered in a global vector F of length N F .This way, the constraint that { F j Γ } provide tractions that belong to  can be written C F = Fd . Consequently, the optimization problem reads

{ FΓ (x)} opt = arg min { FΓ }∈ 𝓁 E 2 CRE u h , qh ({ FΓ }) or Fopt = arg min F,C F= Fd E 2 CRE ( u h , qh ( F) ) . ( 45
)

Solution using PGD

On the one hand, from the PGD representation (35) of qh,m|K ,

qh,m|K (x ref , p geo , { F j Γ }) = J -T (p geo ) 𝛁 ( ∑ Γ⊂𝜕K ∑ j∈ Γ 𝜂 Γ K F j Γ 𝜌 j Γ,m (x ref , p geo ) ) , ( 46
)
the global admissible field qh,m computed with PGD can be written as follows:

qh,m = N(x ref , p geo ) F, ( 47
)
where N(x ref , p geo ) is a parameterized d × N F matrix. On the other hand, the PGD representation of

E 2 CRE (u h , qh,m ({ F j Γ })) reads E 2 CRE (u h , qh,m )= 1 2 ∫ Ω ( qh,m - 𝛁 u h )• -1 ( qh,m - 𝛁 u h ) = 1 2 || qh,m || 2  -∫ Ω qh,m 𝛁 u h + 1 2 || u h || 2  = 1 2 FT M F - FT LU h + 1 2 || u h || 2  ( 48
)
with

M = ∫ Ω N T N and L = ∫ Ω N T B (BU h = 𝛁 u h ).
Then the optimization problem is recast as follows:

Fopt = arg min F,C F= Fd ( 1 2 FT M F - FT LU h ) . ( 49
)
It is solved using a Lagrangian functional that leads to the system:

( MC T C 0 )( Fopt Λ ) = ( LU h Fd ) . ( 50
)
We notice that elementary matrices that compose M, C,andL can be precomputed in the offline phase and defined under a PGD form with separated variables decomposition, from both matrices defined at the element level from (46) and the topology of the mesh. Technically, coefficients of elementary matrices are computed and stored once for all with a parameterized PGD form. From the topology of the mesh and the shape of elements, these coefficients are particularized by means of evaluation of PGD terms and the assembly of global matrices is performed under a sparse format without resorting to any numerical integration in the online phase. The global system (50) of size (N F +N F d )×(N F +N F d ), N F d being the number of constraints imposed in  , is then factorized and solved in this latter phase. The number of PGD modes that are used in this process is defined from a posteriori error estimation at the element level as suggested in Section 3.3.

Once the system ( 50) is inverted, the CRE estimate is directly recovered as follows:

E 2 CRE = 1 2 FT opt M Fopt - FT opt LU h + 1 2 || u h || 2  . ( 51
)
Remark 8. When considering nonlinear problems, the CRE functional may be nonlinear with respect to qh,m so that iterative minimization methods are required. An alternative is to minimize || qh,mq h ||  .

A representation of the obtained HFT-PGD technique, that uses PGD model reduction all along the CRE process and thus leads to a convenient alternative to the classical hybrid-flux approach, is given in Figure 5.

NUMERICAL RESULTS

In this section, we evaluate performance of the proposed approach on two-dimensional (2D) and three-dimensional (3D) numerical experiments. In all studied cases, we consider regular meshes such as these classically given by mesh generators (no badly shaped or distorted elements). Consequently, the geometrical parameters considered in the PGD take values that remain in a convenient interval for which PGD decompositions lead to a negligible truncation error with a reasonable number of modes. In addition, the numerical experiments involve an FE approximation with first-order Lagrange basis functions.

2D thermal problem over a regular physical domain

We perform discretization error estimation on a 2D holed plate Ω, which is discretized by means of a mesh composed of 3-node triangle elements (Figure 6). The plate has length 2H and height H with H = 10, and the circular hole is centered on the plate with radius r = 1. We consider a steady-state thermal problem and homogeneous isotropic material properties with  = I. A prescribed zero temperature is applied on the external boundary, while a normal flux F d = 1 is imposed on the inner boundary, ie, on the hole boundary. Because of symmetries, only one quarter of the plate is studied.

From the associated FE solution, equilibrated tractions may be computed using the first step of the hybrid-flux (or EET) technique.

Details on the PGD solution

Using notations of Figure 2, we compute parameterized solutions 21) with a single fourth-order FE element and 20 PGD modes (m = 20). The domains I x3 =[ 0, 1] and I ȳ3 =[ 0.1, 1] are discretized with 100 points each, after checking that this is sufficient to ensure an accurate description of the evolutions with respect to x3 and ȳ3 . The first 3 PGD modes of 𝜌 1 Γ 12 ,20 are shown in Figure 7. In Figure 8, we represent the PGD approximation of 𝜌 1 Γ 12 for various configurations of parameters x3 and ȳ3 . The computation of this PGD solution is done once for all, in an offline phase and stored for later use.

𝜌 j Γ kl ,m (x ref ,𝛼,x 3 , ȳ3 )=𝛼 ∑ m i=1 𝜓 i (x ref )𝛿 x i (x 3 )𝛿 y i (ȳ 3 ) of (

FIGURE 6

Representation of the 2D domain, associated FE mesh, and applied loading After identifying the PGD parameters 𝛼, x3 ,a n dȳ 3 over each element of the mesh (see Figure 9), an accurate PGD approximation of the admissible flux qm can then be directly evaluated inside each element in an inexpensive online phase; we recall that this method is referred as EET-PGD method.

CRE estimate obtained from EET-PGD

From equilibrated tractions and PGD solutions, we have all ingredients to estimate the discretization error using the CRE method. In Figure 10, we compare local contributions to the CRE estimate 2E 2 CRE , obtained from the EET-PGD technique when computing an admissible flux qh,m , with the following: (1) contributions to the CRE estimate obtained from the classical EET technique when computing an admissible flux qh ; (2) contributions to the exact error || e|| 2  evaluated using a highly refined mesh (overkill solution). One observes similarities between the two CRE estimations, showing up areas where the mesh needs to be refined. These areas are correctly predicted when comparing to the exact error distribution. 13. We observe that choosing m = 7 is enough to reconstruct an admissible flux solution, which is equivalent to the one obtained with the classical EET technique. We also observe that m = 3 enables to capture the complexity of the local problems and to provide for a relevant error estimate, even though it is not mathematically guaranteed.

Speedup obtained using the PGD solution

Eventually, we compare the CPU time required to compute the equilibrated flux field depending on which method is used (see Figure 18 for a graphical representation). All the computations were performed on an Intel Core i5 2.4 GHz with When using the EET-PGD technique, the offline CPU cost to compute the PGD solution is 312 seconds; this solution can then be used in a multiquery context. In the online step, computing the equilibrated flux from the classical EET technique (Cholesky factorization) takes 0.01509 seconds per element (0.0587 s for the whole mesh composed of 42 elements), whereas computing the equilibrated flux from a direct evaluation of the PGD solution takes 0.00426 seconds per element (0.0077 s for the whole mesh). We thus observe a speedup of almost 10 in the second CRE step (construction of equilibrated fluxes in each element), and the global speedup on the whole hybrid-flux technique (with associated CPU cost of 0.0960 s) is about a factor 2.

In Figure 14, we represent this speedup for different levels of refinement of the initial mesh (corresponding meshes are given in Figure 15). The speedup increases as the mesh becomes finer, reaching a speedup of 125 on a 2688 elements mesh for the local recovery, while the overall hybrid-flux technique shows a speedup of magnitude 5 on this same mesh. An additional step would be to reduce the CPU time of the tractions reconstruction, by optimizing implementation, to fully benefit from the use of the PGD technique; this is performed in the following section.

Optimization of tractions

We now implement the HFT-PGD technique presented in Section 4, that aims at optimizing the tractions using PGD. We represent in Figure 16 the effectivity index for EET, EET-PGD, and HFT-PGD techniques, with respect to the number m of PGD modes considered. We observe that the global minimization brought by HFT-PGD largely improves the quality of the error estimate and that its tend to 1 when m increases.

In Figure 17, we give the speedup on the overall CPU cost required to compute the error estimate. We observe that HFT-PGD is particularly interesting when considering a small number of PGD modes, even though it always provides for a convenient speedup compared to EET or EET-PGD. Moreover, we observe that the speedup is stable with respect to m when considering EET-PGD alone.

In the following, we consider decompositions with 3 PGD modes. We represent in Figure 18 the normalized CPU time associated with EET, EET-PGD, and HFT-PGD methods when considering the coarsest mesh. The splitting of CPU time contributions emphasizes that the whole cost is associated with the recovery of (optimized) equilibrated tractions when using HFT-PGD.

We represent in Figure 19 the evolution of the speedup for both EET-PGD and HFT-PGD with respect to the number of elements in the mesh. We observe that there is a constant speedup of about 2 for EET-PGD compared to classical EET, for various numbers of elements. HFT, Hybrid Flux Technique; PGD, proper generalized decomposition and a speedup higher than 30 for HFT-PGD; this latter speedup increases with the number of elements, reaching about 100 for 2500 elements.

We eventually analyze the 3 construction techniques, ie, EET, EET-PGD, and HFT-PGD, in terms of CPU time compared to the CPU time required to perform the FEM computation (Figure 20). We observe that contrary to EET and EET-PGD, HFT-PGD leads to an error estimation procedure which is much cheaper than the FEM computation itself; even though an additional global computation is performed in HFT-PGD, the low computational cost is due to the fact that all matrices are precomputed (no numerical integration in the online phase) and that the storage is optimized (sparse matrices are declared). temperature is applied on the external boundary, while a normal flux F d =-1 is applied on hole boundaries. Exploiting symmetries, one quarter of the structure is studied.

Computation of admissible fluxes using PGD

The FE mesh that is used being made of 3-node triangle elements, PGD solutions associated to an element are the same as in Section 5.1. The values of PGD parameters 𝛼, x3 ,andȳ 3 over each element of the mesh are given in Figure 22.

Considering the EET-PGD, ie, using PGD for given tractions obtained from classical equilibration techniques, we represent in Figure 23 still observe that a small number of PGD modes is necessary to reconstruct an admissible flux solution, which is equivalent to the one obtained with the classical EET technique (m = 3 already enables to capture the complexity of the local problems).

Comparison between EET, EET-PGD, and HFT-PGD

We represent in Figure 25 the effectivity index when computing the CRE error estimate with EET, EET-PGD, or HFT-PGD techniques, with respect to the number m of PGD modes. We still observe that HFT-PGD enables to get an effectivity index close to 1 for m ≥ 3. We give in Figure 26 the speedup obtained using EET-PGD or HFT-PGD on the overall CPU time required to compute the error estimate. In the range of values of m that we consider (m ≤ 10), EET-PGD and HFT-PGD are cheaper than EET while HFT-PGD is cheaper than EET-PGD.

In the following, we consider decompositions with 3 PGD modes. We compare in Figure 27 the elementary contributions to the CRE estimate 2E || e|| 2  (evaluated from an overkill solution). We observe that error representations are similar between EET and EET-PGD, and the one given by HFT-PGD is very close to the exact one.

We represent in Figure 28 the normalized CPU time associated with EET, EET-PGD, and HFT-PGD techniques. The associated speedup provided by EET-PGD and HFT-PGD compared to EET is shown in Figure 29; it is given with respect to the mesh size by considering different levels of refinement of the initial mesh (see Figure 30). We observe that there is a constant speedup of about 2 to 3 for EET-PGD compared to classical EET, and a speedup higher than 10 for HFT-PGD; this FIGURE 33 First 3 proper generalized decomposition modes (from left to right) of 𝜌 1 Γ 124 ,20 with from top to bottom: space functions 𝜓 i (x), and parameters functions 𝛿

x 3 i (x 3 ), 𝛿 y 3 i (ȳ 3 ), 𝛿 x 4 i (x 4 ), 𝛿 y 4 i (ȳ 4 ),and𝛿 z 4 i (z 4 )
latter speedup increases with the number of elements, reaching about 30 for 1500 elements. Eventually, when analyzing performance of EET, EET-PGD, and HFT-PGD in terms of CPU time compared to the CPU time required to perform the FEM computation (Figure 31), we observe EET-PGD has a CPU cost comparable to the initial FEM computation whereas HFT-PGD leads to an error estimation procedure which is much cheaper than the initial FEM computation (speedup of about 8).

3D thermal problem

We now consider discretization error estimation on a 3D holed cubical domain Ω, according to a given mesh composed of 4-node tetrahedron elements (Figure 32). The cube has edge length L = 1, and a spherical hole is centered on the cube with radius r = 0.2. We consider a steady-state thermal problem and homogeneous isotropic material properties with  = I. A prescribed zero temperature is applied on the external boundary, while a flux F d = 1isimposedontheinner boundary, ie, on the hole boundary. Owing to problem symmetries, only one-eighth of the structure is studied.

From the associated FE solution, equilibrated tractions may be computed using the first step of the hybrid-flux (or EET) technique.

Details on the PGD solution

Using notations of Figure 3, we compute parameterized solutions 𝜌 j Γ klt ,m (x ref ,𝛼,x 3 , ȳ3 , x4 , ȳ4 , z4 ) of ( 21) with a single fourth-order FE element and 20 PGD modes (m = 20). The domains associated with each geometrical parameter are discretized with 100 points each, after checking that this is sufficient to ensure an accurate description of the evolutions with respect to these parameters. Furthermore, a high-order SVD (HOSVD) technique is applied to the Jacobian terms appearing in (26) to recover an affine dependency with respect to the parameters (about 1000 modes are kept). The first 3 PGD modes of 𝜌 [START_REF] Verfürth | A Review of a Posteriori Error Estimation and Adaptive Mesh Refinement Techniques[END_REF] Γ 124 ,20 are shown in Figure 33. The computation of this PGD solution is done once for all, in an offline phase and stored for later use.

Performance of HFT-PGD

After identifying the PGD parameters over each element of the mesh (see Figure 34), an accurate PGD approximation of the admissible flux qm can then be directly evaluated inside each element in an inexpensive online phase.

In Figure 35, we compare local contributions to the CRE estimate 2E 2 CRE obtained from the HFT-PGD technique (with m = 20) with those obtained from the classical EET technique. The estimate value is 2.83 and the CPU time is 624 seconds when using the EET technique, whereas the estimate value is 2.67 and the CPU time is 219 seconds when using the HFT-PGD technique (speedup of 2.8 and sharper estimate).

We give in Figure 36 the evolution of the value of the error estimate obtained using the HFT-PGD technique (compared to the classical EET technique) with respect to the number m of PGD modes used. We observe that a relevant error estimate is obtained with HFT-PGD for m ≥ 8 and that a sharp error estimate can be recovered considering m ≥ 16. The speedup provided by HFT-PGD compared to EET is shown in Figure 37; it is given with respect to the mesh size by considering different levels of refinement of the initial mesh, and for different values of m.WeobserveherethatHFT-PGD remains competitive compared to classical EET while m ≤ 10. In addition, comparing CPU time with that of the initial FEM computation (Figure 38) confirms that performance of HFT-PGD in terms of computational cost for this 3D case is affected when choosing m ≥ 10. This may be practically changed by using a better implementation and storage of PGD representations.

3D elastoplasticity problem

As a last numerical experiment, we consider error estimation on 3D structure made of a elastoplastic material (Ludwik law) and according to a mesh composed of 8-node hexahedron elements (see Figure 39), as described in Hild et al. [START_REF] Hild | Toward 4D mechanical correlation[END_REF] It corresponds to the numerical simulation of a tensile test on a nodular graphite cast iron sample. The structure is clamped on its bottom side and submitted on its top side to a given displacement field obtained from a linear interpolation of experimentally imposed boundary conditions (see Hild et al [START_REF] Hild | Toward 4D mechanical correlation[END_REF] ). The other sides are free.

The initial mesh is made of 3 × 5 × 10 H8-elements (mesh 1), but refinements with 6 × 10 × 20 elements (mesh 2), 9 × 15 × 30 elements (mesh 3), and 12 × 20 × 40 elements (mesh 4) are also considered. The external loading is applied by means of 10 steps in the nonlinear iterative algorithm. Using the concept of dissipation error introduced in Section 2.2.3, We now assess performance of EET-PGD and HFT-PGD on this specific application. Because of the regular meshes that are used, PGD parameter values remain in a limited interval so that very few PGD modes are required to obtain an

23 FIGURE 1

 231 FIGURE 1Configuration at the level of a 3-node triangle element, with linear tractions FK|Γ defined on element edges

2 FIGURE 3

 23 FIGURE 3Mapping between reference and parameterized elements for a 4-node tetrahedron element

2 FIGURE 4

 24 FIGURE 4 Flow chart of the hybrid-flux approach without (left) or with (right) the use of EET-PGD in step 2. CRE, constitutive relation error; EET, Element Equilibration Technique; PGD, proper generalized decomposition

FIGURE 5

 5 FIGURE 5Flow chart of the classical EET approach (left) versus the HFT-PGD approach (right) to compute the CRE estimate. CRE, constitutive relation error; EET, Element Equilibration Technique; HFT, Hybrid Flux Technique; PGD, proper generalized decomposition

FIGURE 7 FIGURE 8

 78 FIGURE 7 First proper generalized decomposition modes of 𝜌 1 Γ 12 ,20 : space functions 𝜓 i (x) (top), and parameters functions 𝛿 x i (x 3 ) (center) and 𝛿 y i (ȳ 3 ) (bottom) are represented with respect to the order i,fori = 1, 2, 3 (from left to right)

FIGURE 9

 9 FIGURE 9 Identified values of proper generalized decomposition parameters 𝛼, x3 ,andȳ 3 for each element K of the mesh

FIGURE 10

 10 FIGURE 10 Local contributions to the constitutive relation error (CRE) estimate 2E 2 CRE obtained using either the Element Equilibration Technique (EET) technique (left) or the EET-PGD technique (center), and elementary contributions to the exact error || e|| 2  (right). PGD, proper generalized decomposition

  Choosing m = 20 to compute PGD solutions in the EET-PGD technique may be unnecessary. To analyze this point, we show in Figure 11 values of the effectivity index i eff = √ 2E CRE || e||  with respect to the number m of PGD modes used to evaluate the equilibrated flux qh,m . We also represent in Figure 12 the evolution of the relative error || qh,mqh ||  || qh ||  with respect to m, where qh is the equilibrated flux field constructed with the EET technique. A map of || qh,mqh ||  for m = 1, m = 2, and m = 3 is given in Figure

FIGURE 11

 11 FIGURE 11 Value of the effectivity index with respect to the number of proper generalized decomposition (PGD) modes used to evaluate the equilibrated flux in the constitutive relation error estimate. EET, Element Equilibration Technique

FIGURE 12 FIGURE 13

 1213 FIGURE 12 Evolution of the relative error || qh,mqh ||  ∕|| qh ||  with respect to the number of proper generalized decomposition modes used to compute qh,m

FIGURE 14 FIGURE 15

 1415 FIGURE 14 CPU speedup on the equilibrated flux step for meshes with increasing numbers of elements (speedup on the local recovery alone or on the overall hybrid-flux technique are explicited)

FIGURE 16 FIGURE 17

 1617 FIGURE 16Evolution of the effectivity index with respect to the number of proper generalized decomposition (PGD) modes considered, for Element Equilibration Technique (EET), EET-PGD, and HFT-PGD. HFT, Hybrid Flux Technique

FIGURE 18 FIGURE 19

 1819 FIGURE 18Comparison of CPU times required to compute the error estimate depending on the method which is used. EET, Element Equilibration Technique; HFT, Hybrid Flux Technique; PGD, proper generalized decomposition

FIGURE 20

 20 FIGURE 20Ratio of CPU times for Element Equilibration Technique (EET), EET-PGD, and HFT-PGD with respect to the FEM CPU time,

FIGURE 21 FIGURE 22 FIGURE 23 FIGURE 24 FIGURE 25 FIGURE 26

 212223242526 FIGURE 21Representation of the 2D. domain, associated FE mesh, and applied loading

FIGURE 27 FIGURE 28 FIGURE 29

 272829 FIGURE 27 Local contributions to the exact error || e|| 2  and to the constitutive relation error estimate 2E 2 CRE obtained using Element Equilibration Technique (EET), EET-PGD, or HFT-PGD techniques (from left to right). HFT, Hybrid Flux Technique; PGD, proper generalized decomposition

FIGURE 34

 34 FIGURE 34 Identified values of proper generalized decomposition parameters 𝛼, x3 , ȳ3 , x4 , ȳ4 ,andz 4 for each element K of the mesh

FIGURE 35

 35 FIGURE 35 Local contributions to the constitutive relation error estimate 2E 2 CRE obtained using either the Element Equilibration Technique technique (left) or the HFT-PGD technique (right). HFT, Hybrid Flux Technique; PGD, proper generalized decomposition

FIGURE 36 FIGURE 37 FIGURE 38 FIGURE 39

 36373839 FIGURE 36 Evolutions of the error estimates obtained with Element Equilibration Technique (EET) and HFT-PGD techniques with respect to m. HFT, Hybrid Flux Technique; PGD, proper generalized decomposition

FIGURE 40

 40 FIGURE 40 Local contributions to the CRE estimate 2E 2 CRE obtained using the Element Equilibration Technique (EET) technique associated with meshes 1, 2, 3, and 4 (from left to right)

  The force associated to X, involving the additional yield stress R for isotropic hardening, is denoted Y, and we introduce quantities s =( 𝜎 𝜎, Y) T , ̇ep =(ε p , -Ẋ) T ,a n d ̇ee =(ε e , Ẋ) T , so that the dissipation reads d = ∫

	state equations e e =Λ ( s)= 𝜕Ψ * 𝜕s (or s = 𝜕Ψ 𝜕e e	T 0 ∫ Ω s • ̇ep . The material behavior is described by

  x 3 , ȳ3 , x4 , ȳ4 , z4 ) are local coordinates of nodes 3 and 4 in the coordinate system associated with element K,and (𝜂, 𝜉, 𝜏) are local coordinates in the coordinate system associated with element K ref . The global mapping  thus involves 6 parameters, ie, p geo =(𝛼, x3 , ȳ3 , x4 , ȳ4 , z4 ) and J = 𝛼 3 ȳ3 z4 .Consequently, the problem (21) can be parameterized and defined in a fixed reference element K ref under the form:

2

  CRE (obtained from EET, EET-PGD, or HFT-PGD techniques) with contributions to the exact error FIGURE 30 Sequence of refined meshes used to assess CPU speedup FIGURE 31 Ratio of CPU times for Element Equilibration Technique (EET), EET-PGD, and HFT-PGD with respect to the FEM CPU time, for various numbers of elements. HFT, Hybrid Flux Technique; PGD, proper generalized decomposition

FIGURE 32

Representation of the 3D domain, associated FE mesh, and applied loading

accurate solution over each element. We choose m = 5 in practice, and we report in Table 1 the performance of EET-PGD and HFT-PGD in terms of effectivity index and speedup compared to classical EET, for the 4 meshes described previously. We observe large speedups for both methods, due to the fact that a limited number of PGD modes is considered and that the equilibration procedure has to be performed at each iteration step.

CONCLUSIONS AND PROSPECTS

We presented a general framework that highlights the original and beneficial use of PGD in verification procedures performed by means of the CRE concept. Based on an offline/online strategy, this framework involves PGD representations at the level of each element of the finite element mesh. This drastically decreases the computational cost and technicalities that are essentially associated with the computation of admissible fields. It thus facilitates practical implementation, decreases CPU time, and increases the quality of the CRE error estimate.

As the proposed technique is focused on balance equations alone, it directly to any nonlinear time-dependent mechanical problem. In this framework, we forecast a large speedup in industrial applications due to the fact the construction of admissible fields needs to be performed at each time step and that the speedup increases with number of elements. We thus believe this work the to robust, practical, and real-time methods for controlling FEM simulations performed by Computational Mechanics software.

Additional on the proposed could be investigated: application to higher-order elements on commeshes; (2) consideration as a tool for designing optimal meshes, as well as for mesh refinement on the fly; and (3) comparison with concurrent approaches that exhibit an offline stage such as the patch-based mixed finite element approach. [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF][START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF] All these points will be the topics of forthcoming research works.