José Aguado

Luisa Silva

Keywords: linear interpolations of primal dictionary, Hertz problem

The journey of this thesis began in about three and a half years ago. During this journey, many people who have supported me and my work,

ºÚ a Ú Ú я ।

Elasticity is a property in the bodies having densely populated elemental parts; when sheared by deforming forces, it helps the body recover its original state.

Sridhara Acharya (870 CE -930 CE)

. . . dass zwei elastiche isotrope Körper sich in einem sehr kleinen Theil ihrer Oberfläche benuhren, und durch disen Their einen endlichen Druch der eine anf den andern ausüben . . . Die Fragen, deren Beantwortung uns naturgemäss zunächst obleigt, sind . . . die Frage nach der Form und absoluten Grösse der Druckfigur.

. . . the case of two elastic isotropic bodies which touch each other over a very small part of their surface and exert upon each other a finite pressure, distributed over the common area of contact. . . . The questions that we first seek answers to are . . . What is the form and the absolute magnitude of the curve of pressure ? Heinrich Hertz (1857-1894) 4.11 Projection error of nested level n + 1 snapshots on dictionary of level n . 4.12 Comparison of reconstruction error for the Hertz and ironing problem using the original dual dictionary and DTW-adapted dual dictionary . . .

A.1 Reconstruction instance with spurious snapshot selection by the greedyactive set algorithm if non-penetration condition is applied "hardly" . . .

List of Tables

Introduction

In the manufacturing of engineering structures, process variabilities often lead to deviations from the tolerance specifications on the features of a part or an assembly. These deviations can influence the stresses at the contact interfaces in the assembly, potentially leading to reduced performance or even failure of the assembled structure. Economic constraints often provide the motivation to salvage the "non-conforming" hardware. Depending on the nature of the deviations, a non-conforming part or assembly can be either: accepted, reworked to improve compliance or rejected. The decision is usually based on further analyses or tests on the non-conforming part. Naturally, our understanding of the mechanics of the contact interfaces in assemblies is crucial in this decision-making process.

Analytical and empirical solutions are usually not available for complex mechanical/structural configurations, while experimental testing is cost-prohibitive. Therefore, numerical modelling is often the most practical route to aid the decision-making process. However, the numerical simulations consume a significant amount of time and computational resources in building the model and subsequent resolution, especially for assemblies that include dynamic contacts generating additional complexities. Also, the communication delay between the shop floor and the simulation team adds to the delays. These bottlenecks have been schematically described in Fig. I. Therefore, it is necessary to accelerate assembly simulations in order to comply with desired decision times in the manufacturing cycle.

Real-time simulation tools have been of interest in various decision-making processes across fields like design and health-monitoring, among others. In such applications, fast simulations need to be performed for parametric values that are fed to the simulator in real-time. However, the complexity of a model usually increases exponentially with the number of parameters in the model. Therefore, high-fidelity models are often replaced with surrogate models that offer fast simulations, sometimes referred to as digital twins. Reduced Order Modelling (ROM) a.k.a. Model Order Reduction, is a well-known framework to create such models. ROMs are generally built to efficiently solve a parametric problem in a multi-query context. Therefore, ROM is a promising framework for moving towards efficient simulations of non-conforming engineering hardware.

For a real-time simulation of non-conforming assemblies, the task is two-fold. the modelling of geometric variabilities in individual parts. Shape parametrized ROMs have been explored [START_REF] Lauzeral | A model order reduction approach to create patient-specific mechanical models of human liver in computational medicine applications[END_REF][START_REF] Marconi | A nonlinear reduced order model with parametrized shape defects[END_REF][3] which can be used for parametrizing geometric variabilities. The second, and more challenging, is the application of ROM approach to contact mechanics. As contact mechanics deals with multi-body problems with interaction between the boundaries of each body, it introduces geometric non-linearities. These non-linearities are different in nature compared to, say, non-linear material behaviour or loading conditions which have been widely studied by the ROM community. Therefore, contact ROMs still need extensive research before they can be used in industrial real-time applications, such as the acceleration of the decision-making process discussed previously. Therefore, in this thesis, the current limitations of contact ROMs have been investigated and corresponding mitigation strategies have been proposed. The thesis begins with an introduction to contact mechanics and the state of the art in ROMs applied to contact mechanics (Chapter 1). Then, the limitations of the so-called low-rank approach are investigated (Chapter 2). Strategies to counter these limitations are developed using over-complete dictionaries (Chapter 3) and non-linear interpolations (Chapter 4).

In this chapter, an overview of concepts and state of the art related to contact mechanics and reduced order modelling are presented. Ideas in contact mechanics are first introduced using a two-body problem where inequality constraints apply to the energy minimization statement of the standard mechanics problem. Some of the classical methods to solve inequality constrained problems are discussed. As one of the principal sources of complexity in such problems is contact detection, some methods applicable to contacts between one-dimensional boundaries are introduced. In the second part of the chapter, the concept of Reduced Order Modelling and its application on unconstrained problems is discussed. In the final part of the chapter, the two elements are combined, where the developments in reduced models of contact mechanics are presented.

Overview of Contact Mechanics

Many computational mechanics problems involve interaction between multiple bodies with various types of physics, say thermal or mechanical, where each kind of physical interaction adds a challenge to the modelling process. Contact mechanics, the field that considers mechanical interaction in multi-body problems, finds applications in a wide variety of fields as almost all mechanical system, natural or man-made, displays some sort of contact phenomena. A few examples that fall in this category include: foundation studies in civil engineering, human joints in biomechanics, roller joints in mechanisms, adhesive and frictional contacts in tribology and crash simulations in automotive industry.

Historical works in contact mechanics include studies by that of Hertz [START_REF] Hertz | Ueber die Berührung fester elastischer Körper[END_REF] on analytical solutions of contact between spherical frictionless surfaces of elastic bodies, that of Coulomb [START_REF] Coulomb | Essai sur une application des regles de maximis et minimis a quelques problemes de statique relatifs a l'architecture[END_REF] on frictional contact, and many more. However, the work by Signorini [START_REF] Signorini | Questioni di elasticità non linearizzata e semilinearizzata[END_REF] on the contact problems with elastic bodies, which consisted of the Boundary Value problem with non-penetration constraints, is central to approach of contact mechanics in the computational field. The problem was titled "problems with ambiguous boundary conditions." by Signorini, aptly describing the nature of the problem. Analytical solutions to certain contact problems have been computed in [START_REF] Johnson | Contact Mechanics[END_REF][START_REF] Johnson | Introduction to Contact Mechanics: a summary of the principal formulae[END_REF][START_REF] Goryacheva | Contact Mechanics in Tribology[END_REF]. x 2 will possibly come into contact with point x1 , and hence considered as a contact pair. The distance between the two points in the initial and deformed configuration is given by g(µ, v)(x 2 ) and g(µ, v)(x 2 )k(µ, v; v)(x 2 ) respectively

Parametrized contact mechanics problems

Usually, an engineering problem is described as parametrized if a parametric solution is sought i.e. the solution of the problem is sought for a range of parameter values. For a mechanical problem, the parameters can be any attribute of the mechanical system such as material properties, geometry, loading or boundary conditions, etc. Parametric solutions usually find applications in design optimization, real-time monitoring and inverse problems (see Section 1.2). In this section, parametric mechanical problems with non-penetration constraints on displacement field are discussed. Linear elastic and small deformation problems are considered, keeping the evaluation of internal energy straightforward. The inequality constraints may be a result of the presence of either an obstacle or a second body in the domain, with possibility of contact between different surfaces in the deformed configuration.

For convenience, the contact problems will be referred to as Type-1, referring to a problem involving a deformable body and an immovable obstacle, and Type-2 referring to problems involving two deformable bodies. A schematic of the two types of contact problems is shown in Fig. 1.1. Contact phenomena like friction and adhesiveness are neglected. The resolution of contact mechanics problems depends on computing the distance between surface of bodies involved. As evaluation of distance is not a trivial task because of dynamic contact pairing, different strategies exist based on underlying simplifications. The details to these approaches can be found in [START_REF] Wriggers | Computational Contact Mechanics[END_REF][START_REF] Fischer | Mortar based frictional contact formulation for higher order interpolations using the moving friction cone[END_REF][START_REF] Yastrebov | Numerical Methods in Contact Mechanics[END_REF].

The generic weak form of a mechanical problem with a parameter µ involving contacts (for both Type-1 and Type-2) is described by the following inequality constrained minimization problem:

u = arg min v∈V 1 2 a(µ; v, v) -f (µ; v) s.t. k(µ, v; v)(x) ≤ g(µ, v)(x) on Γ 2 (1.1)
where • a(µ; v, v) is strain energy function and f (µ; v) is the work contribution by the external forces. The operators are defined ∀v ∈ V and ∀µ ∈ P, where V be an appropriate function space for displacement v and P is the parametric space for the parameter µ. Such problems with parametric dependence are usually referred to as parametrized problems. Note that µ may not be necessarily a scalar and the parametric space P can be possibly multidimensional.

• k(µ, v; v)(x) and g(µ, v)(x) are the distance functions that indicate separation between bodies, defined conveniently on surface Γ 2 . The first term indicates the contribution of displacement field v to the distance and the second term indicates the distance in the undeformed configuration. For a Type-2 problem, the distance functions can be expressed as:

k(µ, v; v)(x) = v 2 (x 2 ) -v 1 ( x1 ) • n g(µ, v)(x) = x 2 -x1 • n (1.2)
The pair (x 2 , x1 ) and the normal n are determined on the discrete domain using various methods, some of which are given in upcoming Section 1.1.3.1. For a Type-1 problem, the displacement field of the obstacle v 1 vanishes.

Notation 1.1 :

Notations to distinguish strong and weak forms are adopted from [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]. Linear and bilinear operators that appear in weak formulations are expressed in the format o(•; •, •), where the non-linear dependencies and linear dependencies are separated by the semi-colon ';' . On the other hand, strong form functions are expressed in the format o(•; •, •)(•), distinguished from weak forms using an additional argument which is the geometric position. The same applies to linear forms expressed as o(•; •) and o(•; •)(•), where the linear dependency is on the second argument.

Contact Formulations

To treat the contact problem (1.1) numerically, the contribution of contact to the total energy of the system is computed. This additional energy term is defined differently in each contact formulation, however, the basic idea is to introduce contact forces, directly or indirectly, that will prevent penetrations of contact surfaces. Among numerous formulations available, three of them are introduced in this section. Two of them are based on introducing the Lagrange Multiplier, which computes the contact forces directly, while the other method penalizes the penetrations.

Lagrange Multiplier Method

One of classical methods is an attempt to accurately satisfy the inequality constraints of the contact problem (or any minimization problem) is the Lagrange Multiplier Method (LMM). The resulting optimization problem can be expressed as a saddle point problem of the Lagrangian functional:

(u, λ) = arg minmax v∈V,η∈W + L(v, η) (1.3)
where

L(v, η) = 1 2 a(µ; v, v) -f (µ; v) + b(µ, u; η, u) -d(µ, u; η) (1.4)
where the distance function is expressed by following weak forms:

b(µ, u; η, u) = Γ 2 η k(µ, u; u)(x) ∂Γ d(µ, u; η) = Γ 2 η g(µ, u)(x) ∂Γ (1.5)
where W + is an appropriate non-negative function cone for contact pressure defined on the surface Γ 2 . The inequality λ ≥ 0 admits a physical meaning because λ is equivalent to contact pressure. A negative contact pressure implies traction, in other words, adhesion between contact surfaces which contradicts the simplifying assumptions.

The KKT conditions associated to the optimizer (u, λ) ∈ V × W + of (1.3) can be expressed as follows:

a(µ; v, u) -f (µ; v) + b(µ, u; λ, u) = 0, v ∈ V (1.6a) b(µ, u; η, u) -d(µ, u; η) ≤ 0, η ∈ W + (1.6b) b(µ, u; λ, u) -d(µ, u; λ) = 0 (1.6c)
Notice that the last KKT conditions (1.6c) implies that contribution of contact energy term will be zero at the solution. This makes sense because the non-penetration constraint can prevent the system from reaching the true minima of the energy functional, but it cannot add energy to the system.

Penalty Method

To constrain the solution in the feasible region defined by the inequality constraints, violations of inequality constraints can be penalized by adding a penalty term to the energy of the system. The resulting system is an unconstrained problem, unlike the LMM problem.

u = arg min v∈V 1 2 a(µ; v, v) -f (µ; v) + Γ 2 τ (k(µ, v; v)(x) -g(µ, v)(x)) + 2 ∂Γ (1.7)
where (z) + := z if z > 0 0 otherwise The solution of (1.7) may not necessarily satisfy the inequality constraints exactly, but may be nearly compliant for sufficiently large values of the penalty parameter τ . For very large values of τ , the penalty solution should approach the LMM solution, however, the numerical system are ill-conditioned in this regime of τ .

The penalty parameter τ can be interpreted as a linear spring that resists penetration. Therefore, a natural variation of the penalty method can be built by using a non-linear penalty function where the penetration and the resultant contact pressure have a nonlinear relation. This approach is useful to reduce the penetration effects, but also in cases where a constitutive law relating contact pressure and penetration is known. However, the problem of ill-conditioning still persists when the constitutive law enters a highly stiff regime [START_REF] Wriggers | Computational Contact Mechanics[END_REF]Chapter 5,[START_REF] Signorini | Questioni di elasticità non linearizzata e semilinearizzata[END_REF].

Augmented Lagrange Multiplier Method

To circumvent the ill-conditioning problems of penalty method, Augmented Lagrange Multiplier Method [START_REF] Hestenes | Multiplier and gradient methods[END_REF] problems have been developed. This class of method use a combination of penalty and LMM. The basic idea is to establish contact forces using the penalty method and then transfer these forces, in multiple steps, into the Lagrange Multiplier, thereby reducing illegal penetrations. The advantage of using this method is that the solution will have high compliance to the constraints irrespective of the magnitude of the penalty parameter. At a given step k, the following problem is minimized:

u k = arg min v∈V 1 2 a(µ; v, v) -f (µ; v) + Γ 2 τ (k(µ, v; v)(x) -g(µ, v)(x)) + 2 ∂Γ + b(µ, v; λ k , v) -d(µ, v; λ k )
with the update step:

λ k+1 (x) = λ k (x) + τ (k(µ, v; v)(x) -g(µ, v)(x)) +
and initial conditions:

λ 0 (x) = 0
As λ k does not contain the true Lagrange Multiplier during the intermediate steps, this method still allows some penetration. However, as the Lagrange Multipliers are updated with new contact pressure forces, the penetration reduces in the subsequent steps. With sufficient number of steps, penetration diminishes to very low-levels and the solution approaches that of the LMM solution.

Discrete Formulations

In this section, the usual discretization strategies based on Finite Element (FE) are briefly introduced. One of the foremost problems in contact mechanics is the detection of contact pairs. For large-scale problems, this step can prove to be expensive. Then comes the part of building the finite element (FE) operators that can be used to enforce or verify the constraints. As the contact problems in upcoming chapters will be solved by using the LMM approach, the discussion on discretization schemes will be limited to this method. Discretization of the penalty and Augmented LMM can be carried out in a similar manner. Details of all approaches can be found in [START_REF] Wriggers | Computational Contact Mechanics[END_REF][START_REF] Yastrebov | Numerical Methods in Contact Mechanics[END_REF].

Contact Detection and Pairing

Establishing the region of contact between two bodies is one of the primary challenges in contact problem. In simplified cases of small displacement, conforming meshes can be used, thereby allowing the use of node-to-node contacts. In this case, both points in the contact pair x 2 and x1 are nodes over respective surfaces, as shown in Fig. 1.2. Naturally, the contact pair can be predetermined based on the discretization. Obviously, the use of conforming mesh limits the applicability of this method. For example, it cannot be applied to sliding surfaces. Node-to-segment method is used in cases of non-conforming meshes with 1D boundaries, for e.g., 2D problems with possibility of contact on the edges of bodies. The contact pairings (x 2 , x1 ) are not straightforward to compute as it was in node-to-node method. The contact pair evaluation is one of the sources of complexities in contact problems, especially for large displacement problems as contact pairs can be a strong function of displacement field. Established methods of computing contact pairs include:

• Closest Point Projection: For each node x 2 on Γ 2 , the closest point x1 on Γ 1 is chosen as the pair [START_REF] Hallquist | Sliding interfaces with contact-impact in large-scale lagrangian computations[END_REF]. Naturally, the distance functions are evaluated using (1.2) with n = n1 , the normal at the point x1 . A limitation of this method is the possibility of non-unique closest point in certain geometries [START_REF] Yastrebov | Numerical Methods in Contact Mechanics[END_REF]Chapter 2].

• Ray Tracing: At each node x 2 of Γ 2 , the local normal n 2 is extended to intersect Γ 1 at x1 [START_REF] Poulios | An unconstrained integral approximation of large sliding frictional contact between deformable solids[END_REF]. Hence, distance functions are evaluated using n = -n 2 for this method. This method may not find the nearest point, but avoids the problems of non-uniqueness and projecting on irregular points.

The two methods are shown in Fig. 1.3. The normals n 1 and n 2 are evaluated in deformed configuration. The pairing ( x1 , x 2 ) is hence dependent on the displacement field in both methods. This forces the distance function k to have an implicit non-linear dependence on displacement, apart from the explicit linear dependence that is evident in the expression. The discrete equivalent of the Lagrangian functional in (1.4) can be derived by introducing finite element spaces V h ⊂ V and W h + ⊂ W + :

L h (v, η) = 1 2 v T K(µ)v -v T f (µ) + η T [C(µ, v)v -g(µ, v)] (1.8)
where v ∈ R N h u and η ∈ R N h λ . Here, N h u and N h λ indicate the number of respective finite element degrees of freedoms (dofs). Also, note that the notation u and v are used for both continuous and discrete versions of displacement fields.

Let {N 1 , N 2 . . . N N h u } and {M 1 , M 2 . . . M N h λ
} denote the basis for the FE space V h and cone W h + respectively. Here, for a system with n displacement dofs per node, N i : Ω → R n and M j : Γ → R. The discrete operators K(µ) and C(µ, v) can be defined as:

K(µ) ij = a(µ; N i , N j ) C(µ, v) ij = b(µ, v; M i , N j )
and the vectors f and g are defined as:

f (µ) i = f (µ; N i ) g(µ, v) i = d(µ, v; M i )
The bilinear and linear terms b and d are non-linear w.r.t. the displacement field v and the parameter µ, which is reflected in discrete form as C(µ, v) and g(µ, v) respectively. Using the above operators, the discrete form of the KKT conditions in (1.6) for the optimal solution (u, λ) ∈ R N h u × R N h λ can be written as:

K(µ)u -f (µ) + C T (µ, u)λ = 0 (1.9a) λ ≥ 0 (1.9b) C(µ, u)u -g(µ, u) ≤ 0 (1.9c) λ T [C(µ, u)u -g(µ, u)] = 0 (1.9d)
Apart from the non-linearity of the constraint operators C and g, additional nonlinearities appear because the regions where the inequality constraints are active and inactive are unknown. Active regions correspond to the regions where contact is established, or in other words, constraints satisfy the equality. These non-linearities are potential contributors to the overall computational complexity.

Stability of the saddle point problem:

The stability conditions of mixed formulations usually appear in the form of an inf-sup condition. The inf-sup constant of the bilinear operator b in (1.4) must be positive for every pair of test functions in the two solution spaces. Not every combination of finite-element spaces V h and W h + satisfies the stability condition. In this thesis, linear elements for displacement and node-centered zero order (piecewise-constant) shape functions for contact pressure are used (based on Reference [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]), which generates smooth solutions for two-dimensional problems. More details on inf-sup stability of contact problems can be found in [START_REF] Wriggers | Computational Contact Mechanics[END_REF][START_REF] El-Abbasi | Stability and patch test performance of contact discretizations and a new solution algorithm[END_REF].

Resolution schemes

Even for contact problems involving small deformations with linear elastic constitutive law, geometric non-linearities can still occur if large displacements are involved. This is because the contact pairs change with the displacement field. Even for small displacement problems where contact pairs can be approximately precomputed in the reference configuration, the region of contact cannot be predetermined. Linearization schemes are needed to solve the discrete problem, that resolve both the contact region and the contact pairs in an iterative manner. Two such schemes: Fixed-point and Newton-Raphson iterations are described here. Python notation is used to indicate slicing operations.

Fixed point method

Fixed point schemes are used to find solutions to non-linear equations by posing them in the fixed-point form x = f (x), leading to a fairly simple iterative form x p+1 = f (x p ). For a discrete system with multiple dofs, the linearized form usually looks like A(x p )x p+1 = b(x p ). These linearized forms can be iterated until convergence, without computing derivatives or gradients of the non-linear equation.

Fixed point methods for an inequality constrained problem can be implemented using active set approach [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Luenberger | Linear and Nonlinear Programming[END_REF], in which the constraints are divided into active and inactive sets. For the discretized contact problem, the balance form (1.9a) can be used in the fixed-point method:

K(µ)u + C T (µ, u)λ = f (µ) (1.10)
Given the current estimate of active set A, the activation state of the inequality constraints can be written as:

λ[A c ] ← 0 (1.11a) C A (µ, u)u -g A (µ, u) = 0 (1.11b)
where

C A (µ, u p ) = C(µ, u p )[A, :] g A (µ, u p ) = g(µ, u p )[A]
and A c is the complement of set A. These conditions apply because each entry of the vectors C(µ, u)ug(µ, v) and λ correspond to a nodal constraint and a nodal Lagrange Multiplier respectively. As the active set contains the list of active nodes, it can be used for filtering the entries of constraint operators and enforcing the equality conditions i.e., where contact between bodies is established. The displacement field and the active Lagrange multipliers at iterative level p+1 can be obtained using the linear system obtained from the equilibrium condition (1.10) and the active constraint (1.11b)

K(µ) C T A (µ, u p ) C A (µ, u p ) 0 u p+1 λ p+1 A = f (µ) g A (µ, u p ) (1.12)
and the discrete contact pressure can be updated as:

λ p+1 [A] ← λ p+1 A λ p+1 [A c ] ← 0
The sign of entries in the expression C(µ, u)u-g(µ, v) and in the vector λ correspond to satisfaction or violation of inequality constraints at the respective node's support. The active set for the next iteration is obtained by verifying the constraints:

A i ← i C(µ, u p )u p+1 -g(µ, u p ) i ≥ 0 and λ i > 0 (1.13)
Moreover, the contact pairs are also be updated using the latest displacement field and the new operators C(µ, u p+1 ) and g(µ, u p+1 ) are computed. This process is continued until convergence.

Newton-Raphson

The Newton-Raphson iterative method (a.k.a. Newton's method) is a derivative-based iterative method for root finding and optimization problems. For finding the roots of an equation f (x) = 0, the Newton-Raphson update at an iterative level p can be written as ∆x p+1 = -f (x p ) f (x p ) , where ∆x p+1 = x p+1x p . For a discrete system with residual b(x) -A(x)x, the Newton-Raphson update step can be written as:

A T (x)∆x p+1 = b(x p ) -A(x p )x p
where A T (x) is the tangent matrix of the residual.

The Newton-Raphson method for contact problems is based on computing the tangent matrix K c T of the residual of equilibrium equation (1.10) and then computing incremental updates. The contribution to the tangent matrix K T comes from of the contact part of the energy functional or the equilibrium equation. It consists of the second order terms that appear due to variations in contact pairing. The computation of tangent matrix for the non-linear system is quite detailed and can be found in Reference [START_REF] Wriggers | Computational Contact Mechanics[END_REF]Chapter 9,[START_REF] Wriggers | Computational Contact Mechanics[END_REF]. The iterative update at a level p + 1 can be written as:

K(µ) + K c T (µ, u p ) C T A (µ, u p ) C A (µ, u p ) 0 ∆u p+1 ∆λ p+1 A = f (µ) g A (µ, u p ) - K C T A (µ, u p ) C A (µ, u p ) 0 u p λ p A
The stiffness matrix K(µ) appears in addition to the contact tangent matrix because the stiffness matrix is the tangent of the elastic energy functional. After computing the increment, the state vector can be updated as:

u p+1 λ p+1 A = u p λ p A + ∆u p+1 ∆λ p+1
A This is followed by the active set update in (1.13).

Reduced Order Modelling for parametric problems

In recent applications of numerical methods to engineering design, decision-making processes are increasingly dependent on simulation results [START_REF] Guratzsch | Structural health monitoring sensor placement optimization under uncertainty[END_REF][START_REF] Modesto | Proper generalized decomposition for parameterized helmholtz problems in heterogeneous and unbounded domains: Application to harbor agitation[END_REF][START_REF] Mosavi | Engineering Design and Decision Making Models[END_REF][START_REF] Bigoni | Simulation-based anomaly detection and damage localization: An application to structural health monitoring[END_REF]. Various configurations of the product can be explored and compared by optimization algorithms to develop an optimal configuration using modern simulation tools [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF][START_REF] Nguyen | A review on simulation-based optimization methods applied to building performance analysis[END_REF][START_REF] Woldemariam | Numerical Simulation-Based Effect Characterization and Design Optimization of a Micro Cross-Flow Turbine[END_REF]. For the integration of simulation tools into the decision-making process, simulation time for each configuration becomes a primary factor in determining the swiftness of the process. In this context of multi-query simulations, a high computational complexity of traditional simulation methods might discourage comprehensive exploration of the design space for economic reasons.

Traditional methods of simulations were restricted to the design phase of engineering developments due to the heavy computational costs. On the other hand, modern simulation tools capable of real-time simulations have promising applications not only in the design phase but also during the life-cycle of engineering products, for example, for operations and maintenance [START_REF] Hartmann | Model Order Reduction a Key Technology for Digital Twins[END_REF]. The idea behind such applications is to create real-time simulations whose predictions can be used to initiate corrective or preventive actions. At the heart of the development of real-time simulations lies the approach of Reduced Order Modelling (ROM).

Real-time simulators need to mimic the relevant physics of the system being studied taking into account the parameters that are fed to the simulator in real-time. Highfidelity numerical methods usually have very high computational complexity for such applications. The real-time query can be addressed by solving a high-dimensional problem, in which the real-time parameters are treated as coordinates of the problem (like space and time). Such a numerical model is usually referred to as a parametrized model and is usually very expensive to resolve using standard numerical methods. ROM is a framework for reducing the complexity of high-dimensional problems. ROMs have found successful applications in a wide variety of problems such as design optimization [START_REF] Mcbane | Component-wise reduced order model lattice-type structure design[END_REF], uncertainty quantification [START_REF] Freitas | Fast and accurate model reduction for spectral methods in uncertainty quantification[END_REF], inverse problems [START_REF] González | Proper generalized decomposition based dynamic data driven inverse identification[END_REF], optimal control [START_REF] Kärcher | Certified reduced basis methods for parametrized elliptic optimal control problems with distributed controls[END_REF], real-time monitoring [START_REF] Modesto | Proper generalized decomposition for parameterized helmholtz problems in heterogeneous and unbounded domains: Application to harbor agitation[END_REF], and others.

Various methodologies have been adopted to create ROMs to address different bottlenecks. Some are based on reducing the number of dofs of the system, whereas others focus on the efficient evaluation of integrals over a small part of the domain. However, a common feature among most ROM methods has been the computation of a training set of snapshots i.e. a set of high-fidelity solutions based on which a reduced model is built. A brief description of a few classes of reduction methods is given here:

• Low-rank methods: These methods primarily focus on finding a low-rank solution subspace. Usually, this is done by computing the basis of the subspace using Singular Value Decomposition (SVD) or Proper Orthogonal Decomposition (POD). Low-rank methods are further detailed in Section 1.2.1.

• Hyper-reduction methods: Primary aim of hyper-reduction methods is to reduce the number of unknowns as well as quadrature points to compute weak form terms by generating a reduced integration domain (RID), which is generated by selecting appropriate dofs [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF]. This feature of hyper-reduction is used in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF] to reduce the effort in the computation of internal variables of non-linear constitutive models, where the RID is made up of regions containing maximum energy of POD basis vectors. On the other hand, the quadrature rules for RID in Reference [START_REF] Hernández | Dimensional hyper-reduction of nonlinear finite element models via empirical cubature[END_REF], are computed using the least-square procedure so that all training samples are exactly integrated. Cheaper evaluation of weak forms over the RID leads to an efficient model for non-linear problems.

• Non-intrusive methods: Some methods have been developed that do not need access and modification of the high-fidelity source codes, and are termed as nonintrusive methods. These are highly useful if proprietary source codes or complex constitutive models are used. Also, it can be applied to different kind of physics without changing the non-intrusive framework. As the numerical model and its operators are not accessible, many non-intrusive methods rely on fitting an interpolation in the parametric domain. Radial basis functions (RBF) have been popularly used in this context. In these works, a POD-based solution subspace is computed and an RBF-based interpolation for the reduced coefficients is developed for the POD subspace [START_REF] Audouze | Reduced-order modeling of parameterized pdes using time-space-parameter principal component analysis[END_REF][START_REF] Audouze | Nonintrusive reduced-order modeling of parametrized time-dependent partial differential equations[END_REF][START_REF] Xiao | Non-intrusive reduced-order modelling of the navier-stokes equations based on rbf interpolation[END_REF]. Another approach is that of sparse subspace learning where a predefined set of hierarchical parametric basis functions with special quasi-interpolative properties, allowing enrichment of the reduced model at minimal cost [START_REF] Borzacchiello | Non-intrusive sparse subspace learning for parametrized problems[END_REF]. The CUR decomposition method takes a different approach for generating the parametric interpolation. A set of selected spatial and temporal snapshots are used to build a surrogate model with coefficients that fit the snapshots in the training set [START_REF] Le Guennec | A parametric and non-intrusive reduced order model of car crash simulation[END_REF].

• Projection-based a priori method : The Proper Generalized Decomposition (PGD) treats the parameters of the problem as coordinates and yet circumvents the curse of dimensionality by using separated representations for solving the PDE [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids[END_REF]. This is accomplished by means of a greedy approach, where the current solution is updated a with rank-1 tensor product of separated functions in each step. The rank-1 updates are computing using the separated form of PDE, one direction at a time, where each direction corresponds to one or more coordinates of the problem, for instance, position, time, or a parameter, or a combination of these. This is described as an a priori method because a parametric solution is built and stored, resulting in quick accessibility of the solution. Efficient computation of separated form of non-linear terms has been implemented in the PGD framework using cross-approximations [START_REF] Aguado | Tensor Representation of Non-linear Models Using Cross Approximations[END_REF].

Even though these approach have different names and methodologies, all of them, directly or indirectly, assume the low-rank hypothesis. In the next section, low-rank methods are described for generic unconstrained problems.

Low-rank methods

In many numerical problems, the underlying structure lies in a subspace whose dimensionality is quite small compared to the full dimensionality of the solution space used in the numerical model. Examples include thermal problems, transport problems, elasticity and fluid flows among many others [START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF][START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF]. Low-rank methods can be used to reduce the complexity of such problems. The idea behind low-rank ROMs is to split the cost of computation into two stages.

• The offline stage, where most of the computational complexity is resolved, consists of computing a set of solutions in the parametric space i.e. the training set, and then computing the reduced subspace that defines the span of training set. The idea is that space spanned by the training set snapshots is a good approximation to the true solution subspace. This is usually done by computing a reduced basis (RB) corresponding to the span of the training set. This step is typically performed only once to create the reduced model.

• The online step is performed every time a new query is submitted, in which the ondemand solution is computed using the reduced model which is cheaper to evaluate than the high-fidelity model.

Using notation from the previous section, let the high-fidelity (FE) solution space be denoted by V h and the reduced subspace be denoted by V r = span(Φ), where Φ denotes the corresponding RB. In online stage, solutions in the subspace V r are sought to approximate the high-fidelity solution. In discrete sense, Φ is a matrix with each column corresponding to a basis vector. Then, the approximation in the subspace can be expressed as:

u ≈ Φ u u ∈ R N r u
All reduced degrees of freedom and reduced operators will be indicated by a hat • For constructing the RB Φ in offline stage, a finite subset of the parametric space P tr ∈ P is explored and an approximate subspace V r is extracted from the set of highfidelity snapshots from the training set. The number of degrees of freedom associated to the reconstruction problem in online stage is same as the cardinality of basis Φ, and it influences the efficiency of online phase. This is the key to low-rank methods, as it is based on the assumption that the intrinsic structure of the system behaviour lies in a low-rank subspace. If this assumption does not hold, it may be difficult to generate a RB that approximates well the system behaviour in the entire parametric space.

Next, two well known methods of finding the low-rank subspace: the Proper Orthogonal Decomposition and the Reduced Basis Methods are described briefly, as these methods form the backbone of many reduced models [START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF][START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF]. This is followed by the necessity and tools for efficient computation of terms that are non-linear w.r.t. the parameter and/or the solution.

Proper Orthogonal Decomposition

In Proper Orthogonal Decomposition (POD) framework, high-fidelity snapshots are generated over the training set P tr are collected and arranged as columns of a matrix. Let S tr = {u s } N S s=1 ∈ R N h u ×N S be a matrix containing the training set snapshots. Here, N S is the number of snapshots and N h u is the number of finite element dofs. Typically, N S N h u for two reasons: first is that the snapshot computation can be expensive and second is the fact that discovery of a low rank subspace only needs a relatively small number of snapshots well distributed in the parametric domain. The RB can generated using the Singular Value Decomposition (SVD) or the Principal Component Analysis (PCA) methods. The SVD can be expressed as:

S tr = ΦΣΨ T
where Φ ∈ R N h u ×N S and Ψ ∈ R N S ×N S are unitary matrices containing orthonormal vectors. Σ ∈ R N S ×N S is a diagonal matrix containing singular values.

The decomposition given by the SVD has many interesting properties. The columns of Φ form a orthonormal basis for column-space of S tr and therefore can be used to approximate the solution subspace. Moreover, the columns of Φ, called as singular vectors, are arranged in the decreasing order of contribution to the total energy in S tr . Their contribution is indicated by the singular values in Σ, which can be used to truncate the basis Φ. In fact, the singular values {σ i } N S i=1 are arranged in decreasing order, where σ i is the i-th diagonal value of Σ, exhibit many interesting properties in relation to the matrix S tr , some of which are listed below:

• The square of singular values sum up to the total energy in the snapshot matrix, i.e. ||S tr || F = Ns i=1 σ 2 i

• The energy contribution by first n singular vectors to the snapshot matrix is given by n i=1 σ 2 i . • If Φ and Ψ contain the first n columns of Φ and Ψ, and also Σ ∈ R n×n contains the first n singular vectors {σ i } n i=1 in its diagonal; then the approximation

S tr ≈ Φ Σ Ψ T
is the best n-rank approximation in the sense of Frobenius-norm, with an approximation error

N S i=n+1 σ 2 i .
A related property is that Φ and Ψ form the best n-rank subspaces of the column and row space, respectively, of the matrix S tr in sense of Frobenius-norm.

The above properties can be used for truncation of the basis Φ by applying a tolerance δ to the fraction of energy contribution of first n singular vectors. The first n singular vectors of Φ that contribute more 1δ fraction of energy are kept. This number n can be calculated using:

min n 1 - n i=1 σ 2 i Ns i=1 σ 2 i < δ
The truncation leads to a subspace of the full solution space, whose rank is a function of tolerance δ. This allows us to regulate the dimensionality but at the cost of losing diminished accuracy. Generation of the reduced model using the reduced basis Φ will be shown in Section 1.3.1.

Reduced Basis Methods

Even though POD methods offer a robust technique to compute the RB, they do not offer a strategy to sample the parametric space. On the other hand, the Reduced Basis Methods (RBM) are usually equipped with error estimators that allow greedy sampling of the parametric space [START_REF] Prud'homme | Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods[END_REF][START_REF] Prud'homme | A Mathematical and Computational Framework for Reliable Real-Time Solution of Parametrized Partial Differential Equations[END_REF]. The idea of greedy methods is to arrive at the objective by using the locally best way forward at each step, which is typically done by solving a minimization problem. An error estimator ∆(µ) is an analytical error bound for the true reconstruction error as follows:

||u(µ) -Φ u(µ)|| 2 ≤ ∆(µ) ∀µ ∈ P
The error estimator is computed using the discrete residual r h ) of the high-fidelity system. For an unconstrained discrete problem with the residual:

r h (u, µ) = f (µ) -K(µ)u
the error bound ∆(µ) can be expressed as:

∆(µ) = K -1 (µ)r h (u, Φ u(µ)) 2 ≤ K -1 (µ) 2 ||r h (u, Φ u(µ))|| 2 = 1 σ min (K(µ)) ||r h (u, Φ u(µ))|| 2
Further details pertaining to efficient computation of the residual norm and the least singular value can be found in [START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF]Chapter 3]. The error estimator has negligible cost of evaluation compared to the true error, as the true numerical solution u is unknown. Therefore, at given greedy level k, the estimated best value of the parameter to be added is given by µ k+1 = arg max ∆(µ) and is cheaper to compute. Thus, the training step is enriched with a new snapshot u(µ k+1 ) and the RB Φ is updated using methods like Gram-Schmidt process. This avoids the computation of the basis from scratch at each step k. The process can be stopped when max ∆(µ) is under a certain tolerance. The RBMs are advantageous in terms of reliability and efficiency of the reduction process, though the performance of the process, of course, depends closeness of the error estimator with the true error (see [START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF]Chapter 7]).

Computation of non-linear terms

Operators such as K(µ) and f (µ) can be a non-linear function of µ and therefore, their computation in each iterative step of a linearized reduced problem can prove to be expensive. Despite the introduction of low dimensional spaces V r , the complexity of computing these operators depends on the number of dofs in high-fidelity model. Commonly used methods to efficiently compute non-linear terms include the computation of affine decomposition using the Empirical Interpolation Method (EIM) and the computation of an reduced integration domain (RID) using hyper reduction methods. The affine decompositions of the non-linear terms can be written as:

K(µ) ≈ N K i=1 α K i (µ)K i f (µ) ≈ N f i=1 α f i (µ)f i
where the scalar values α K i and α f i are the only quantities dependent on µ. Therefore, the cost of online computation for each query of µ is reduced to the cost of computing these scalar values.

It is possible that simple dependencies on µ admit affine decompositions that are easier to compute (see elastic energy computation in [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]). However, this is not the case in general and also there are many cases where operators show dependence on the solution u itself, for e.g. large deformation problems. In such cases, an affine decomposition is even more valuable as the resolution of the reduced problem will involve an iterative procedure demanding evaluation of non-linear terms in each iterative step. Usually, the EIM framework is used for computing the affine decomposition in such cases [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF]. EIM based affine decomposition is also equipped with so-called magic points in the non-parametric domain. Evaluation of the non-linear term only at the magic points is sufficient for computing the coefficients α K i and α f i . The cost of computation is reduced drastically as the parameter independent modes K i and f i can be computed offline. For further details on computation of EIM based affine forms, see [START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF]Chapter 10].

Similar reduction in computation cost can also be achieved by hyper-reduction methods, where the non-linear terms are evaluated only a selected subset of the quadrature points in the entire integration domain. Recalling the property of standard Galerkin approach in finite elements, each row in matrix K(µ) (and each entry in f (µ)) corresponds to integration over the local support of a particular node. The hyper-reduced evaluation of non-linear operators can be written as PK(µ) (and Pf (µ)). The rectangular matrix P each row containing only one non-zero entry equal to 1, at a position corresponding to a node in the RID [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF]. In other words, weak forms are computed only on supports surrounding the nodes in the RID. The RID has be constructed by selecting specific dofs or quadrature points in the entire domain. Various methods have been applied to this end, including free energy indicators [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF], optimized cubature methods [START_REF] Hernández | Dimensional hyper-reduction of nonlinear finite element models via empirical cubature[END_REF][START_REF]Reduced order and sparse representations for patient-specific modeling in computational surgery[END_REF], and Discrete Empirical Interpolation Method [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF].

ROMs for parametrized contact mechanics problems

ROMs have been applied extensively to unconstrained problems such as thermal and mechanics problems [START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF][START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF] and also to mixed problems with equality constraints such as the incompressibility condition in Stokes and Navier-Stokes flow problems [START_REF] Quarteroni | Numerical solution of parametrized navier-stokes equations by reduced basis methods[END_REF]. Application of ROMs to variational problems with inequality constraints (also referred to as variational inequalities) has been more recent [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF][START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF][START_REF] Haasdonk | A Reduced Basis Method for Parametrized Variational Inequalities[END_REF][START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF][START_REF] Bader | An empirical interpolation approach to reduced basis approximations for variational inequalities[END_REF][START_REF] Burkovska | Reduced basis methods for pricing options with the black-scholes and heston models[END_REF]. Inequality constraints appear in mechanical problems with obstacles or multi-body mechanical problems where there is a possibility of contact between bodies and obstacles, or with other bodies. Moreover, the region of contact is unknown a priori. Although the scope of this thesis is limited to contact mechanics, several other applications of variational inequalities are found in porous media flow problems [START_REF] Oden | Theory of variational inequalities with applications to problems of flow through porous media[END_REF], cavitation problems in lubrication systems [START_REF] Kinderlehrer | An Introduction to Variational Inequalities and Their Applications[END_REF], anti-plane frictional problems [START_REF] Sofonea | Variational Inequalities with Applications[END_REF] and even in financial trading problems [START_REF] Burkovska | Reduced basis methods for pricing options with the black-scholes and heston models[END_REF].

As seen in the Section 1.1.2, inequality constraint problems are often posed in mixed form with the KKT conditions that force the Lagrange Multipliers to be non-negative. Enforcing this condition in ROMs is non-trivial as traditional methods of computing RBs do not preserve the non-negative nature of the input information.

One of the first works on reducing the contact mechanics problem [START_REF] Haasdonk | A Reduced Basis Method for Parametrized Variational Inequalities[END_REF], proposed the idea of using the contact pressure snapshots directly to define a non-negative subcone. Both displacement and contact pressure snapshots were generated in a greedy fashion using error estimators. Fundamental aspects of the reduced problem, such as the existence and uniqueness of the solution and inf-sup stability were also explored.

Compression of snapshots to create a reduced basis was studied in [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF], where Non-Negative Matrix Factorization (NNMF) was used to compute a basis with user-specified cardinality but does not provide any means to specify truncation tolerance for the reduced basis. Also, error estimators for greedy sampling of parametric space are developed.

Algorithms to sort a precomputed set of snapshots in order of importance to create a compact basis are explored using projection methods in [START_REF] Burkovska | Reduced basis methods for pricing options with the black-scholes and heston models[END_REF] and [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF], using an Angle-Greedy and Cone-Projected Greedy (CPG) procedures respectively. The former does not take into account non-negative restrictions, whereas the latter uses a cone projection involving non-negative coefficients. This makes CPG more efficient in capturing the contact pressure subcone.

The hyper-reduction approach in [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF] defines a subdomain of the contact problem containing the most important points. The reduction is achieved with the usage of a POD basis for displacement and resolution of the weak form on the reduced integration domain. Reconstruction of displacement solution on the full domain is relatively straightforward using the POD basis. On the other hand, the reconstruction of contact pressure involves solving a non-negative least square problem using the snapshots.

A Craig-Bampton based resolution of the contact problem was discussed in [START_REF] Manvelyan | An efficient model order reduction scheme for dynamic contact in linear elasticity[END_REF] where the reduction of the displacement field was achieved using the Krylov subspace method. The Lagrange Multiplier method was not reduced under the assumption that the number of contact dofs remains small.

All of these works are based on Lagrange Multiplier approach to solve the inequality constrained reduced problem. An exception to this trend is [START_REF] Bader | An empirical interpolation approach to reduced basis approximations for variational inequalities[END_REF] where penalty approach is used. Instead, EIM is used for efficient computation of barrier functions that are used to restrict the solution in the feasible domain.

Most contributions build reduced models on problems with static contact pairs, meaning contact pairs do not change with the state of the system. However, [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF] considers dynamic contact pairs, with node-to-segment formulation. To efficiently evaluate non-linearities due to dynamic pairing, they use Empirical Interpolation Method to define an affine decomposition of distance functions.

Low-rank approach for contact problems

As seen above, many of the ROM approach to contact mechanics are based on the low-rank approach. As the validity of low-rank hypothesis will be discussed in detail in Chapter 2, a generic low-rank approach is detailed here, based on ideas of [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF][START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF]. The reduction of the displacement and contact pressure fields is discussed, followed by defining the reduced problem using the reduced spaces for the two fields. As the contact problem is a mixed problem, the displacement and contact pressure field will also be referred to as primal and dual fields Reduction of displacement: The displacement field u for a contact mechanics does not need any special treatment. Any methods that can be applied for reducing a mechanics problem without contacts can also be used for the displacement. Left singular vectors of displacement snapshots were used in many references in ROMs of contact mechanics. The same will be the case in this thesis unless specified, and the primal RB will be derived using POD.

Reduction of contact pressure: Computing the dual basis is more complicated than the primal basis. The dual field i.e. the Lagrange multiplier must satisfy the nonnegativity constraints. In case of the contact problems without adhesive and cohesive surfaces, this constraint admits a physical meaning as negative contact pressure cannot be admitted, as discussed in Section 1.1.3.2. An orthogonal basis generated using projection based methods cannot satisfy such constraints, as it is surely bound to contain negative entries. To ensure such constraints, one way is to define a subcone for the dual field instead of a subspace. A function cone, unlike a function space, is spanned by a set of non-negative basis functions and non-negative coefficients. In fact, W + and W h + in Section 1.1.3 are also cones. By extension, a subcone, unlike a subspace, must be equipped by a non-negative RB and must be spanned by non-negative coefficients.

W r + = span + (Θ) = N i=1 λ i θ i , λ i ≥ 0
where Θ and {θ i } N i=1 are the dual RB and dual RB functions. λ i s are the dual reduced dofs. Once the reduced subspace V r and subcone W r + are available, the reduced problem can be generated from (1.6), replacing the continuous spaces V and W + with the reduced spaces V r and W r + . The reduced KKT conditions on the solution ( u, λ) ∈ R N r u × R N r u can be expressed as:

K(µ) v -v T f (µ) + C T (µ, v) λ = 0 (1.14a) λ ≥ 0 (1.14b) C(µ, u) u -g(µ, u) ≤ 0 (1.14c) λ T ( C(µ, u) u -g(µ, u)) = 0 (1.14d)
where the discrete reduced operators are built by introducing the reduced basis functions φ and θ

K(µ) ij = a(µ; φ i , φ j ) C(µ, v) ij = b(µ, v r ; θ i , φ j ) f (µ) i = f (µ; φ i ) g(µ, v) i = d(µ, v r ; θ i ) with v r = Φ v
It is evident that the reduced form (1.14) has the same structure as the discrete form generated using finite elements in (1.9), as both equations are built using a Galerkin formulation and by introducing discrete solution spaces, albeit of different dimensionality.

Non-negative nature of the dual RB Θ forbids the use of orthogonal decompositions of Lagrange Multiplier snapshot matrix Λ. To this end, non-negativity preserving decompositions have been explored by various authors. Non-negative matrix factorization (NNMF) method [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] can be used to decompose a non-negative matrix, such as Λ into two low-rank non-negative matrices W, H, such that Λ ≈ WH. The left-hand matrix W is used as the dual RB in [START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF].

Another possibility is to directly use the snapshot vectors as the basis vectors, instead of computing the NNMF. This was the approach of [START_REF] Haasdonk | A Reduced Basis Method for Parametrized Variational Inequalities[END_REF], where the possibility of non-unique dual solution is also discussed, as snapshot vectors are not guaranteed to be linearly independent. Also, the number of dual dofs increases with the number of snapshots, preventing optimal reduction of the system. A greedy snapshot selection method based on the criteria of maximizing the volume of the reduced cone, namely the Cone-Projected Greedy (CPG) algorithm, was proposed in [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]. This algorithm creates a dual RB by greedily selecting snapshots from the snapshot matrix, attempting to create a more compact basis in comparison to using the full snapshot matrix. The algorithm is based on projection Π Θ of a vector λ on a vector cone Θ,

Π Θ (λ) :≈ Θα , where α is arg min γ ||λ -Θγ|| ∀γ ≥ 0 (1.15)
The greedy algorithm evaluates the error between each snapshot and its projection on the cone spanned by previously selected snapshots, and then adds that snapshot with maximum cone projection error. The process continues until the cone projection error is within a set tolerance. [START_REF] Burkovska | Reduced basis methods for pricing options with the black-scholes and heston models[END_REF] had also proposed a similar algorithm called Angle-Greedy algorithm, but the projection error is calculated based on computation of the angle between the candidate snapshot vector and the space, and not the cone, spanned by the previously selected snapshots i.e. it does not place non-negativity constrain on the coefficients γ in (1.15). This is not the best way of selecting snapshots for variational inequality problems since solutions must be sought in the reduced cone, and not the entire reduced space.

Remark 1.1 : Constraint satisfying primal subspace in flow problems

If we shift our gaze to other types of constrained ROMs, an interesting treatment of the mixed problems appears in incompressible flow problems. In such problems, the velocity snapshots are divergence-free and the divergence operator in the mass conservation equation is linear; offering the possibility to compute a divergence-free subspace for the velocity field. Consequently, the reduced incompressible flow problem is as an unconstrained problem, as all candidates in velocity subspace naturally satisfy the constraints [START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible navier-stokes equations: rigorous reduced-basis a posteriori error bounds[END_REF][START_REF] Liberge | Reduced order modelling method via proper orthogonal decomposition (pod) for flow around an oscillating cylinder[END_REF]. Apart from naturally satisfying the constraints, some elegant properties of vector algebra a result in disappearance of bilinear term that includes pressure. Unfortunately, it is not straightforward to extend the same idea to contact mechanics problems (or variational inequalities), i.e. to compute a primal subspace that satisfies the non-penetration (or the inequality) constraints. This is because linear combinations of non-penetrating snapshots do not satisfy the non-penetration condition in general. As a consequence, the inequality constraints need be enforced explicitly while solving the reduced problem, in general.

a The product of a divergence-free velocity field and the gradient of pressure vanishes due to the orthogonality of the solenoidal and irrotational vector fields Construction of non-linear operators: For an efficient ROM, the construction of nonlinear operators must also be cheap. The construction of operators C(µ, u) and g(µ, u) in each iterative step can prove expensive. Affine decompositions of the nonlinear constraint operators using the EIM are discussed in [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]. The affine form was achieved by computing the EIM based interpolation of the distance functions k(µ, v, •)(x) and g(µ, v)(x) (recall Notation 1.1):

k(µ, v; •)(x) ≈ N k p=1 α k p (µ, v)κ p (•)(x) g(µ, v)(x) ≈ Ng q=1 α g q (µ, v)γ q (x)
where the interpolation is exact at the set of magic points {x k i } and {x g i }. Note that κ p is linear w.r.t. its argument which can be any member in V. Therefore, the coefficients can be computed so that they satisfy the following conditions:

k(µ, v; •)(x k m ) = N k i=1 α k p (µ, v)κ p (•)(x k m ) ∀m = {1, 2, . . . N k } g(µ, v)(x g n ) = Ng q=1 α g q (µ, v)γ q (x k n ) ∀n = {1, 2, . . . N g } (1.16)
Therefore, the distance function needs to be evaluated at most N k + N g points, instead of the whole contact surface. Also, κ i (•) and γ i are computed just once using the training set and the corresponding discrete operators can also be built during the offline stage using the definition of weak forms in (1.5).

C p ij = Γ 2 θ i κ p (φ j )(x) ∂Γ g q i = Γ 2 θ i γ q (x) ∂Γ (1.17)
where b p and d q correspond to the weak forms of individual EIM modes κ p and γ q that are independent of the state of the system. Thus, during the online stage, the non-linear operators can be approximated as:

C(µ, v) = N k p=1 α k p (µ, v) C p g(µ, v) = Ng q=1 α g q (µ, v) g q
where only α k p and α g q are evaluated efficiently in online stage using (1.16). Another approach for efficient computation of non-linear operators is the hyperreduction approach [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF]. In this approach, a small subset of the full domain Ω A ⊂ Ω 1 ∪ Ω 2 is computed that serves as the reduced integration domain (RID). The RID Ω A consists of the most "important" integration points computed by performing EIM on the displacement RB. Then, integration of the operators C(µ, v) and g(µ, v) is carried out on a subset of the potential contact surface, computed using the intersection of the RID and the contact surface Ω A ∩ Γ 2 . Thus, the distance functions are computed only on a small subset of the contact surface.

Scope of the thesis

The current literature consists of various methodologies to reduce the contact mechanics problem, especially to the Lagrange Multiplier approach. In this thesis, the following contributions are made:

• As seen in this chapter, the reducibility of contact problems has not been sufficiently investigated in current literature. This is particularly important as contact mechanics problems exhibit a local nature in terms of the contact zone. As the region of contact is dynamic for many contact problems, significant changes in the contact zone can prove consequential to reducibility. Therefore, in Chapter 2, we explore the validity of so-called low-rank hypothesis and robustness of the reduced models using various metrics. The chapter concludes that the contact pressure field lacks linear separability, an essential feature for reducibility.

• The lack of low-rank structure implies that the offline phase might be expensive as small training sets will be unable to explore the underlying structure sufficiently.

Even if the user is willing to compute large training sets, it is still not possible to compute a robust low-rank subspace for the contact pressure. In Chapter 3, we focus on algorithms to handle dictionaries containing a large number of snapshots using ideas from sparse regression methods.

• As the linear inseparability of contact pressure field is established, we explore the idea of non-linear interpolations between contact pressure snapshots in Chapter 4.

Interpolation is done in a warped space generated using Dynamic Time Warping algorithm. We show that reconstructions of contact pressure can be more reliably performed with smaller training sets, demonstrating the potential of non-linear dimensionality reduction methods for contact mechanics problems.

| Limitations of low-rank approach to contact problems

Strategies for reducing complex numerical models to facilitate building parametric models requires the acquisition of underlying physics using precomputed data. A brief overview of these methods was given in Section 1.2. Commonly used ROM methods are based on the paradigm of the so-called linear dimensionality reduction. The foundation of this framework rests on the low-rank hypothesis, which states that for a high-dimensional quantity of interest, often a representative low-dimensional space exists and can be extracted using using linear transformations. The successful applications of ROM methods like POD, based on linear dimensionality reduction methods like PCA, to many parametric models can be attributed to the underlying low-dimensional behaviour of these systems. In Sections 1.2 and 1.3, the wide variety of low-rank applications were mentioned, including constrained problems such as incompressible flows and contact mechanics problems. Reduced models for contact mechanics are important because the inequality constraints are geometrically non-linear, thereby introducing additional computational challenges (see Sections 1.1.2 to 1.1.3). However, the assumption of low-rank behaviour has implications on reduced models of contact mechanics problems, which are a class of variational inequality problems.

This chapter focuses on applicability of the low-rank hypothesis to the variational inequality problems, specifically contact mechanics problems. The effect of this hypothesis on various contact mechanics problems are demonstrated. A generic methodology of low-rank approach in contact mechanics problems was discussed in Section 1.3. In this chapter, the limitation of low-rank approach will be demonstrated on reduced models of various contact problems. It will be shown that the contact pressure field is linearly inseparable and therefore does not satisfy low-rank hypothesis. To endorse this argument quantitatively, the so-called validation metrics will be introduced to used to endorse the arguments about linear inseparability using specific numerical examples. Finally, con-cluding arguments about linear inseparability and a few perspectives on circumventing the limitations due to lack of low-rank structure are given.

Low-rank contact problem using active set method

The non-linear reduced problem in (1.14) is based on the low-rank hypothesis, as it introduces the low-rank primal and dual bases. This approach is similar to [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF][START_REF] Haasdonk | A Reduced Basis Method for Parametrized Variational Inequalities[END_REF][START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF]. Two sources of non-linearities exist in this problem: the resolution of active/inactive constraints and the dynamic contact pairs. Irrespective of the linearization algorithm utilized, both non-linearities depend on the current estimate of the displacement field at a given linearization step.

Here, (1.14) is solved using the fixed-point iterations and the active-set method applied to the reduced constraints. The linearized problem at a fixed-point iteration level p can be expressed as follows:

K C T A (µ, u p ) C A (µ, u p ) 0 u p+1 λ p+1 A = f g A (µ, u p ) (2.

1)

The subscript A indicates the active set of reduced contact constraints, i.e. the constraints that satisfy the equality and force the solution to lie on the boundary of the feasible region. Note that (2.1) is similar to the fixed-point system in (1.12). The fixed point algorithm stated in Alg. 2.1, was implemented in Python using NumPy package [START_REF] Harris | Array programming with NumPy[END_REF] for vectorial operations.

Construction of operators:

The implication of dynamic contact pairing between two bodies of a contact problem is that the non-linear operators that enforce the nonpenetration constraint must be updated for each new estimate of displacement field. On the other hand, for generating an efficient ROM, the construction of nonlinear operators must be inexpensive. This is not the case in Alg. 2.1, where the full order operators C(µ, u) and g(µ, u) are constructed in each iteration (Step 5). Affine decompositions of the nonlinear constraint operators using the EIM are discussed in [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]. EIM decomposition of the constraint operators further splits computational complexity into offline and online stages, by selecting the so-called "magic points" that are a small subset of all potential contact nodes where the nonlinear terms are evaluated in the online phase. The hyper-reduction method also permits computation of operators using fewer points [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF]. This is crucial to solve problems where the contact pairs are strongly dependent on the displacement field. However, the efficient construction of constraint operators is not covered in this work, and the focus is maintained on the non-negative dual space and its low-rank approximation. Build constraint operators C(µ, u p ), g(µ, u p ) using FEM

6:

Project constraint operators on RBs C(µ,

u p ) = Θ T C(µ, u p )Φ g(µ, u p ) = Θ T g(µ, u p ) 7:
Filter-out rows that are not in active set:

C A (µ, u p ) = C(µ, u p )[A, :] g A (µ, u p ) = g(µ, u p )[A] 8:
Solve system (2.1)

9: Set λ p+1 [A] ← λ p+1 A and λ p+1 [A c ] ← 0 10:
Update active constraints set

A[i] = False if λ p+1 i < 0 True if ( C u p+1 -g ≥ 0) i violations of (1.14b)
and (1.14c) 11: end while 12: Reconstruct u = Φ u and λ = Θ λ 13: Output: u, λ

Illustrative case

In this section, the Type-1 problem of a 1D elastic rope-obstacle problem from [START_REF] Haasdonk | A Reduced Basis Method for Parametrized Variational Inequalities[END_REF][START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF][START_REF] Bader | An empirical interpolation approach to reduced basis approximations for variational inequalities[END_REF] is considered. The geometry of the problem enjoys further simplification of static contact pairs. In other words, each point on the elastic rope can come in contact with a unique point on the obstacle, thereby stripping the distance functions k and g of their non-linear dependence on the displacement field u.

The parametrized model of the elastic rope with obstacle function can be expressed as:

ν∇ 2 u(x) = f on x ∈ [0, 1] u(0) = u(1) = 0 u(x) ≥ g(x, γ) on γ ∈ [-0.5, 0.5] (2.2)
where the parameterized obstacle function is defined as:

g(x, γ) = -0.2(sin(πx) -sin(3πx)) -0.5 + γx (2.3)
and the quantities ν = 30 and f = 250 are independent of the parameter γ. Snapshots are generated in the training set, P tr ∈ P consisting of 10 equidistant points in the parametric space P defined in (2.3). The deformation and contact pressure snapshots are shown in Fig. 2.1. The primal RB is built using the POD approach and dual RB is created using the cone-projected greedy algorithm of [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]. The reconstruction errors are shown in Fig. 2.2 for the two cases of full and truncated dual RB. As expected, the points in the training set show very low reconstruction error, and the points outside the training set show a moderate error. The only exception at γ = -0.35 where reconstruction with full RB is highly accurate because the contact area does not change in the regime γ ∈ [-0.4, -0.3] and the training set has two snapshots in this region. When dual RB is truncated, some points in the training set achieve the same error level as that of points outside the training set. These points correspond to the same dual snapshots that were eliminated by the truncation procedure of the CPG algorithm. Another relevant observation is that the number of active dofs of the reduced dual solution λ is small for all reconstruction cases. For points in the training set, λ has exactly 1 active dof, whereas for the points outside the training set, it has a maximum of 2 active dofs. This is due to the fact that the dual basis is composed primarily of snapshots and does not undergo any compression like in the case of the primal basis generated by POD. In general, any efficiently generated reduced basis has its vectors arranged in order of importance. In case the basis consists of left singular vectors, generated by SVD, the decreasing importance of basis vectors is indicated by the singular values. This is a very useful indication of the rank of the subspace. Even though the dual basis is not computed using SVD, it can still be useful to compute singular values and assess their decay. This evolution of singular values, seen in Fig. 2.3, shows that the decay is slow compared to the primal variable. Also, the discussions in previous sections highlight the fact that dual basis Θ cannot be truncated without a significant loss of accuracy. This is also evident from Table 2.1, even a truncation tolerance as high as 0.2 leads to truncation of just four vectors in dual RB.

The above observations are explainable by the hypothesis that the dual snapshots are highly "inseparable" and do not lie in a low-dimensional subcone. Each contact pressure snapshot in Fig. 2.1 is linearly independent of the other due to difference in the contact area. Hence, if any snapshot is picked randomly from this figure, it cannot be reasonably approximated by any linear combination of the rest of the snapshots. In other words, each snapshot vanishes outside the contact zone, and this zone is almost "unique" for each snapshot. For a particular snapshot to be reasonably approximated, other snapshots with similar contact zones are needed. Therefore, the elimination of snapshots from the dual basis leads to a drastic increase in the reconstruction error for the corresponding parametric values in the training set. The dual variable, therefore, not only needs special treatment due to its positivity constraints but is also highly inseparable because of sensitivity to contact position and area. The inseparability is also evident in the non-zero pattern ("spy" plot) of dual snapshots visualized in Fig. 2.4, where most rows have an almost unique sparsity pattern. It is unlikely that a subspace whose members show varying sparsity patterns will admit a low-dimensional behaviour. More numerical examples are explored to study the lack of low-rank behaviour in the Section 2.4.

Metrics for low-rank methods

In this section, we introduce the validation metrics [START_REF] Davies | Learning shape: optimal models for analysing natural variability[END_REF], that provide a quantitative evaluation of the robustness of low-rank models. Validation metrics defined in [START_REF] Davies | Learning shape: optimal models for analysing natural variability[END_REF] identify three different measures for this purpose:

• Compactness: As the name suggests, this metric is a measure of dimensionality of the reduced space. Reference [START_REF] Davies | Learning shape: optimal models for analysing natural variability[END_REF] defines compactness as the squared sum of the first m singular values of the snapshot matrix. This metric is similar to the decay of singular values studied in the previous section in Fig. 2 The Frobenius norm is chosen in (2.4), because it is a natural choice for a POD/SVD basis that contains the best rank-wise approximations in the sense of Frobenius norm, but other norms can be also be used. This is actually equivalent to the original definition of compactness in [START_REF] Davies | Learning shape: optimal models for analysing natural variability[END_REF], because in the case Ψ consists of the left singular vectors, compactness corresponds to cumulative energy of singular values.

For the dual field, the metric is given using two kinds of projections: orthogonal and CPG. Though the dual CPG metric is more relevant to the inequality constrained problems, the dual orthogonal metric is also included as it is more intuitive due to correspondence with the singular values.

• Generalization Ability: This metric measures the ability of the reduced basis to approximate parametric instances that are not in the training set. It is computed by using a leave-one-out approach on the training set to generate a m-rank basis. Then, the reconstruction error of the eliminated snapshot is used to compute generalization ability as given by (2.5),

G(m) = 1 N N i=1 s i -Π Ψ(Str\s i )[:m] (s i ) ||s i || (2.5)
where Ψ(S tr \s i ) is the basis created using the snapshot subset S tr with s i removed (i.e. S tr \ s i ). The symbols S tr , Π and Ψ follow the definition in (2.4).

• Specificity: The subspace spanned by the reduced basis is expected to contain elements similar to that of the training set. This metric measures the extent of dissimilarity between randomly picked elements from the subspace and the snapshots in the training set. To compute specificity, a set of random shapes in the associated subspace are generated using random coefficients R = {Ψα i } N i=1 , where each vector α i is drawn randomly within a predefined range, and compared to the closest snapshot in the training set. The mean value of this error is defined as Specificity.

S(m) = 1 N N i=1 min s j ∈Str ||s j -Ψ[: m]α i || ||s j || (2.6)
The vectors {α i | α i ∈ R m } are drawn randomly in the range defined by µ α ±σ α ⊂ R m . The vectors µ α and σ α contain the element-wise mean and standard deviation of vectors {γ k } m k=1 , where γ k is the reduced coordinate of the training set snapshot k. (e.g. if Ψ is an orthogonal basis, the vectors γ k are the columns of matrix Γ := Ψ T S tr ). In the specific case where the basis Ψ is the Θ, the CPG basis, only non-negative entries are allowed in {α i }. The symbols S tr , Ψ follows the definition in (2.4). High-fidelity model: A finite element model of the geometry in Fig. 2.5a is created using a quad mesh. Each half-cylinder has 513 nodes and 466 elements of which 78 elements lie of the potential contact surface i.e. the semi-circular edges. Linear shape functions are used to discretize the displacement field, while contact pressure is discretized using piece-wise constant shape functions centred at the surface nodes of the quad mesh (collocation method, as described in [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]). For surface integrals, a single-point gauss quadrature centred on the node is used. This formulation gives mostly smooth contact pressure profiles, except near the peak pressure.

Numerical Examples

Reduced Model: Snapshots are generated in the training set, P tr ∈ P consisting of 12 equidistant points in the parametric space. The full rank and truncated dual RB (using CPG) is considered, as given in Table 2.2, whereas the primal RB is not truncated. The full RB has 12 basis vectors for both primal and dual fields, whereas the truncated dual RB has 8 vectors. Note that dual RB truncates only 4 out of 12 basis vectors for a high truncation tolerance of 0.05.

The reconstruction errors are shown in Fig. 2.6 for full and truncated dual RBs. As expected, the points in the training set show very low reconstruction error, and the points outside the training set show a relatively high error. When the dual RB is truncated, points corresponding to truncated dual snapshots have a high error level, in the same order of error for points outside the training set. Like the illustrative case in Section 2.2, the non-zero pattern in snapshots in Fig. 2.7 also show a unique sparsity pattern for each snapshot. Also, the comptactness shown in Fig. 2.8 displays a slow decay of the dual orth and dual cone, indicating a high-rank behaviour of the dual variable. The slightly slower dual cone curve than the dual orth was expected because of the nonnegativity constraint that appears in approximation using the cone of the dual basis. These observations again reinforce the proposition that the dual basis does not admit a low-rank behaviour. Generalization ability and Specificity metrics introduced in Section 2.3 are plotted in Fig. 2.9. The generalization ability for primal field is much better than the dual field, another reflection of separability issues with the dual field. On the other hand, primal and dual specificity are about the same order, though the primal specificity is consistently lower.

Ironing problem

A Type-2 problem where the contact area changes quite significantly, namely the ironing problem [START_REF] Fischer | Mortar based frictional contact formulation for higher order interpolations using the moving friction cone[END_REF], is considered. The ironing problem consists of an iron block pressed against a flat slab and moved along the length of the slab (Fig. 2.10). The problem is simplified with two more assumptions: the first is that iron moves slowly enough that the problem can be considered to be quasi-static and the second is that surfaces are frictionless. The horizontal position d x of the iron is taken as the parameter for the reduced model. In this problem, as the potential contact surface on the slab is quite larger than the actual position of contact, the contact pressure snapshots display large changes in contact position, shown in Fig. 2.11. The same is reflected in spy pattern of the contact pressure snapshot matrix, shown in Fig. 2.12 Details of FE model: A finite element model of the ironing problem is created using a structured quad mesh for both the iron and the slab. Like the Hertz problem, the displacement field is approximated using linear shape functions and contact pressure using piecewise-constant shape functions centred at the surface nodes. To better demonstrate the inseparability issues, two meshes: coarse and fine, are considered for iron and slab.

Iron Slab

Coarse Mesh 30 × 30 20 × 100 Fine Mesh 60 × 60 40 × 200

The compactness metric, which is equivalent to decay of truncation error as discussed before, is shown in Fig. 2.13 for a coarse and fine mesh, created using 128 snapshots. The contrasting part of the two graphs is the green curve representing the dual orth compactness, as this curve vanishes at rank 100 for the coarse mesh, but not for the fine mesh. This is natural, since coarse mesh has only 100 nodes on the surface of the slab and thus allowing a maximum rank of 100, whereas the fine mesh has around 200 nodes. Since the orthogonally computed dual (dual orth) compactness is also a measure for determining the true rank of the snapshot matrix, the green curve naturally vanishes at a rank of 100. However, the most interesting part of this figure is the blue curve representing the truncation computed using cone-projection. Even for the coarse mesh, whose true rank is bounded by the number of surface nodes, i.e. 100, the dual cone compactness doesn't decay well beyond this rank. This is another interesting behaviour of the low-rank approach since the dual solution is sought in a subcone that involves positivity constraints, further degrading compactness.

Apart from compactness, another way to demonstrate the dimensionality of the dual subcone is by using the projection error of nested training sets, which will be referred to as nested error subsequently. To compute the nested error, snapshots are computed in a nested set of points in the parametric space, so that n-th level training set is a subset of the (n + 1)-th level training set. In each nested level, 2 n + 1 points are uniformly distributed in the parametric space [0, L]. This means, every alternate point in (n+1)-th level is the mid-point of two consecutive points in n-th level. The nested error H(n) can then be defined as:

H(n) = S n+1 tr -Π Ψ n (S n+1 tr ) F S n+1 tr F (2.7)
where the symbols Ψ and S tr carry same meaning defined in Section 2.3. ||•|| F indicates the Frobenius norm. The RB computed using the n-th level training set is indicated by Ψ n , and snapshot matrix for (n + 1)-th level is indicated by S n+1 tr . Nested training sets where each level contains nearly twice the number of points compared to previous level are created and the dual RB for each level is computed. The slope of nested errors computed using orthogonal and cone projections indicate that the offline stage will require a very large number of snapshots to explore the dual subcone, which is another indication of its high-dimensionality, as shown in Fig. 2.14. The effects of mesh size are similar to those seen in compactness, as nested error using dual cone projection does not decay even after training set is full rank in case of coarse mesh. Generalization ability and Specificity metrics are plotted in Fig. 2.15. The generalization ability shows a similar trend as that of nested error, as they are similar quantities measuring reconstruction errors outside the training set. The primal basis outperforms its dual counterparts in this metric, as expected. The specificity metric is quite higher for dual quantities than the primal ones, which is because random shapes generated using the RB are quite different.

Conclusions and Perspectives

In this chapter, the limitations of low-rank ROMs to contact mechanics problems were discussed. Demonstrations of this approach show moderate reconstruction errors for parametric points outside the training set. To explore the possible sources of this error, we focus on the key assumption of the low-rank hypothesis i.e. the solution subspace is low-dimensional. It turns out that the Lagrange Multiplier associated with the inequality constraints, i.e. the contact pressure, does not admit a low-rank subspace. Qualitative and quantitative assessments are provided to support this argument. The linear inseparability of dual subspace/subcone is demonstrated using various metrics that show high-dimensionality of contact pressure compared to the relatively low-dimensionality of the displacement field.

To circumvent the limitations due to the lack of low-rank structure, several potential approaches need to be explored. One possibility is the usage of dictionary-based approximation with a relatively large number of snapshots spanning the dual subspace efficiently. Another possibility is resorting to either a non-linear transformation; potentially transforming the dual snapshots in a way that they lie in a low-dimensional subspace. contact mechanics problems

Motivation

The contact mechanics problem has a feature, the contact zone, that exhibits a local nature. The contact zone is a function of the state of the system, and hence, may strongly depend on the loading, geometry, and other physical parameters. In general, for variational inequalities, only a subset of the constraints are usually active. An inequality constraint is referred to as active only if the current state of the problem satisfies the equality. Due to the last KKT condition in (1.14d), the Lagrange multipliers λ assume a non-zero value only for active constraints i.e. the part of surface where contact is established. The strong dependence of the active zone on the state of the system leads to the linear inseparability of the Lagrange multiplier field. It has been shown in Chapter 2, that the local characteristics of contact pressure can render subspace learning ineffective and cause difficulties in creating efficient reduced models of contact problems. The separability issues may cause the following problems:

1. Complexity of subspace learning: It is desirable to keep the complexity of the offline stage at affordable levels. Usually, this is achievable because a relatively low number of snapshots compared to the dofs of the high-fidelity problem are sufficient to create a reliable reduced basis, thanks to the low-rank nature of many resolution problems. However, given the inseparability of the contact pressure field, a large training set may be necessary to discover the entire solution subspace to a reasonable accuracy.

2. Lack of low-rank subspace: Discovery of a low-rank subspace is key to creating a low-cost reduced model, as the number of reduced dofs scales with the rank of the discovered solution subspace. While low-rank methods for contact problems [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF][START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF] define a contact pressure subcone, it is shown in Chapter 2 using various metrics that compression of contact pressure data into a low-rank subcone or subspace leads to loss of useful information. Therefore, despite the creation of a large training set, as mentioned in the previous point, extracting a low-rank subspace might not be possible.

The degree of these effects depends on the nature of the contact problem. In this chapter, the approximation using a dictionary instead of a low-rank subcone is explored, addressing the second point as the dictionary is made of snapshots without applying any compression, like NNMF or CPG. In effect, the high complexity of the offline stage is accepted as a compromise for creating a reliable reduced model. The ideas behind over-complete dictionaries are expanded in the coming sections.

Sparse methods with over-complete dictionaries

The primary motivation for using a dictionary-based approximation is to mitigate the limitations due to linear inseparability, discussed in Section 3.1. The word "mitigate" is explicitly used here to indicate that the strategy developed in this chapter does not resolve the inseparability issue at a fundamental level, rather it attempts to limit the effects of inseparability. The intention is to accommodate the cost of computing snapshots over a large training set and focus mainly on addressing the point 2 of Section 3.1 using an algorithm that can handle large dictionaries efficiently. A matrix containing the large set of snapshots is referred to as an over-complete dictionary D. The term "over-complete" is used loosely in this context, as it is impossible to sufficiently sample a truly inseparable subspace for a given approximation error target and is simply used to imply that the dictionary has numerous snapshots, possibly larger than the number of contact dofs in the high-fidelity problem.

The dictionary-based approximation can be visualized in an abstract n-dimensional space in Fig. 3.1. Typically, in cases where low-rank hypothesis is not satisfied, the underlying physics of the problem assume a manifold structure. A manifold can be thought of as a low-dimensional entity that does not satisfy the properties of a linear space. The idea expressed in this figure is that a manifold, like a hypersurface with non-zero curvature, can be approximated better using a large cloud of points lying on it. If the point cloud is rich enough, it is more likely that this approximation will be better than than low-dimensional subspace which is a hyperplane1 . It is this large cloud of points that is referred to as over-complete dictionary in this chapter.

In dictionary-based approximation, unlike the low-rank approach [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF][START_REF] Balajewicz | Projection-based model reduction for contact problems[END_REF], no compression of the contact pressure snapshots is applied. The dual solution can simply be restricted to the column-cone2 of the dictionary D, which might be reasonably rich to approximate the true contact pressure subcone. In the online phase, only a few snapshots will be selected to estimate the contact pressure field. This is justified because not all information in the dictionary is necessary for the reconstruction of a particular instance and hence, the dual reduced dofs will admit a sparse solution. This feature is also seen Illustration of low-rank vs. dictionary-based approximation. Snapshots (circles) in an n-dimensional system lie on a manifold of dimension lower than n. The low-dimensional subspace computed using low-rank methods has a dimension n. A dictionary-based approximation of the queried point (star-shaped) can be computed as linear combination of the nearest snapshots (shown by a line connecting the nearest candidates). The low-rank approximation is simply the projection of the queried point on the low-dimensional subspace.

in case of the low-rank model generated using snapshots selected by CPG. Therefore, enough motivation exists to use sparsity enforcing methods to choose a few columns of a dictionary D to approximate the dual solution. Keeping these ideas in mind, sparse regression techniques are briefly introduced at this stage and a few of these methods from current literature are discussed. Following this, their application to resolution problems will be briefly discussed.

Sparse regression in data approximation problems

Sparse regression techniques have been under development in recent decades across various fields, especially in signal processing and statistics. The main idea behind these techniques is to approximate a signal with the few most suitable elements of a precomputed dictionary of signals [START_REF] Elad | Sparse and Redundant Representations[END_REF], or to develop the simplest predictive model [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF]. A typical sparse regression can be written as the following optimization statement:

min||α|| p s.t.||x -Dα|| < ε (3.1)
where x is the signal being approximated, D could either be a dictionary of signals or a basis that possibly admits sparse representation, ||•|| p is a sparsity inducing norm, and ε is the tolerance on the approximation error. It is also possible to write other forms of sparse regression problems where the approximation error is minimized and p-norm term is constrained to stay under a specific tolerance of sparsity, or the minimization of a penalty form with a weighted sum of the two terms.

It is common knowledge in the communities using sparse regression methods that the 0 -norm is a measure of cardinality, hence (3.1) with p = 0 is the natural choice in a theoretical sense. However, this problem is NP-hard in general, and therefore, it is more common to use p -norm with p ∈]0, 1]. A detailed explanation of sparsity inducing properties of p -norms with p ∈]0, 1] is given in [START_REF] Mairal | Sparse modeling for image and vision processing[END_REF]. A few sparse regression methods are briefly detailed here:

1. LASSO [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF]: The specific case of (3.1) with p = 1 is known as the LASSO is one of the most widely known sparse regression methods. The 1 -norm, even though not a true measure of cardinality, has a constant "driving force" or more technically speaking, a gradient towards the origin that is constant in magnitude everywhere. This drives many entries of α to zero, unlike the 2 norm that can reduce the magnitude of α entries, but not set them to zero. Also, the 1 -norm allows solving (3.1) in polynomial time.

2. Dantzig Selector [START_REF] Candes | The Dantzig selector: Statistical estimation when p is much larger than n[END_REF]: This method allows selection of a sparse α using a modified version of (3.1), that can be converted into a linear programming problem after with algebraic manipulations. This is achieved with p = 1 and the ∞ norm for the residual. (3.1) is recast as:

min||α|| 1 s.t. D T (x -Dα) ∞ < ε
which can be expressed in a linear programming form as:

min i γ i s.t. -γ < α < γ -ε1 < D T (x -Dα) < ε1
where γ is a dummy variable obtained by opening the absolute operator α .

3. Orthogonal Matching Pursuit [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF]: The OMP method aims to approximate the solution of (3.1) with p = 0 using a greedy technique. At each greedy step k, the OMP algorithm searches for the column of dictionary D that is the least orthogonal to the current residual (x -Dα k-1 ) and adds to the set of indices previously selected. This process is carried out until the residual is below the tolerance ε. This is done by projecting the residual vector on each column and selecting the largest projection.

I ← I ∪ arg max j D[:, j] T (x -Dα k-1 )
where I is the current set of selected indices and columns of dictionary D are normalized. Then the new set of coefficients α k are computed using the least square solution:

α k I = (D T I D I ) -1 D T I x with α I = α[I] and D I = D[:, I].
Update of I is locally optimal, but the above least-square update of α k , the OMP solution is optimal w.r.t to the currently selected subset D I of the dictionary. This is in contrast to the predecessor of OMP, the Matching Pursuit [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF], which used locally optimal updates for α k at the current iteration.

FOCUSS [70]:

The FOCUSS method finds a sparse approximation using an iterative process that starts from a fully dense solution and progresses towards "localized energy solutions". The FOCUSS method is broadly based on the following ideas:

• For a full rank dictionary of size, say m × n (m < n), the system Dα = x is under-determined. For such a system, the closed form solution for (3.1) with p = 2 is given by α = D † x, where D † is the Moore-Penrose pseudo-inverse. Note that this solution is not sparse in general, as 2 -norm does not have sparsity inducing properties.

• A modified minimization problem can be defined by replacing the ||α|| 2 with a weighted norm W -1 α 2 , or more generically W † α 2 if W is singular. The authors of [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm[END_REF] argue that by changing W, every possible solution of the under-determined system Dα = x can be obtained.

• The basis of the algorithm lies in iterating towards a sparse solution using a weight W k which induces sparsity. For the iteration k, the weight is chosen as W k = diag(α k-1 ). The trick lies in the fact that algorithm ends up minimizing the following weighted norm:

W † k α 2 2 = i α k i α k-1 i 2
due to which smaller entries of α k-1 tend to diminish further.

On combining these ideas, the FOCUSS algorithm ends up with quite a simple implementation, given in Alg. 3.1. Usually, the algorithm is initialized with the α 0 containing all non-zeros. In [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm[END_REF], the ||α|| 2 minimizing solution α 0 = D † x was used for initialization.

W k = diag(α k-1 )
5:

α k = W k (DW k ) † x 6:
k ← k + 1 7: end while

Application of sparse methods to unconstrained resolution problems

Sparse regression methods, though originally developed for application on pure data approximation or data compression problems, have seen recent interest in applications towards physics problems [START_REF] Brunton | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF][START_REF] Wang | Predicting catastrophes in nonlinear dynamical systems by compressive sensing[END_REF][START_REF] Schaeffer | Sparse dynamics for partial differential equations[END_REF]. The methods that have been discussed previously are also meant for data approximation. At this point, application of such methods on unconstrained resolution problems is discussed. A resolution problem can be framed as residual minimization problem in order to use dictionary-based approximation and sparse regression. Consider a discretized resolution problem with the state vector v and residual r(v). The basic idea is to approximate the solution v ≈ D v using the elements of dictionary such that the residual of the numerical model r(D v) is below a certain tolerance, while minimizing the p -norm of coefficient vector, as shown in (3.2). Sparsity is imposed on the coefficients of the dictionary, meaning that only a few dictionary columns are selected.

min|| v|| p s.t. ||r(D v)|| < ε (3.2)
Dictionary-based sparse approximation for resolution problems has been discussed in [START_REF] Balabanov | Randomized linear algebra for model reduction-part II: minimal residual methods and dictionary-based approximation[END_REF], where the solution space is set to the span of dictionary columns. They also discuss the application of random sketching to residual evaluation, which provides an efficient way of solving the minimization problem. This is similar to the approach of randomized SVD proposed in [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. The idea is to project the residual on a low-dimensional subspace in colsp(D). 3 The matrix B, containing a basis corresponding to a random combination of dictionary columns, is used for this projection in (3.3). This essentially reduces the complexity of enforcing the residual constraint. Given a dictionary D of size m × n, a smaller matrix that has information from all columns of D can be generated by multiplying a matrix R of size n × L with L min (m, n). After orthogonalization, this matrix B can be used for projecting the residual r.

min|| v|| p (3.3a) s.t. B T r(D v) < ε (3.3b) where B = orth(DR) , R ∈ U n×L [0,1] (3.3c)
The number L can be tuned to not only control the complexity of the minimization problem, but also to crudely constrain sparsity of the solution. For L n, colsp(B) ⊂ colsp(D) is low-dimensional and hence, the projection B T r is a "weaker" evaluation of r, which allows for a sparser v at the price of permitting residuals orthogonal to B.

Another related approach in which snapshots are directly used for reducing the model without any explicit compression is the CUR decomposition [START_REF] Le Guennec | A parametric and non-intrusive reduced order model of car crash simulation[END_REF]. In this method, not only columns, but also rows from the snapshot matrix are used for developing a parametric regression. It works by algorithmically selecting a column and row subset from the snapshot matrix, and using them for developing a regression model. A contrasting feature of the CUR approach is the use of regression without invoking the partial differential equation or its weak form, resulting in a non-intrusive method.

Diffusion-reaction problem

A transport problem of a quantity v with two-dimensional domain of size 3 × 3 units, divided into a grid with cells of size 1×1 is considered. The cells have a piecewise uniform diffusivity value. A three-dimensional parametrization is considered with diffusivity values ν 1 , ν 2 and ν 3 , as shown in Fig. 3.2. The reaction coefficient σ is considered to be uniform across all sectors.

-ν(x, y)∇ 2 v(x, y)

+ σv(x, y) = f v(x, 0) = 0 v(x, 3) = 0 v(0, y) = 0 v(3, y) = 1 (3.4)
where ν(x, y) describes the piecewise uniform diffusivity values.

The 2D domain is discretized into a structured 15 × 15 mesh. The diffusivity of each block ν i ∈ (1, 150) is considered for the parametric model. The reaction coefficient is chosen as σ = 10. The reaction is included mainly to add some local effects to an otherwise globally diffusive problem. A dictionary of 320 snapshots is created using randomly sampled diffusivity values. Fig. 3.3 shows some sample snapshots from the training set. Fig. 3.4 shows reconstruction errors using the dictionary-based approximation methods using (3.3). FOCUSS and OMP are used for sparse regression. As expected, the training set snapshots are almost accurately reconstructed with a single non-zero coefficient. In the validation set, reconstruction errors are observed in the order of 10 -2 to 10 -1 for sparsity target L = 10 and errors improve about an two orders for L = 30. The outcome of sparsity for each L, shown using a cross (×), is almost the same as the parameter L for the validation set. This is due to the algorithm trying to select as many vectors permitted by the user. Now that the dictionary methods have been demonstrated on unconstrained methods, we would like to extend these methods to inequality-constrained problems. However, it is not straightforward to apply random sketching to problems involving inequality constraints. The projection of inequality constraints on a randomized matrix like B creates linear combinations of constraints, which makes it harder to identify which constraints are being violated and which are not. In upcoming sections, two dictionary-based approximation schemes will be studied, but the projection of residual on B is applied in only one of them.

Dictionary methods for contact problems

Now that the feasibility of dictionary-based approximations have been established in an unconstrained problem, we move towards inequality constrained problems. To this end, two dictionary-based approximation methods for contact mechanics problems have been formulated, namely the greedy active-set approach and convex-hull approximation approach. To provide a broad overview of the features of these methods, Table 3.1 gives a list of features of these methods. The important differences between the two methods are bases types and the manner of enforcing constraint. These details will be discussed along with their formulations and application to contact problems are discussed in this section.

A greedy active-set method for dual dictionary element selection

To solve the contact problem using a dictionary of contact pressure snapshots, a greedy method inspired from the OMP and active set method is devised. From this point onwards, D d indicates a dictionary of dual snapshots of the contact problem. Consider the following sparse problem: (3.5b), (3.5c) and (3.5d) are based on the same principle as the reduced KKT conditions (1.14), but applied over a dictionary D d rather than a low-rank dual basis. Comparing (3.5) with the generic form of sparse regression statement (3.2), the balance equation is imposed as an equality constraint, instead of an inequality constraint where the residual norm is limited to certain tolerance. Imposing limits on the residual norm is usually more relevant when using random sketching, like in (3.3) when the system is under-determined. Also, it is not straightforward to use standard sparse regression algorithms in presence of the KKT inequality constraints (3.5c) and (3.5d). Equality constraint can be used as long as the system is determined, which is currently the case in (3.5b). In (3.5c), τ indicates a small penetration that has been permitted in the algorithm, which is necessary to ignore spurious penetrations created by the low-rank displacement basis Φ. Unless specified, the τ is set to a value in the same order of truncation tolerance δ (see Appendix A). The primary notion behind greedy algorithms is to enrich the approximation by a rank-1 term in each step. For standard sparse regression problems like (3.2), the OMP algorithm starts with an empty coefficient vector and adds a single non-zero entry in one step. The position of the new non-zero entry is chosen by searching the vector in the dictionary with the highest correlation to the residual. The greedy active-set algorithm devised here performs the greedy enrichment in an almost similar way but with a modification, i.e. the orthogonality with the violation of inequality constraint (3.5c) is used instead of the residual of system (3.5b). More precisely, the algorithm starts with a zero λ and in each enrichment step, the dictionary column with highest correlation to the violations of the non-penetration condition is added in the greedy enrichment. Thus, enrichment index of the dictionary p add can be computed as:

min λ 0 (3.5a) s.t. Φ T KΦ Φ T C T D d,I D T d,I CΦ 0 u λ I = Φ T f D T d,I g (3.
p add := arg max p D d [:, p] T CΦ u k -g +,τ (3.6) where (z) +,τ i := z i if z i > τ 0 otherwise
, which is a hard-thresholding operator.

(3.6) can be computed more efficiently using the reduced operators C(µ, u k-1 ) and g(µ, u k-1 ):

p add := arg max p C u k -g +,τ [p] (3.7)
where

C(µ, u k-1 ) = D T d C(µ, u k-1 )Φ g(µ, u k-1 ) = D T d g(µ, u k-1 )
However, to ensure that inequality constraints (3.5d) that prohibits non-negative pressures, it is necessary to not only enrich but also to eliminate terms that violate the constraints in each greedy step, which is simply done by eliminating the largest negative λ i.e. setting the largest negative entry to zero. The elimination index p rem is computed as:

p rem = arg min p λ k - [p] (3.8)
where (z) - i := z i if z i < 0 0 otherwise Though the motivation for this algorithm was based on the OMP to perform the dual dictionary element selection, the greedy active-set algorithm ends up being quite similar to the fixed-point active-set method. In this approach, only a maximum of one vector is activated in a given iteration thereby maintaining a small active set in intermediate iterations. The algorithm is given in Alg. 3.2, where Steps 8 to 10 contain the enrichment and the elimination process. In each iterative step, only one of the two operations is performed. Due to the architecture of the algorithm, the only steps that scale with the size of the dictionary are enrichment/elimination Steps 8 to 10. The size of mixed system in Step 6 is not significantly influenced by the size of active set during the intermediate iterations, thanks to the small sizes of active set I which is closely related to the sparsity of λ.

Application to the Hertz problem

This section deals with the implementation of the greedy active-set algorithm for the Hertz contact problem discussed in Section 2.4.1 of previous chapter (Fig. 2.5a). The parametrization also remains same, i.e. the imposed displacement on the top cylinder Build reduced constraint operators C(µ, u k-1 ) and g(µ, u k-1 )

6:

Solve the linear system (3.5b) Reduced model: Primal and dual snapshots on three training sets of sizes 12, 30, 60 and 120 are computed in the offline stage. The parametric points in each training set P tr are distributed uniformly in parametric space P =]0, 0.3]. The primal basis is computed using SVD and truncated using cumulatively normalized singular values using a tolerance δ. For the dual variable, a dictionary consisting contact pressure snapshots is used instead of a reduced basis. A spy plot of the dual dictionary for the training set of size 30 is shown in Fig. 3.5 for visualization purposes. A validation set P val with 119 points, which are the mid-points of the 4th training set is used to compute reconstruction errors.

Two instances of contact pressure reconstructions by the greedy active-set algorithm for parametric points outside the training set are shown in Fig. 3.6. The dictionary snapshots selected by the algorithm for the reconstruction are also shown. As expected, only a few dictionary elements, two in this case, are selected. The two selected snapshots are the ones whose contact position is similar to the required contact position for the specific instance. It indicates the greedy enrichment and elimination algorithm is selecting the correct snapshots for reconstruction. Mean reconstruction errors over the validation set vs. size of the training set are shown in Fig. 3.7a, along with the rank of primal basis. The reconstruction errors improve by more than an order of magnitude as dictionary size is increased from 12 to 120 for the case δ = 10 -10 . However, such gains with over-complete dictionaries diminish for the case δ = 10 -6 , especially the primal error improves only slightly. This is probably because at δ = 10 -6 , the primal truncation starts to become dominant source of error, rather than just the linear inseparability of contact pressure which seems to the case for δ = 10 -10 .

The computational time is computed as the total time taken for computation of Alg. 3.2, excluding the time taken in Step 5. The operator construction time is excluded because their efficient construction is not a part of this work. Mean computational time on the validation set (Fig. 3.7c) shows an initial increase with the dictionary size, but later settles down, but it is evident that the computational time is mostly influenced primarily by the number of iterations (Fig. 3.7b). This can also be seen in Fig. 3.7d, where the computational time per iteration is more or less flat w.r.t. the size of dictionary. The small increases in computational time per iteration with dictionary size can be explained due to the increase in the dimension of primal basis Φ, especially the case of δ = 10 -10 where the primal rank is strongly dependent on the training set size. Detailed plots of reconstruction errors over the validation set P val are shown in Fig. 3.8 for primal truncation tolerance δ ∈ {10 -10 , 10 -8 , 10 -6 }. Each curve corresponds to one of the training sets. Naturally, a larger training set is expected to have a lower (or at least the same) reconstruction error level. This anticipation is marginally satisfied in the case of δ = 10 -10 , but not in the case of δ = 10 -6 , where there are certain regions where the dual error for larger training sets is higher than a smaller training set. This happens due to the spurious penetrations discussed previously (and in Appendix A) and implies that even after relaxing the constraint by value of τ (see (3.5c)), some spurious selections are still possible for higher values of δ. Also, for all cases, the errors tend to be larger in general near d = 0.0, as the contact is sensitive to the parameter d around this region.

Sparsity pattern of the dual reduced dofs λ (i.e. the coefficients of the dual dictionary D d ) resulting from the greedy active-set approach for the training set of size 30 is shown in Fig. 3.9. As the snapshots in D d are arranged in increasing order of d, most of the columns of the sparsity pattern have two close non-zeros, indicating that nearest two snapshots to the targeted reconstruction were chosen by the algorithm. Spurious selections are seen in some rows, but the values of corresponding coefficients are quite small and do not have significant impact on solution.

Application to the Ironing problem

Unlike the Hertz problem, other contact problems might exhibit larger changes in the contact zone, like the ironing problem from Section 2.4.2 (Fig. 2.10). The effects of the strong variations of contact zone on the inseparability of the contact pressure were also discussed. Therefore, over-complete dictionaries might be more significant in generating an effective reduced model.

The parametrization of the horizontal displacement d x of the iron over the slab is considered for the reduced model. The vertical displacement d y = 0.3 is a fixed value. Given length of the slab L, the parametric space is defined by

P := {d x | 0 ≤ d x ≤ L}.
To study the influence of size of training sets on the reconstruction quality, training sets P tr are built in a nested manner which allows studying the evolution of reconstruction errors as a function of training set size i.e. size of the over-complete dictionary. Each nested level has uniformly distributed points over the parametric space P. The number of points is nearly doubled in the next nested level which also contain the points from the current level. Therefore, at the n-th nested level, 2 n + 1 snapshots are computed, and the contact pressure snapshots constitute the dual dictionary D d . The primal basis, however, is truncated by filtering the singular values by a fixed tolerance δ, and therefore its rank must be less than 2 n . Training sets of nested levels n = 3 to n = 7 are used to build the primal RB and the dual dictionary. The validation set P val contains all points from nested level n = 8 that are not present in level n = 7. Fig. 3.10 shows two instances of reconstruction of contact pressure along with the dictionary elements that were chosen by the algorithm. The algorithm selects snapshots that are closest to the required contact position. This can also be seen in Fig. 3.11, where the sparsity of dual dofs shows a clear selection pattern. Unlike the sparsity pattern of Hertz problem (Fig. 3.9), the ironing problem does not display spurious snapshot selections; possibly due to the stronger dependence of contact position of the ironing problem; unlike Hertz problem where contact pressure curves were centered around the same position. Therefore, the activation process is less ambiguous and less sensitive to primal truncations for the case of ironing problem. Fig. 3.12 shows the evolution of reconstruction (relative) errors, primal basis rank, number of iterations and computation time over the validation set P val vs. training set nested levels (log 2 (#P tr )). Means of the reconstruction related quantities over the validation set is shown for primal truncation tolerances of 10 -6 , 10 -8 and 10 -10 . The dual error improves, although slowly, with increasing size of dictionary, but the primal error shows signs of stagnation near level 6 and 7 and even earlier for the case of δ = 10 -6 . Errors for δ = 10 -8 and 10 -10 are almost the same, possibly because the inseparability of contact pressure and the finite size of dual dictionary are important sources of error. Conversely, in case δ = 10 -6 , the primal truncation seems to have more influence as primal error stagnates early, similar to the observations in case of the Hertz problem.

Moreover, the number of iterations and consequently the computation time seems to increase significantly with the dictionary size, in contrast to the Hertz problem where this relation was comparatively flat. A possible explanation for this behaviour could be that the larger dictionaries of the ironing problem, coupled with stronger parametric dependence of contact pairs, forces the greedy active-set to perform more iterations to find the appropriate dual dictionary candidates and contact pairs. Mean computation time per iteration shows trends that are a bit similar to the primal rank i.e. both primal rank and computational time per iteration are flat for δ = 10 -6 , whereas both are monotonically increasing for δ = 10 -10 . This was an expected observation since a large part of the system (3.5b) is related to the dimensionality of the primal basis Φ. The size of the dual dictionary, even though not truncated, does not seem to have any significant influence on the computation time per iteration.

A non-penetrating convex hull approach for monolithic dictionaries

As discussed in Remark 1.1, an interesting feature of mixed ROMs for incompressible flow problems is the computation of divergence-free RB. The advantage of this feature is that incompressibility constraint does not need to be explicitly enforced in the reduced problem. Moreover, orthogonality of a divergence-free velocity subspace and the irrotational pressure gradient term causes the pressure terms in the reduced problem to completely vanish (see [START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible navier-stokes equations: rigorous reduced-basis a posteriori error bounds[END_REF][START_REF] Liberge | Reduced order modelling method via proper orthogonal decomposition (pod) for flow around an oscillating cylinder[END_REF] for details). Computation of a divergence-free RB does not require any special treatment, as the divergence-free property of the training set snapshots is easily preserved in the resulting reduced subspace.

In this section, we attempt to extend a similar notion of constraint satisfying solution space to contact mechanics. But the approach is quite different compared to the incompressible flow problem because the nature of constraint is different. The incompressibility constraint, an equality constraint, applies globally over the domain; whereas the non-penetration constraint, an inequality constraint, applies locally on the contact surfaces and can be either active or inactive locally. Also, the displacement and contact pressure do not possess any special mathematical properties that cause the pressure term to vanish, as in the case of incompressible flow. In fact, the approach presented in this section, a reduced solution "space" is defined that naturally satisfies the non-penetration condition, but does not result in elimination of the contact pressure.

Unlike the greedy active-set algorithm previously explored which used POD basis for displacement and a dictionary for contact pressure, the following approach utilises a monolithic dictionary whose columns are snapshots consisting of both displacement and contact pressure. The reason for using a monolithic dictionary will be become evident in the formulation of this approach. Moreover, this method allows the residual to be projected on a low-dimensional space like it was done in the unconstrained problem in Section 3.2.2. It also permits the utilisation of sparse regression methods discussed in Section 3.2.1, provided the non-negativity constraint can be enforced. In this work, the non-negative version of the FOCUSS algorithm is used.

Consider a contact problem where the contact pairs are constant and independent of the parameter of the problem. The discrete inequality constraint can be expressed using matrices that do not change with the displacement field. The feasible region K in the displacement solution space can be defined as:

K := {u ∈ V | Cu -g ≤ 0} (3.9)
As the operators C and g are linear, it is easy to demonstrate the convexity of K. This would mean that a convex combination4 of non-penetrating snapshots would naturally satisfy the inequality constraint. The implication here is that if the reduced model explores the convex hull5 of the training set snapshots, the inequality constraint need ) is linear, convex and is made of segments corresponding to each element of the constraints Cu = g. In a high-dimensional problem, each segment of the constraint can be thought of as a hyperplane. The convex hull, which is simply a convex enclosure of snapshots, can also be thought of as a high-dimensional set.

not be explicitly enforced. This notion is visually shown in Fig. 3.13, where the convex hull of training set snapshots does not violate the convex non-penetration constraint. This property of convexity can be exploited to build efficient reduced models, if the following hypothesis holds:

Hypothesis 3.1 : Convex Subset hypothesis
Given a parametrized high-dimensional inequality-constrained problem with a convex feasible set K, its solutions lie in a low-dimensional convex subset inside K.

The Convex Subset hypothesis is an extension of the idea behind the low-rank hypothesis to inequality-constrained problems with convex feasible regions. This hypothesis allows us to explore a solution set limited to the convex hull defined by the snapshots without worrying about the non-penetration conditions. Note that the convexity of feasible region assures that any solution in the convex hull of displacement snapshots still lie inside the feasible region, independent of the Hypothesis 3.1. Therefore, in this approach based on the convex hull of the dictionary, the displacement field candidates built using convex combinations of the primal dictionary D p will be explored:

u ≈ D p u s.t. 1 T u = 1 u ≥ 0
As the non-penetration constraint will be satisfied naturally by all candidates in the convex hull of dictionary snapshots, a dictionary-based approximation can be computed using coefficients u that satisfy equilibrium equations. Therefore, the objective is to minimize the following residual:

r( u, λ) = K(µ)D p u + C T (µ, u)D d λ -f (µ)
As mentioned before, a monolithic dictionary will be used in this problem, in which each column consists of contact pressure snapshot stacked under the displacement snapshot. This means the coefficients u and λ will be same and will be denoted by α here onwards. Use of monolithic dictionary not only reduces the number of unknowns, but also makes the problem of minimizing the residual ||r(α, α)|| well-defined. Moreover, the residual can be projected on a low-dimensional subspace of the colsp(D p ), described by the matrix B containing its basis, while also adding sparsity constraints to the unknown α, like in (3.3). Thus, the problem can be stated as:

min ||α|| p (3.10a) s.t. B T r(α, α) < ε (3.10b) 1 T α = 1 (3.10c) α ≥ 0 (3.10d)
As the constraint (1.9c) is not imposed directly in this approach, the non-penetration condition is satisfied only because of the convexity of the feasible region. Thus, the equality condition C A (µ, u)u = g A (µ, u) for active constraints is also not imposed explicitly. Consequently, the complementary slackness KKT condition (1.9d) may not be satisfied exactly. This will be seen in the upcoming numerical examples.

The nnFOCUSS Algorithm: To solve the convex hull approximation problem (3.10), the non-negative version of the FOCUSS, the nnFOCUSS algorithm [START_REF]Reduced order and sparse representations for patient-specific modeling in computational surgery[END_REF][START_REF] Manucci | Sparse data-driven quadrature rules via p -quasi-norm minimization[END_REF] is used. The nnFOCUSS algorithm works by computing an appropriate relaxation parameter for each iterative update that maintains non-negativity. After computation of new coefficients at Step 5 of the Alg. 3.1, the relaxation step of nnFOCUSS algorithm is performed as follows:

if min(α k ) < 0 then ∆α = -(α k -α k-1 ) - α k ← α k-1 + min α k-1 ∆α ∆α

Element wise division end if

The algorithm must be initialized using a non-negative coefficients α 0 , which is computed using the non-negative least squares [START_REF] Lawson | Solving Least Squares Problems[END_REF] solution (which is not sparse, in general). The rest of the algorithm is same as standard FOCUSS. The problem (3.10) can be solved by plugging the following inputs in the nnFOCUSS algorithm:

D = B T K mono 1 T , x = B T f (µ) 1 , α 0 = nnls(D, x)
where K mono is the operator for the monolithic residual r(α, α):

K mono = K(µ)D p + C T (µ, u)D d
nnls indicates the non-negative least squares. An open source implementation of nnls provided in SciPy package [START_REF] Virtanen | SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF] is used. In this thesis, analytical validation for the convexity of feasible set and Hypothesis 3.1 will not be given. Instead, the following tests, that apply the leave-one-out approach on the training set snapshots, are proposed to analyze if these conditions hold within a reasonable error: Test A. Even if the operators C(µ, u) and g(µ, u) are non-linear and the feasible region K(µ) shows a dependence on the parameter µ, it may be reasonable to assume a nearly convex feasible set if the contact pairs do not change drastically. For such cases, the convexity of feasible set can be checked numerically using the training set snapshots and the corresponding non-linear operators. Given the constraint, C(µ, u)ug(µ, u) ≤ 0, the cross penetration (CrPen) error of the k-th dictionary element d across all constraints encountered in training set is computed as follows:

CrPen [k] = max [C(µ, d)d -g(µ, d)] +,0 ∀ d ∈ D p (3.11)
with d = D p [:, k] and D p = D p \d. If all values of CrPen are sufficiently small, it would be reasonable to assume the feasible set K(µ) to be convex. The notation [•] +,0 has the same meaning as in (3.6).

Test B. The Hypothesis 3.1 can be verified by approximating the left out snapshot with the convex hull of the rest of the snapshots using least-square criteria. The approximation error, referred to here as the convex hull least square (CHLS) error, of the leave-one-out approach can be used to assess the existence of hypothetical low-dimensional convex set. Each dictionary element d is approximated using the convex hull of D p and the error is recorded:

CHLS [k] = min α D p α -d 2 ||d|| 2 s.t.1 T α = 1 (3.12)
A high value of CHLS might occur due to lack of low-dimensional convexity or due to scattered training set data. The CHLS can be visualized using a lowdimensional problem shown in Fig. 3.14, where a point mass moves is permitted to move between a V-shaped obstacle. When leave-one-out approach is applied to a set of five snapshots, CHLS is high when the snapshot in the corner is left out. On the other hand, a low value of CHLS assures that it is possible to reconstruct reasonably using convex hull of the dictionary.

Illustrative example: elastic rope-obstacle problem

The convex hull approximation is first demonstrated in a problem where constraint operators are constant, so that convexity of feasible space is assured. The elastic ropeobstacle problem form Section 2.2 is considered, but with a different parametrization. In this case, the elasticity ν(x) of the rope is a function of parameter γ 1 .

ν(x)∇ 2 u(x) = f on x ∈ [0, 1] u(0) = u(1) = 0 u(x) ≥ -0.2(sin(πx) -sin(3πx)) -0.5
where

ν(x) = γ 1 if x < 0.5 30 otherwise on γ 1 ∈ [10, 50] (3.13)
Training set: γ 1 ∈ P tr = {10, 15, 20 . . . 45, 50}

Validation set: γ 1 ∈ P val = {12.5, 17.5, 22.5 . . . 42.5, 47.5}

The monolithic dictionary with 9 snapshots corresponding to training set defined above is computed in offline stage. The constraint is linear and independent of the parameter, and therefore, the feasible region is convex, which eliminates the necessity of Test A. Hence, only the Hypothesis 3.1 is tested using leave-one-out approach (Test B). The test results are shown in Fig. 3.15, indicating that hypothesis holds to a reasonable level.

The low-rank matrix B on which the residual will be projected is computed using left singular vectors of the truncated SVD B ← svd(D p , δ), with tolerance δ = 10 -7 . The reconstruction errors and resultant sparsity of the solution in the training set and a validation set is shown in Fig. 3.16. The training set is reconstructed within near numerical precision. In the validation set, the algorithm chooses the training set points that are nearest to γ 1 to create the best possible reconstruction with the given dictionary.

As discussed earlier, the KKT conditions, except the non-negativity of contact pressure, is not explicitly imposed. However, in this particular case with perfectly convex feasible region, one would expect that at least the non-penetration would be satisfied accurately. However, the linear combinations are not perfectly convex as the condition (3.10c) is satisfied only to a certain precision. The observed orders of these quanti- ties during reconstruction are given in Table 3.2. The validation set reconstructions in fact show penetration and non-zero complementary slack in the order similar as that of the deviation from convex combinations i.e. |1 T α -1|. Primal and dual sparsity are naturally equal as the dictionary is monolithic.

Demonstration on over-complete dictionaries of the Hertz problem

Hertz problem with same parametrization and training sets from Section 3.3.1.1 is solved using the convex hull approximation. For projecting the residual on a low-rank subspace, the matrix B is computed in two different ways:

B ← svd(D p , δ) or B ← orth(D p R)
where R is a matrix with random entries from a uniform distribution (3.3c). As discussed in Section 3.2.2, this matrix is used to create a smaller matrix by randomly combining dictionary columns.

As the constraint operators C and g show non-linear dependence w.r.t. u, it is necessary to test the convexity of the feasible set K(µ) (Test A). All CrPen errors ( CrPen ) were found to be quite small, in the order of O(10 -7 ) for all combinations of d and d in (3.11). Test B is performed to check the Hypothesis 3.1, with results shown in Fig. 3.17. The reconstruction errors for the reduced model can be expected to be more than the CHLS shown in the figure. For large dictionaries, it makes sense to compare the performance of the method using various δ and sizes of R. The mean reconstruction errors over the validation set vs size of dictionaries can be seen in Fig. 3.18. The use of randomized matrices R naturally generates more error than truncated SVD, as the truncated SVD selects the best representative subspace for the given training set. Similar to the greedy active-set method of Alg. 3.2, the improvement in error flattens out for large dictionaries. On the other hand, the number of iterations and computation time seems to be nearly independent of the dictionary size.

The non-convex case of the Ironing problem

The convexity test for the feasible set (Test A) is performed on the Ironing problem. Since, the contact zone changes significantly, the test performed using the training set snapshots fails as expected. The CrPen error was found in O(10 -1 ). Therefore, the ironing problem cannot be solved using the current formulation of convex hull approach. The non-convexity can also be demonstrated by a convex combination that violates the non-penetration condition significantly, as shown in Fig. 3.20. For the two snapshots shown in this figure, the slab surface is deformed at different positions. The convex combination where the two displacement field snapshots were equally weighted, generated a deformed configuration with the slab being deformed at the two locations corresponding to each of the snapshots. But no deformation was generated at the position of the iron, causing penetration at that position. In this case, local behaviour of the displacement field leads to the non-convexity of the ironing problem.

Discussions

The application of dictionary methods was proposed to mitigate the challenges posed by linear inseparability of the contact pressure. As the lack of local information prevents the creation of a robust RB, over-complete dictionaries have the potential to preserving local information computed in the offline stage, thereby providing a better chance at resolving the position of contact more effectively. The same is demonstrated by building reasonably accurate reconstructions using dictionary methods.

For the reconstructions of Hertz problem using the greedy active-set algorithm, the reconstruction errors decay consistently when richer training sets are used, and the computational time per iteration seems to be fairly constant, but change slightly with the primal rank.

The case of ironing problem where the contact position changes significantly is also considered. However, the decay of reconstruction errors vs. size of dictionary curves were found to be flatter compared to the Hertz problem. This is probably due to stronger dependence of contact position on the parameter of the problem. Though the dictionary size didn't have any significant influence on the computational time per iteration, the number of iterations increased significantly for larger dictionaries in case of ironing problem.

Comparing the two dictionary methods: the greedy active-set algorithm and the convex hull approximations, the reconstruction errors of the Hertz problem were better using the greedy active-set algorithm (see Figs. 3.7a and 3.18). This observation can be possibly attributed to the use of the monolithic dictionary in case of the convex hull approximations; which forced common interpolation coefficients for the displacement as well as contact pressure. In particular, the convex hull approach shows large errors for small training sets, as the convex hull of a small set of snapshots might be a very small subregion in the entire solution space.

Another possible source of error in convex hull approximations is that the snapshots in the dictionary may not be the best ones to explore the low-dimensional convex subset. On the other hand, the computational time (and number of iterations) for the convex hull approach seems to be more independent of dictionary size, compared to the greedy active-set algorithm (see Figs. 3.7b, 3.7c and 3.19), thanks to the elimination of active set selection process.

In case of the convex hull approach, the reconstruction error seems to be more sen-sitive to the size of dictionary than the truncation tolerance/rank of the orthogonal basis B. Also, for the greedy active-set method, the sensitivity to truncation seems to be low between 10 -10 and 10 -8 . Moreover, in both approaches, the computational time per iteration is nearly flat when the rank is low. Thus, one can infer that that large dictionaries and moderate truncation tolerances / rank of B are an optimal choice. One could argue that the number of iterations increases significantly with size of dictionary in case of the ironing problem solved using the greedy active-set method (Figs. 3.12b and 3.12c), however, this is related to the inherent non-linearity of the contact problem. The number of iterations could possibly be reduced by using smarter initialization of the contact pairs according to the queried parametric value rather than the initializing with reference state. For the convex hull approach, the primary disadvantage is the assumption of convexity of the feasible region, limiting the applicability to convex problems. The non-convex nature of ironing problem meant that convex hull approximation method could not be used in this case.

Perspectives

Efficient construction of constraints

In the Alg. 3.2, Step 5 can be the most expensive step if full operators are built from scratch and then projected on the reduced bases. The constraint operators C and g show non-linear dependence on the displacement field u and therefore must be reconstructed in every iteration. An efficient construction of the constraint operators for low-rank methods is discussed in [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF] using the Empirical Interpolation Method (EIM), where affine decompositions of the distance functions are computed offline. This enables efficient construction of constraint operators in online phase by sampling distance functions at a few points on the contact surface, followed by computing the coefficients of the affine form.

The distance function k(µ, v, v) of (1.2) when expressed in discrete form (as discussed in [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]) turns out to be a sparse vector. The spy pattern of this vector can be a strong function of the parameter when contact area changes significantly. This is caused by the dependence of the distance function on the contact pairs established using the current configuration. Hence, it is bound to exhibit some form of linear inseparability just like contact pressure. Therefore, it is possible that linear inseparability issues might also appear in generating an effective affine decomposition.

Application of sketching methods to dictionary-based approximation

As discussed briefly in Section 3.2.2 and in detail in [START_REF] Balabanov | Randomized linear algebra for model reduction-part II: minimal residual methods and dictionary-based approximation[END_REF], random sketching can be used to project the residual of a dictionary-based approximation problem on a randomly sampled subspace of the dictionary. By controlling the size of the subspace, one can control the cost of solving the dictionary-based approximation problem, relatively independent of the dictionary size. The difficulties in applying these techniques to contact pressure dictionary appear due to the inequality constraints. If random sketching methods are applied, the inequality constraints are linearly combined in random fashion leading to loss of information on violation (penetration) and satisfaction (non-penetration) of constraints.

Smart snapshot selection for Convex hull explorations

RB methods equipped with error estimators sample the parametric space efficiently. Similarly, it might be possible to extend the same notion to efficiently sampling the convex feasible set. The convex hull approximations might be more robust if the snapshots lie on the edges/vertices of the low-dimensional convex subset presumed in Hypothesis 3.1. For e.g. convex hull approximation for the low-dimensional problem in Fig. 3.14, it is important to include the snapshot near the corner of the feasible region in the training set for a robust reconstruction. However, deducing the corners in a high-dimensional problem may not be geometrically intuitive as the low-dimensional problem in the figure. pressure manifold

In Section 3.1, the philosophy of dictionary-based approximation methods for contact problems was discussed. The principal argument in favour of over-complete dictionaries was the lack of low-rank structure due to linear inseparability of the contact pressure field, which has two consequences: (a) a large number of snapshots might be necessary to truly discover the underlying structure of the contact pressure field, and (b) linear model order reduction methods cannot capture local effects of contact pressure in a low-rank subspace. The focus of Chapter 3 was on efficient handling large dictionaries of contact pressure snapshots in the online phase, but does not address the high offline cost of creating such dictionaries. This chapter focuses on the use of non-linear dimensionality reduction methods as means to improve reconstruction errors without paying the large computational cost of creating over-complete dictionaries.

Linear Subspaces vs. Low-dimensional Manifolds

Many of the dimensionality reduction methods work on the principal assumption that the underlying physics occupies a subset of the solution space defined in the full model. At this point, the question arises about the nature of the aforementioned subset of solution space. Low-rank methods assume that this subset is a low-dimensional linear subspace 1 of the full solution space. The assumption that such a low-dimensional subspace exists implies that the solution of any point in the particular parametric space may be approximated using a linear combination of a limited number of precomputed basis vectors that naturally lie in this subspace. As seen in Chapter 1, methods like Proper Orthogonal Decomposition (POD) are capable of extracting a basis for the linear subspace that describes the underlying physics. However, in some problems the parametric solutions occupy a low-dimensional manifold, but not necessarily linear subspace. Consequently, it is not possible to compute a robust low-rank basis, as shown for the contact pressure field in Chapter 2. Nonetheless, there exists a paradigm called non-linear dimensionality reduction that aims at computing the underlying structure by finding the associated non-linear behaviour. A subclass of non-linear dimensionality reduction known as manifold learning is quite extensively used in the analysis and visualisation of high dimensional datasets [START_REF] Bregler | Nonlinear image interpolation using manifold learning[END_REF][START_REF] Ghodsi | Dimensionality reduction a short tutorial[END_REF]. The fundamental idea behind manifold learning is similar to low-rank methods and is known as the manifold hypothesis [START_REF] Fefferman | Testing the manifold hypothesis[END_REF]. Roughly speaking, a manifold can be described as an n-dimensional space that is globally non-linear in general but locally resembles a linear vector space.

The manifold hypothesis assumes that high-dimensional data sourced from practical applications often lie in low-dimensional manifolds. The low-rank hypothesis is a special case of the manifold hypothesis involving the assumption of the low-dimensional manifold being a linear subspace. Manifold learning methods rely on finding the non-linear lowdimensional manifold within the high-dimensional data. A well-known application of manifold learning is the case of a data set containing gray-scale images of an object rotated by various angles. Linear interpolation methods like PCA fail to reduce this type of data, but manifold learning methods can extract the inherent one-dimensional space [START_REF] Bregler | Nonlinear image interpolation using manifold learning[END_REF][START_REF] Ghodsi | Dimensionality reduction a short tutorial[END_REF].

To better approximate the manifold structure of a parametrized problem, partitioning of the parametric domain [START_REF] Dihlmann | Model reduction of parametrized evolution problems using the reduced basis method with adaptive time-partitioning[END_REF] and a local reduced basis methods [START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF] have been explored. Interpolation in the tangent space to the Grassmann manifold structure have been used to construct reduced models with a wide parametric range [START_REF] Amsallem | Interpolation method for adapting reduced-order models and application to aeroelasticity[END_REF].

One of the most commonly applied non-linear dimensionality reduction techniques in the parametrized ROMs is the kernel methods [START_REF] Haasdonk | Greedy Kernel Approximation for Sparse Surrogate Modeling[END_REF][START_REF] González | kPCA-Based Parametric Solutions Within the PGD Framework[END_REF][START_REF] Díez | Nonlinear dimensionality reduction for parametric problems: A kernel proper orthogonal decomposition[END_REF]. Kernel interpolation methods along with greedy algorithms to create nested training sets to create efficient surrogate models are discussed in [START_REF] Haasdonk | Greedy Kernel Approximation for Sparse Surrogate Modeling[END_REF]. The kernel-PCA (kPCA) approach can also be used for the computation of intrinsic parameters that are not explicitly included in the mathematical model [START_REF] González | kPCA-Based Parametric Solutions Within the PGD Framework[END_REF]. An application of kPOD to an advection-diffusion problem with inseparable nature can be found in [START_REF] Díez | Nonlinear dimensionality reduction for parametric problems: A kernel proper orthogonal decomposition[END_REF], where the effectiveness of kernel-based method is shown using the higher concentration of variance in the first singular values of kPOD compared to the case of linear POD.

A known drawback of kernel-based methods and some other manifold learning techniques is that backward mapping, i.e. reconstruction of a high-dimensional entity from low-dimensional data is not as trivial as it is done in linear dimensionality reduction methods such as PCA. Usual the reconstruction is an approximation to the local tangent space of the manifold, computed using a linear combination of snapshots with weights computed by employing a minimization of a discrepancy functional criteria [START_REF] García-González | A kernel Principal Component Analysis (kPCA) digest with a new backward mapping (pre-image reconstruction) strategy[END_REF]. An improved approximation based on a quadratic tangent space is discussed in [START_REF] Díez | Nonlinear dimensionality reduction for parametric problems: A kernel proper orthogonal decomposition[END_REF].

Hyperbolic systems are often associated with moving discontinuities, which also leads to lack of a linear subspace. An approach involving the transformation of the physical variable has been used in such problems to separate the discontinuous and continuous behaviour [START_REF] Bansal | Model Order Reduction Framework for Problems with Moving Discontinuities[END_REF]. A related approach aimed at generating transformed snapshots thereby aligning the discontinuities has been demonstrated on Burger's equation [START_REF] Welper | Transformed snapshot interpolation[END_REF]. As it will be seen, the ideas given in this chapter have similar motivations of aligning the contact pressure snapshots, but without transforming the physical variable.

Contact pressure as a non-linear manifold

In Chapter 2, qualitative and quantitative arguments about the linear inseparability of contact pressure have been made. The lack of low-dimensional linear subspace, i.e. a low-rank subspace for contact pressure was demonstrated. In this section, we discuss the idea of contact pressure lying in a low-dimensional non-linear manifold. Looking back at Fig. 2.5b and Fig. 2.11, the snapshots in respective figures have similarities that are evident intuitively. Given these similarities, there is a strong possibility that these snapshots exist in a low-dimensional non-linear manifold.

This intuitive hypothesis can be supported by Fig. 4.1. Contact pressure snapshots of the ironing problem are shown by the lines as a function of the parameter d (the horizontal displacement of the iron), and the green translucent surface is generated using the surface plot of the same data. One can intuitively see that the contact pressure is a non-linear function of just one parameter, which happens to be also the physical parameter d in this case. The fallacy of low-rank methods is that they attempt to perform a linear interpolation in direction d, when in fact the relationship is non-linear. The idea that will be exploited in this chapter is that an appropriate non-linear transformation applied on this surface in a way that allows interpolation along the "aligned" direction shown in the figure. This will possibly improve the linear separability and reducibility.

Dimension warping

The core idea in manifold learning methods (such as kPCA) is that linearly inseparable data can be appropriately mapped into a higher dimension where it is relatively more separable. The mapping to high-dimensional space has to be essentially non-linear, since a linear mapping will not resolve the inseparability issues.

Consider a non-linear mapping from the spatial domain of contact surface Γ to a warped domain Γ. The warped domain Γ need not admit a physical interpretation.

ϕ(µ, •) : Γ → Γ ζ ← ϕ(µ, x) (4.1)
where ζ is a coordinate in Γ. The notation x(µ, ζ) is used to indicate x that maps to ζ.

A qualitative argument can be made regarding an effective linear behaviour in Γ space. The mapping ϕ must warp the spatial dimension so that the regions with nonzero contact pressure must be mapped to a common warped region, and the gap (noncontact) regions must be mapped to other regions in warped space. In other words, let the warped space be a union of disjoint sets, Γ = Γ c ∪ Γ /c , then for every pair (x, ζ) the mapping ϕ must satisfy the following:

ϕ(µ, •) : Γ → Γ ζ ∈ Γ c if λ(µ, x) > 0 ζ ∈ Γ /c if λ(µ, x) = 0 (4.2)
where Γ c and Γ /c are the regions where contact and no-contact zones are mapped respectively. Any function f defined on Γ has an equivalent f defined on Γ domain and can be related to each other with the following relation:

f (ζ) = f (x)
Remark 4.1 : Preserving the split of contact and gap zones

Extending the above qualitative argument to the case where the contact pressure snapshots have multiple disjoint regions of contact and gap zones, each subregion should be mapped to corresponding subregions in warped space. For e.g., in the ironing problem, each snapshot can be split into three zones: one contact zone and two gap zones. Thus, there should be three disjoint zones in warped space corresponding to each zone of the original space.

Dynamic Time Warping

It is not straightforward to compute the mapping ϕ(µ, •) and the warped space Γ with the desired properties discussed in the previous section. Algorithms that identify regions of similarities between two signals can prove to be useful to generate mappings that improve separability. One such method, known as Dynamic Time Warping (DTW), is an algorithm built to measure similarity between two discrete signals containing similar but displaced features. The algorithm is named so because it was built for measuring similarities in temporal sequences involving speech or motion [START_REF] Rabiner | Pattern-Comparison Techniques[END_REF][START_REF] Müller | Dynamic Time Warping[END_REF]. Other methods exist that deals with sequence similarities, such as elastic matching [START_REF] Talmi | Template matching with deformable diversity similarity[END_REF]; but they are not studied in this thesis. In this section, we illustrate the application of DTW to generate mappings that satisfy the conditions specified in (4.2). 2 Though only unidimensional DTW has been explored here, it is extendable to higher dimensions [START_REF] Shokoohi-Yekta | Generalizing DTW to the multi-dimensional case requires an adaptive approach[END_REF].

The DTW algorithm works on two given signals and computes similarities based on distances between pairs of points. It matches each point on the first signal to at least one point on the second, and vice-versa. A detailed explanation of DTW can be found in [START_REF] Rabiner | Pattern-Comparison Techniques[END_REF] and Python implementations can be found in the dtw module [START_REF] Giorgino | Computing and visualizing dynamic time warping alignments in r: The dtw package[END_REF]. The mapping ϕ(µ 1 , x) and ϕ(µ 2 , x) generated by DTW is described by indices of the discretization of the coordinate x. This is illustrated using a simplified example.

An illustration of DTW: Consider a set of n discrete points x in [0, 1] and two discrete vectors λ k ∈ R n , k ∈ {1, 2}. For a simple illustration, the vectors are chosen such that each has a peak (with value 1.0) at different positions surrounded by vanishing areas (values are given in Table 4.1a). The DTW alignments (computed using [START_REF] Giorgino | Computing and visualizing dynamic time warping alignments in r: The dtw package[END_REF]) are given by indices I k ∈ R n . The index at I 1 [m] is matched to the index at I 2 [m] of respective λ vectors (in Python notation). Usually n ≥ n, where the equality holds if alignment happens to be a one-to-one mapping.

The alignments are given in Table 4.1b and are shown graphically in Fig. 4.2. One can see that the alignments connect not only the peaks of the two vectors, but also the vanishing areas on the either sides as well. The alignments I k are equivalent to ϕ(µ k , •) of (4. Though a detailed explanation of the DTW algorithm is beyond the scope of this thesis, the major features of uni-dimensional DTW are listed here:

1. The first (and the last) points of the two input vectors are always connected, even 4.4b. Notice that the alignments were able to connect respective contact and no-contact zones of both curves. Also, in the warped spatial dimension, respective contact and no-contact zones of both the warped snapshots are the same. In other words, condition (4.2) is satisfied. Also, looking back at the "aligned" direction in Fig. 4.1, we can say that the warping has rendered this direction orthogonal to the warped coordinate ζ, thereby making linear combinations useful. The mapping generated by DTW is linear w.r.t. individual snapshots. However, the mapping is not only different for each snapshot, it depends on the pair of snapshots taken into account. These properties mean that the mapping is essentially nonlinear.

Interpolating in warped space

Although the DTW algorithm is built to compute similarities between signals, we intend to use it as an interpolation tool. Interpolations using DTW have been previously performed to reconstruct missing frames in fast motion detection [START_REF] Almog | Spatial de-interlacing using dynamic time warping[END_REF]. The objective of performing DTW on pairs of contact pressure snapshots is to use the DTW alignments and warped space for exploring the contact pressure manifold without generating intermediate snapshots directly with FE solver. Fig. 4.5 shows a demonstration of an interpolation in the original and the warped spatial dimension. Though interpolating along the alignments in the original dimension looks intuitive (Fig. 4.5a), only the interpolation in warped space is a straightforward linear operation (Fig. 4.5b). Coefficients α and 1α are used as weights for each snapshot in this operation, and the case for α = 0.5 is shown in the figure. In principle, the value of the parameter µ that corresponds to the interpolated snapshot is not known.

As interpolation is performed in the warped dimension, it means that a backward mapping is also needed to transform the interpolated curve back into the original spatial dimension. However, the mapping ϕ(µ, •) is a function of the parameter and the mapping for the interpolated curve needs to be computed. Formally, the interpolations on the two snapshots, and the mappings for parametric values µ 1 and µ 2 can be written as:

λ int (α, ζ) := α λ(µ 1 , ζ) + (1 -α) λ(µ 2 , ζ) (4.3a) x int (α, ζ) := α x(µ 1 , ζ) + (1 -α) x(µ 2 , ζ) (4.3b) ϕ int (α, x) := αϕ(µ 1 , x) + (1 -α)ϕ(µ 2 , x) (4.3c)
where λ int and ϕ int are interpolated contact pressure and mapping respectively and α is the linear interpolation parameter between the two snapshots. The interpolated contact pressure in original and warped space can be related as:

λ int (α, x) = λ int (α, ϕ -1 int (α, ζ)) = λ int (α, ζ) (4.4)
Once DTW alignments are available, (4.3) and (4.4) can be easily implemented in a short routine dtw interp to perform the interpolation.

Notation 4.1 :

There is a slight abuse of notation in (4.3), (4.4) and in interpolated quantities here onwards. The first argument of quantities • int (•, •) is the interpolation parameter α and not the physical parameter µ.

Notation 4.2 :

Interpolated quantities • int (•, •) also depend implicitly on µ 1 and µ 2 , which are not shown in the arguments to simplify the notation.

dtw interp is given in Alg. 4.1. It takes two discrete contact pressure snapshots λ k ∈ R n , k ∈ {1, 2} and the corresponding spatial coordinate x ∈ R n i.e. there are n discrete dofs on the contact surface. The DTW package computes alignments and a warped space of dimension n. Indices I k provide the mapping to the warped space for respective λ k . The warped vectors λ k ∈ R n can be linearly combined to create the interpolated vector λ int ∈ R n and the corresponding spatial coordinate x int ∈ R n . The interp1d function in Step 5 is borrowed from the Python module scipy.interpolate [START_REF] Virtanen | SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF]. This function is used for changing the discretization from x int to x using linear piecewise interpolations. Steps 3 and 4 corresponds to (4.3a) and (4.3b) and Step 5 corresponds to (4.4).

Algorithm 4.1 dtw interp

1: Input: x, λ 1 , λ 2 , α 2: I 1 , I 2 = dtw(λ 1 , λ 2 ) I 1 and I 2 map indices to warped dimension 3: λ int = αλ 1 [I 1 ] + (1 -α)λ 2 [I 2 ] λ j [I j ] is the warped vector λ j 4: x int = αx[I 1 ] + (1 -α)x[I 2 ]
interpolating in warped spatial coordinate 5: λ int ← interp1d( λ int , x int )(x) linear interpolation from x to x 6: Output: λ int The application of dtw interp is demonstrated on the illustrative example of Table 4.1 in Section 4.2.1. The interpolated quantities for this example are shown in Fig. 4.6.

Remark 4.3 : Interpolation using the alignments I k

In Alg. 4.1, the computation of λ int from λ int can be implemented in an alternate way. The interpolation of the mapping (4.3c) can be computed by interpolating between I 1 and I 2 to generate a set of interpolated indices. The steps would read as:

I int = (αI 1 + (1 -α)I 2 ) + 1 2 λ int [I int ] ← λ int
same as applying inverse mapping ϕ -1 int involving a round off operation to the nearest integer, expressed here using the greatest integer function • and the offset of 1 2 . This is equivalent to using the inverse mapping ϕ -1 in (4.4). However, this is approach is not followed in Alg. 4.1, as the round off operation is sensitive to truncation at machine precision level, especially for α = 0.5, leading to improper interpolation. Instead, x int that corresponds to the warped coordinate ζ is used for interpolating back to x. Also, the resulting I int contains non-unique entries, leading to non-injective inverse map ϕ -1 int . The quality of interpolation using DTW on the ironing problem is quite good as the snapshots have a similar contact pressure profile, with only the contact zone being translated in the x-direction. DTW-based interpolations of the elastic rope-obstacle problem and the Hertz problem are also shown in Figs. 4.7 and 4.8. In these cases contact zone does not translate but magnifies. DTW is able to interpolate the contact zone correctly, but the interpolated pressure profile has some aberrations near the peak pressure of the Hertz case. This is explained by the concentration of the DTW alignments near the peak i.e., many points from λ(µ 2 ) are mapped to the peak of λ(µ 1 ).

DTW-based enrichments for contact mechanics

As mentioned at the beginning of the chapter, the idea is to use manifold exploration to reduce the offline cost. Since the main cause of linear inseparability is the geometrically local nature of contact pressure field, the ability to scale and translate local effects across the domain of interest seems promising. Also, the displacement field does not have linear separability issues and a robust linear subspace seems to be readily available (Fig. 2.14) even if parametric space is relatively unexplored in offline stage. In an optimistic scenario, if the warped interpolations can reduce the size of training set to the same level necessary for capturing the displacement field, the bottleneck posed by linear inseparability will be effectively removed.

We propose the use of DTW as a warped interpolation tool that can generate intermediate contact pressures that haven't been explored in the offline stage. This is an inexpensive way of enriching the contact pressure dictionary; since the computation cost of warped interpolations is negligible compared to generating proper snapshots. Fig. 4.9 shows a demonstration of such an interpolation for three different values of α. It is no-ticeable that the peak position of contact pressure of the interpolated curves relative to the peaks of the original snapshots seem to have a nearly linear relation with α, though this is a particular feature of the ironing problem and may not apply in general. For an effective enrichment, the interpolation should not be done on every pair of contact pressure curves. Rather, each pair of curves should be reasonably similar in some sense. Similarity can be computed using a measure of distance which can be either the L 2 norm of the difference or the distance measure defined by the DTW algorithm. Snapshots that can be considered similar according to the distance measure will be referred to as neighbours here onwards. For each pair of neighbours, the DTW alignments are computed and interpolations are performed with a set of M parameter values {α m } M m=1 . This process is shown in Alg. 4.2, where neighbours of currently active vectors are computed using the distances between vectors of the dual dictionary in a matrix Dist. The entry Dist[i, j] has the distance between D d [:, i] and D d [:, j], and hence, Dist is a symmetric matrix with main diagonal containing only zeros.

The effectiveness of DTW-based enrichments can be shown using the compactness metric. Recall that the "dual cone" compactness metric defined in Section 2.3 can be used for as a measure of linear separability. The comparison of this metric for the original and enriched CPG dual basis is shown in Fig. 4.10. The compactness of the enriched basis was computed by performing 19 equispaced DTW interpolations (i.e. with α m = 0.05m and M = 19 as input for Alg. 4.2) between each pair of neighbours in the original dual basis. The compactness of the enriched basis is considerably lower than the original, indicating that DTW enrichment has introduced the information that was missing in the cone of the original basis. if pair (i, j) is not repeated then (i, j) is an unordered pair

Remark 4.4 : Indefiniteness of the DTW kernel

A more appropriate way of testing the separability of given data in warped space is by observing the decay of singular values (equivalent to the dual orth compactness in Chapter 2). However, computation of singular values requires that the DTW mapping generates a symmetric positive definite (SPD) Gramian matrix G ij = λ i , λ j , where each entry of G contains an inner product of a pair of warped snapshots. Such a comparison would have been possible if SPD kernels are used for computation of the Gramian matrix, as done in [START_REF] Díez | Nonlinear dimensionality reduction for parametric problems: A kernel proper orthogonal decomposition[END_REF] for an advection-diffusion problem. However, instances of Gramian matrix computed using DTW kernel being indefinite are known [START_REF] Lei | A study on the dynamic time warping in kernel machines[END_REF]. The Gramian matrix computed for this case also resulted in negative eigenvalues, hence the dual orth compactness cannot be compared.

The richness of D adapt d w.r.t D d can also be compared using the nested errors defined for the ironing problem (see Section 2.4.2, Fig. 2.14). Let us recall that the nested error (2.7) was defined as the projection error of (n + 1)-th level basis on n-th level basis. The comparison of nested errors with and without DTW interpolation is shown in Fig. 4.11. For each level n, the dictionary neighbours are interpolated using Alg. 4.2 with I = {1, 2, 3, . . . #P tr }, N = 2, M = 1 and α 1 = 0.5. The ironing problem has a peculiar property where the position of contact pressure peak is more or less same as the parameter d. Therefore, the chosen N, M and α 1 form a near best case scenario of the approximation because it suits the distribution of nested training set points i.e. (n+1)-th level contain mid-points of the points in n-th level. Thus, a considerable decrease in the nested errors is seen for smaller training sets. Naturally, the gains in dual orth error due to DTW decreases as the training set becomes larger, as the snapshots in large training sets are closer to the target solutions. Projection error of nested level n + 1 snapshots on dictionary of level n. Note that solid dual orth and dual cone lines are similar to the ones in Fig. 2.14. The dual orth points not seen in the figure drop to the order of 10 -14 (nearly machine precision) However, in a general scenario, the value of an optimum α may not be trivial. Therefore, to arrive at an appropriate value of α, it is necessary to explore multiple DTW-based interpolations corresponding to different α's. This is the reason a set of M values {α m } is supplied to Alg. 4.2. The enrichment process can be strategised in two different ways:

• Number of interpolation points (M ): For a given set of snapshots generated in offline phase, warped interpolation is performed between every pair of neighbouring snapshots with a certain resolution of in ]0, 1[. For a uniformly distributed points, α m = m/(M + 1). The interpolations in Fig. 4.9 is such a case with M = 3.

• Number of enrichment levels (L): It is possible to continue the enrichment process further by creating multiple levels of interpolated vectors. Suppose, the original dictionary contains two vectors λ 1 and λ 2 and interpolation is performed with M = 1 and α 1 = 0.5. Thus, an interpolated vector, let us call it λ 1.5 , can be computed using DTW interpolation in the first level. In a second level, further interpolation can be performed on the pairs {λ 1 , λ 1.5 } and {λ 1.5 , λ 2 } to generate λ 1.25 and λ 1.75 .

Further, these enrichments can be done during either the online or offline stage, depending on the desired distribution of the computational cost. If done in the offline stage, warped interpolations between every possible pair of neighbours need to be computed. On the other hand, in the online stage, interpolations can computed on a smaller set of vectors, i.e. the active dictionary vectors and their neighbours. In this thesis, the enrichment in the online stage is studied, as the focus is on the enhancement of the reconstruction accuracy using DTW and the study of computational cost is reserved for future study. However, a preliminary comment can be made in this aspect, that the DTW-based enrichment is not a computationally intensive task. As the Alg. 3.2 generates a small active set, the number of interpolations performed is not large. The resulting algorithm is given in Alg. 4.3, with a parameter L denoting the number of enrichment levels. In the numerical results that follow, the parameters used are L = 5 and M = 4 with the interpolation points {0.2, 0.4, 0.6, 0.8} M m=1 . Moreover, the improvement in reconstruction is seen more in the dual field rather than the primal. This is an expected observation as any improvement in primal reconstruction due to dual dictionary enrichments is caused indirectly through the mixed system (3.5b).

On the other hand, there is no improvement in the reconstruction of Hertz problem (Fig. 4.12a). One likely reason could have been that the contact pressure snapshots are relatively close, unlike the ironing problem, and therefore the DTW-based interpolations are unable to offer any significant new information. However, this contradicts the observations in Fig. 4.10a where improvements in dual cone compactness is observed. A more likely reason could be, that to improve in reconstruction of Hertz problem, the primal dictionary also needs enrichment rather than just the dual dictionary. However, artificial enrichment could not be applied to the primal dictionary (as it needs a multidimensional DTW implementation), and the dual dictionary enrichments alone are ineffective at improving the errors. This reasoning is more consistent with observation that reconstruction error improves only with increase in training set size and not with DTW-based enrichments. Another possible source of error is the distorted DTW-based interpolation of the contact pressure snapshots, as shown in Fig. 4.8. Though, in reality, the distortion is not as pronounced as shown in the figure, since the interpolation done in Alg. 4.3 is carried out between snapshots that are closer.

Conclusions

The DTW enrichments proved to be useful in improving contact pressure reconstructions for a problem with higher variations in contact area (ironing problem), especially when the dictionary is not rich enough to resolve the local effects. The improvements obtained using DTW diminish with increasing training set size, as the local information related to contact pressure start appearing in the dictionary directly. For problems with less variations in contact pressure (Hertz problem), DTW does not provide any improvements in reconstruction errors.

Overall, it is fair to say that non-linear interpolations can prove useful to train models for problems with large variability in contact positions, especially if the complexity of high-fidelity problem makes it difficult to compute a large training set. This is a strong motivator to test other non-linear reduction techniques such as kPOD on such problems. However, whether it is trivial to enforce inequality constraints in non-linear frameworks or not is yet to be explored.

Persepctives

The 3D case

Extension of DTW based interpolation to 3D problems (with 2D contact surfaces) may not be straightforward. Though the DTW algorithm was originally built for a 1D signals, but a 2D version has been developed [START_REF] Shokoohi-Yekta | Generalizing DTW to the multi-dimensional case requires an adaptive approach[END_REF]. However, the 2D DTW is limited to structured meshes, thereby limiting its application. Even if the geometry of the contact . Reconstructions performed using Alg. 3.2 and Alg. 4.3 for respective dictionaries. DTW enrichments seem to be more impactful in the ironing problem than Hertz problem. Primal and dual errors are computed using H 1 and L 2 norm, respectively. surfaces allow discretization using a structured mesh, it may not be the optimal mesh for generating the snapshots using finite elements. It may however, be possible to compute the FE snapshots using a generic mesh and then interpolating the solution on a structured mesh during a post-processing, thereby making it suitable for DTW-based interpolations.

Efficient evaluation of non-linear terms

To solve the contact problem using the dictionary method, it is necessary to compute the operators C and g, which are dependent on the dual dictionary D d . However, if DTW based interpolation is used, the dictionary D d is dynamic. This can create additional difficulties in devising affine decompositions of these operators. If EIM approach from [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF] is used, computation of the independent terms C p and g q in (1.17) becomes dependent on the dual basis/dictionary. An efficient way of handling such situations is necessary. A straightforward but less efficient way might be to build an affine decomposition of high-fidelity operators C and g, and update the reduced operators whenever the operators or the dictionary is updated.

Non-linear interpolations of primal dictionary

As seen in the case of Hertz problem, the dual dictionary enrichment is not effective at improving reconstruction errors. Also, in case of the ironing problem, mainly the dual reconstruction errors show improvement than the primal one. This observation opens up the possibility of enriching the primal dictionary using non-linear methods. This will also be useful to generate enriched monolithic dictionary for convex hull based approximation.

Conclusions and Perspectives

Conclusions

In Chapter 2, the robustness of low-rank hypothesis in the context of contact mechanics was analyzed. The premise of this analysis was based on the local nature of the contact pressure field which leads to an inherent lack of linear separability. As the underlying physics exhibited a high intrinsic dimensionality, reconstructions performed using ROMs had a moderate error, as expected. The robustness of contact pressure reduced bases computed using SVD and CPG methods were shown to have a slower decay of error (i.e. higher values of compactness) than the displacement field. The lack of low-rank structure had also negatively influenced the generalization ability of the reduced basis.

The linear inseparability of the contact pressure makes it harder to find a robust lowdimensional structure in the framework of linear dimensionality reduction. In Chapter 3, it is shown that a large number of snapshots can provide some improvements to the approximate solution subspace that is exploited by reduced models. This is achieved by using large dictionaries of snapshots, and the reconstructed solution is restricted to a sparse linear combination of these snapshots. It is seen that dictionary-based approximations are quite useful for contact problems with large changes in contact position like the ironing problem, as local effects are quite significant in these problems.

Two dictionary-based approximation strategies were devised. The first is based on the active-set approach implemented greedily, and the second is based on exploring the convex-hull of the snapshots to reconstruct the solution under the condition that the feasible (non-penetrative) region of the displacement space is convex. The computational time for these methods seems to be primarily dependent on the number of iterations for convergence, which may or may not change significantly with the size of the dictionary depending on the contact problem. In terms of reconstruction accuracy, the greedy active-set approach with separate dictionaries for displacement and contact pressure outperforms the convex-hull approximations with monolithic dictionaries. Moreover, the convex-hull approximations could not be applied to non-convex problems such as the ironing problem.

Non-linear dimensionality reduction aspects are particularly interesting for a problem lacking low-rank structure and are explored in Chapter 4. As the computational costs of generating a large dictionary of snapshots can be very expensive, non-linear methods are used to reduce the need for a large number of contact pressure snapshots. Dynamic Time Warping (DTW) based non-linear interpolations were used to synthesize snapshots that were not present in the training set, thereby enriching the dictionary without paying the large cost of full simulations. This is found to be beneficial in cases of large variations in contact position.

Perspectives

1. The problems arising from linear inseparability of contact pressure have been addressed using over-complete dictionaries and DTW-based enrichments in this thesis. However, distance functions that are used to compute the non-penetration conditions (i.e. the inequality constraints) are likely to have similar local effects and consequent linear inseparability issues. Also, the current methods of computing affine decompositions of distance functions (like EIM) are based on low-rank approximations requiring linear separability. Therefore, similar strategies might be required for the efficient construction of the reduced inequality constraint operators.

2. Random sketching methods were applied only to the convex hull approach as the method was free of inequality constraints. It would be interesting if the random sketching could be applied to other dictionary-based strategies as well, including the greedy active-set method.

3. The dictionary method based on the convex hull approach is currently built for monolithic dictionaries, which forces the reduced coefficients to be the same for both primal and dual variables. An approach with two uncoupled dictionaries must be explored for potential improvements in reconstruction quality. Also, an efficient sampling strategy for the offline phase is necessary, where snapshots placed at corners and edges of the presumed low-dimensional convex subset are sampled with priority.

4. As DTW-based non-linear interpolation was found useful in particular cases with large variations of contact positions, application of other non-linear methods like kPCA/kPOD should be explored. These methods might be more relevant for threedimensional contact problems with two-dimensional contact surfaces, as multidimensional DTW is limited to structured meshes.

5. Currently, non-linear interpolation is only applied to contact pressure. It would be interesting to see if non-linear interpolation of displacement snapshots provides any advantages. DTW-based enrichment of displacement snapshots, along with the contact pressure, might be beneficial for creating monolithic dictionaries which need the same number of primal and dual snapshots. Title: Low-rank and sparse approximations for contact mechanics Keywords: model order reduction, contact mechanics, dictionary-based approximation, nonlinear dimensionality reduction Abstract: Typical strategies for reduced modelling of contact mechanics use low-rank approximations. The underlying hypothesis is the existence of low-dimensional subspace and subcone for the displacement and contact pressure fields, respectively. However, the contact pressure exhibits a local nature, as the contact position can vary with parameters like loading or geometry.

The adequacy of low-rank approximations for contact mechanics is studied and alternative routes based on sparse regression techniques are explored. It is shown that the local nature of contact pressure leads to a loss in linear separability, thereby limiting the accuracy of low-rank methods.

Subsequently, approximations using overcomplete dictionaries to mitigate the inseparability issues is investigated. Two strategies are devised using sparse approximation methods. One is based on an active-set approach where the dictionary elements are selected greedily, and another is based on convex hull approximations where the non-penetration constraints need not be explicitly enforced. Lastly, the non-linear dimensionality reduction framework is explored. The snapshot set computed in the training phase is enriched at a low complexity using non-linear interpolations, thereby reducing the burden of creating overcomplete dictionaries in the offline phase.

3. 6

 6 Greedy active-set reconstruction of contact pressure for certain parametric instances of the Hertz problem . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Reconstruction errors, primal rank, number of iterations and computation time vs. training set size for the Hertz problem . . . . . . . . . . . . . . . 3.8 Detailed primal and dual reconstruction errors for the Hertz problem . . . 3.9 Sparsity of dual reduced dofs for the validation set of Hertz problem . . . 3.10 Greedy active-set reconstruction of contact pressure for certain parametric instances of the ironing problem. . . . . . . . . . . . . . . . . . . . . . . . 3.11 Sparsity of dual reduced dofs selected by greedy active-set method for the ironing problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Reconstruction errors, primal rank, number of iterations, and computation time vs. training set size for the Ironing problem . . . . . . . . . . . . 3.13 Illustration of non-penetrating property of convex hull . . . . . . . . . . . 3.14 Illustration of the convex subset hypothesis test with a point mass moving in a convex feasible region. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Testing the convex subset hypothesis for the elastic rope-obstacle problem 3.16 Convex hull reconstruction errors for elastic rope-obstacle problem. . . . . 3.17 Testing the convex subset hypothesis for the Hertz problem . . . . . . . . 3.18 Convex hull reconstruction errors for the Hertz problem. . . . . . . . . . . 3.19 Mean number of iterations and computation times for Hertz problem using convex hull exploration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Demonstration of the non-convex nature of the ironing problem. . . . . .4.1 Visualization of contact pressure non-linearity in ironing problem . . . . . 4.2 DTW illustration: Vectors λ 1 and λ 2 and their DTW alignments . . . . . 4.3 DTW illustration: λ and x vs. warped coordinate ζ . . . . . . . . . . . . 4.4 Alignments and warped snapshots computed by DTW between two contact pressure snapshots of ironing problem . . . . . . . . . . . . . . . . . . 4.5 Interpolation along DTW alignments and in warped space . . . . . . . . . 4.6 Interpolations of the illustrative example for α = 0.5 in the spatial indices (corresponding to x) and the warped space ζ . . . . . . . . . . . . . . . . 4.7 Linear and DTW-based interpolations of contact pressure for the elastic rope-obstacle problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Linear and DTW-based interpolations of contact pressure for the Hertz problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Interpolation of contact pressure snapshots using DTW for the ironing problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Comparison of the compactness for original and warped snapshots . . . .

2. 1 1

 11 Truncation tolerances and ranks for both bases for the rope-obstacle problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Truncation tolerances and ranks for both bases for Hertz problem. . . . . 3.1 Characteristics of the two dictionary-based approximation methods for contact mechanics problems . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Observed values of various compliances for the elastic rope-obstacle problem with parameter γ 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Input variables values and the observed DTW alignments for the illustrative example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Nomenclature u, v Continuous/Discretized Displacement λ, η, λ, η Continuous/Discretized Lagrange Multiplier (Contact pressure) K Feasible region of an inequality constrained problem L Lagrangian functional V, W + Displacement solution space and Contact pressure solution cone a, f Bilinear and Linear forms of energy functional b, d Bilinear and Linear forms of distance function in weak form Displacement and contact pressure shape functions r(•) Residual of a discretized problem • h Discrete version of the entity • C, g Discretized operators of bilinear and linear terms k and g K, f Discretized operators of bilinear and linear terms a and f K mono Operator to generate the residual with monolithic dictionary Γ c , Γ /c Surface in contact, and out of contact (gap) d, d Element in a dictionary, Leave-one-out candidate in a dictionary • r Reduced space of • B An orthogonal basis of a given snapshot matrix or its randomized combination D Dictionary, Monolithic Dictionary R Matrix containing random entries between [0, 1] S tr Training set snapshots N r u , N r λ Number of reduced dofs for u and λ Π • Projection operator on the subspace/subcone of basis • D p , D d Primal and Dual Dictionaries Φ, Θ Primal and Dual reduced bases φ, θ Primal and Dual reduced bases vectors/functions • Reduced operators and dofs ε CHLS Convex hull least square error for the test of convex subset hypothesis ε CrPen Cross penetration error in convexity test of the feasible region C(m) Compactness of rank-m G(m) Generalization ability of rank-m basis H(n) Nested error between n-th and n + 1-th level nested training set S(m) Specificity of rank-m basis α, γ Coefficient vector multiplying with columns of a matrix/basis/dictionary x Vector/Signal to be approximated α Parameter for interpolating in high-dimensional space • int Quantity • interpolated in warped space D adapt p Dictionary adapted/enriched using non-linear interpolations ϕ Non-linear transformation to a higher-dimensional warped space/domain ζ Co-ordinate in the warped space • Quantity • in warped space I k DTW alignments computed for a vector λ k 0, Matrix or Vector consisting of zeros or ones of appropriate size •[:, S], •[S] Slicing of a vector/matrix • using indices in set S (Python style) δ Truncation tolerance for an orthogonal basis ε Machine precision, Tolerance on residual (z) - z if z < 0, 0 otherwise (z) +,τ z if z > τ , 0 otherwise • † Moore-Penrose pseudo-inverse

Figure I :

 I Figure I: Decision-making process using numerical simulations. Additional time delays can be caused by communication and simulation. In an ideal situation: Communication delay + Simulation time Manufacturing time

Figure 1 . 1 :

 11 Figure1.1: Kinematic description of mechanical problems with possibility of contact. Points x 2 will possibly come into contact with point x1 , and hence considered as a contact pair. The distance between the two points in the initial and deformed configuration is given by g(µ, v)(x 2 ) and g(µ, v)(x 2 )k(µ, v; v)(x 2 ) respectively

Figure 1 . 2 :

 12 Figure 1.2: Schematic of node-to-node contact pairing for conforming surface meshes

( a )Figure 1 . 3 :

 a13 Figure 1.3: Schematics of node-to-segment methods for finding contact pairs

Algorithm 2 . 1

 21 Online phase 1: Input: Queried value of parameter µ 2: Given: Primal basis Φ and dual basis Θ Reduced operators K, f (possible to build offline) 3: Initiate boolean array A with one random element set to True. 4: while u not converge do 5:

  , γ) (a) Deformation u(x, γ) and obstacle g(x, γ) snapshots Contact pressure snapshots

Figure 2 . 1 :

 21 Figure 2.1: Sample snapshots for elastic rope-obstacle problem with γ ∈ P tr ⊂ [-0.5, 0.5].

Figure 2 . 2 :

 22 Figure 2.2: Primal reconstruction (H 1 ) errors for points in and outside the training set for the rope-obstacle problem.

Figure 2 . 3 :

 23 Figure 2.3: Decay of singular values (cumulative and normalized) for the training set of the elastic rope-obstacle problem. A similar figure is also given in Reference [48, Figure 4]

Figure 2 . 4 :

 24 Figure 2.4: Sparsity pattern of dual snapshots for the rope-obstacle problem. Each row corresponds to a snapshot

2. 4 . 1

 41 Hertz problemThis section deals with the reduced model of a Type-2 contact problem, namely the Hertz contact problem, of two half-cylinders loaded against each other, as shown in Fig.2.5a. A loading parametrization problem is solved, where the imposed displacement d ∈ (0, 0.3) is the parametric space, with R 1 = R 2 = 1.0. The contact pressure snapshots for various imposed displacement values are shown in Fig.2.5b.

  d = 0.275 d = 0.25 d = 0.225 d = 0.2 d = 0.175 d = 0.15 d = 0.125 d = 0.1 d = 0.075 d = 0.05 d = 0.025 (b) Contact pressure snapshots

Figure 2 . 5 :

 25 Figure 2.5: Hertz problem: two half cylinders loaded against each other and the resulting contact pressure snapshots. Displacement d, imposed on the top cylinder, is treated as parameter in the reduced model.

Figure 2 . 6 :

 26 Figure 2.6: Primal reconstruction (H 1 ) errors for the Hertz problem.

Figure 2 . 7 :Figure 2 . 8 :

 2728 Figure 2.7: Sparsity pattern of dual snapshots for Hertz problem

Figure 2 . 9 :

 29 Figure 2.9: Validation metrics for Hertz problem. The bars indicate ±1σ interval of corresponding metrics

Figure 2 . 10 : 5 Figure 2 . 11 :

 2105211 Figure 2.10: Ironing problem: Iron block is pressed against the flat slab by a displacement d y and moved horizontally. The horizontal displacement of the iron block d x ∈ [0, L] is treated as the parameter in the reduced model

Figure 2 . 12 :

 212 Figure 2.12: A typical sparsity pattern of dual snapshots for ironing problem

Figure 2 . 13 :

 213 Figure 2.13: Compactness for reduced bases of ironing problem.

Figure 2 . 14 :

 214 Figure 2.14: Projection error of nested level n + 1 snapshots on dual RB of level n.

Figure 2 . 15 :

 215 Figure 2.15: Validation metrics for the ironing problem (for fine mesh). The fill regions indicate ±1σ interval of corresponding metrics.

Figure 3 . 1 :

 31 Figure3.1: Illustration of low-rank vs. dictionary-based approximation. Snapshots (circles) in an n-dimensional system lie on a manifold of dimension lower than n. The low-dimensional subspace computed using low-rank methods has a dimension n. A dictionary-based approximation of the queried point (star-shaped) can be computed as linear combination of the nearest snapshots (shown by a line connecting the nearest candidates). The low-rank approximation is simply the projection of the queried point on the low-dimensional subspace.

Algorithm 3. 1 FOCUSS 1 : 2 :

 112 Inputs: D, x, α 0 Initialize k = 1 3: while α not converge do 4:

Figure 3 . 2 :

 32 Figure 3.2: Two-dimensional domain for the transport problem (3.4)

Figure 3 . 3 :

 33 Figure 3.3: Snapshots of the transport problem (3.4) for randomly sampled parametric values

Figure 3 . 4 :

 34 Figure 3.4: Reconstruction errors (bars) of the transport problem using dictionary-based approximations and number of non-zero entries (×) on training and validation sets. L indicates the target sparsity

  the active set and D d,I = D d [:, I]

Algorithm 3 . 2 2 :

 322 Greedy active-set algorithm 1: Input: Queried value of parameter µ Given: Primal basis Φ and dual Dictionary D d Reduced operators Φ T KΦ and Φ T f can be built offline 3: Initialize: k = 0, I = ∅ 4: while u and I not converge do 5:

7 : 8 : 13 :

 7813 Set λ[I] ← λ I and λ[I c ] ← 0 I c is the complementary set Compute the largest constraint violations and enrich/eliminate accordingly: if λ k ≥ 0 : then p add = arg max p CΦ u kg +,τ I ← I ∪ {p add } 9: else: p rem = arg min p λ k - [p] I ← I \ {p rem } Reconstruct u = Φ u and λ = D d λ 14: Output: u, λ d ∈ (0, 0.3). The finite element model also remains the same.

Figure 3 . 5 : 5 × 10 - 3 Figure 3 . 6 :

 35510336 Figure 3.5:A spy pattern of dual dictionary with 30 elements for the Hertz problem. The number of rows, 79, is same as the number of dofs on the potential contact surface. The dictionary columns are arranged in the increasing order of loading parameter d.

Figure 3 . 7 :

 37 Figure 3.7: Evolution of (a) Mean reconstruction (relative) errors and primal rank, (b) Mean number of iterations, (c) Mean computation time and (d) Mean computation time per iteration with training set size for the Hertz problem. The mean reconstruction error is computed by taking mean of reconstruction errors over the validation set P val . Primal and dual errors are computed using H 1 and L 2 norm, respectively. Curves correspond to different truncation tolerance δ for the primal basis, for values 10 -6 , 10 -8 and 10 -10 . Note that computation time excludes time for construction of non-linear operators.

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: Detailed primal and Dual reconstruction errors for the Hertz problem using the training sets 1,2 and 3. Validation set is common to all curves, consisting of 60 equispaced points in P. Each subfigure corresponds to a different truncation tolerance δ for the primal basis. Primal and dual errors are computed using H 1 and L 2 norm, respectively. Cross-marks × indicate points where the greedy active-set algorithm did not converge to the defined tolerance 10 -5

4 Figure 3 . 10 :Figure 3 . 11 :Figure 3 . 12 :

 4310311312 Figure 3.10: Greedy active-set reconstruction of contact pressure for certain parametric instances of the ironing problem. The dictionary nested level, dictionary elements chosen by the algorithm and the reconstruction errors are also shown.

Figure 3 . 13 :

 313 Figure3.13: Illustration of non-penetrating property of convex hull. The non-penetration constraint in (3.9) is linear, convex and is made of segments corresponding to each element of the constraints Cu = g. In a high-dimensional problem, each segment of the constraint can be thought of as a hyperplane. The convex hull, which is simply a convex enclosure of snapshots, can also be thought of as a high-dimensional set.

Figure 3 . 14 :

 314 Figure 3.14: Illustration of the Test B with a point mass moving in a V-shaped obstacle forming a convex feasible region. Two leave-one-out instances and the convex hull of the corresponding D p are shown. CHLS is a measure of distance between the left out snapshot and the convex hull of D p . The leave-one-out instance 1 has a higher CHLS than the instance 2.

Figure 3 . 15 :

 315 Figure 3.15: Test B: testing the Hypothesis 3.1 for the given training set of elastic rope-obstacle problem (3.13). Test A for convexity of feasible set is not necessary as the constraints are linear.

γ 1 ∈ 3 ) 3 . 2 :

 1332 P tr γ 1 ∈ P val |1 T α -1| O(10 -13 ) O(10 -3 ) Penetration O(10 -16 ) O(10 -4 ) Complementary Slack O(10 -15 ) O(10 -Table Observed values of various compliances for the elastic rope-obstacle problem with parameter γ 1 .

Figure 3 . 16 :

 316 Figure 3.16: Convex hull reconstruction errors for elastic rope-obstacle problem with parameter γ 1 using B ← svd(D p , δ = 10 -7 ). Crosses × indicate achieved sparsity (w.r.t. right y-axis).Primal and dual sparsity are naturally equal as the dictionary is monolithic.

Figure 3 . 17 :

 317 Figure 3.17: Test B: testing the Hypothesis 3.1 for the various training sets of the Hertz problem

Figure 3 . 18 :

 318 Figure 3.18: Convex hull reconstruction errors for the Hertz problem using different dimensions of B.The penetration in all reconstruction cases is observed to be in order of 10 -8 . The case of rank(B) = 15 and #P tr is not possible, hence not shown in (b). Primal and dual errors are computed using H 1 and L 2 norm, respectively.

Figure 3 . 19 :

 319 Figure 3.19: Mean number of iterations and computation times for Hertz problem using convex hull exploration with different dimensions of B. Note that computation time excludes time for construction of non-linear operators.

Figure 3 . 20 :

 320 Figure 3.20: Demonstration of the non-convex nature of the ironing problem. Two snapshots from training set of the ironing problem and a convex combination of the snapshots that violates the non-penetration condition

Figure 4 . 1 :

 41 Figure 4.1: Visualization of contact pressure non-linearity in ironing problem. The dimension x ∈ Γ is the coordinate on contact surface, and d ∈ P is the coordinate on the parametric domain. The arrow shown in the figure points to a direction in which snapshots are "aligned" along which the snapshots are, possibly, linearly separable.

  3), the mapping to warped coordinate ζ. The warped snapshots λ k = λ k [I k ] and the warped spatial coordinate x k = x[I k ] of dimension n are shown in Fig. 4.3. The DTW warped coordinate ζ can be simply expressed as set of indices {1, 2 . . . n}.

Figure 4 . 2 :

 42 Figure 4.2: DTW illustration: Vectors λ 1 and λ 2 and their DTW alignments (dotted). Note that the vectors are offset in y-direction for illustration purposes.

Figure 4 . 3 :

 43 Figure 4.3: DTW illustration: λ and x vs. warped coordinate ζ

Figure 4 . 4 :

 44 Figure 4.4: (a) Alignments and (b) warped snapshots computed by DTW between two contact pressure snapshots of ironing problem for parametric values µ 1 = d 1 and µ 2 = d 2 .Notice that alignments in original space are not bijective. y-offset between the two curves is given only for illustrative purpose. Values on ζ axis are not given as they do not carry a physical significance, and also because DTW computes warped space on indices of the discretization rather than physical dimension.

Figure 4 . 5 :

 45 Figure 4.5: Interpolation along (a) DTW alignments and (b) interpolation in warped space, for parametric values µ 1 = d 1 and µ 2 = d 2 . Notice that unlike in (a), the interpolation in (b) is a linear operation.

  λint and λint vs. ζ

Figure 4 . 6 :

 46 Figure 4.6: Interpolations of the illustrative example for α = 0.5 in the spatial indices (corresponding to x) and the warped space ζ

Figure 4 . 7 :

 47 Figure 4.7: Linear and DTW-based interpolations of contact pressure for the elastic ropeobstacle problem.

Figure 4 . 8 :

 48 Figure 4.8: (a) Linear and (b) DTW-based interpolations of contact pressure for the Hertz problem. To explain distorted interpolation, DTW alignments is also shown in (c).

Figure 4 . 9 :

 49 Figure 4.9: Interpolation of contact pressure snapshots using DTW for the ironing problem. Interpolated curves for α ∈ {0.25, 0.5, 0.75} are shown

Algorithm 4 .2 adapt dictionary 1 : 6 : 7 :

 4167 Input: D d , I, {α m } M m=1 α m ∈ ]0, 1[ 2: Parameters: Dist, N 3: Initialize: D adapt d = D d [:, I], i = 0, j = 0 I is the current active set 4: for i ∈ I do 5: nbrs ← indices of the least N positive entries of Dist[i, :] Initialize D int d = {•} for j ∈ nbrs do nbrs contains neighbouring indices 8:

∪

  D d [:, j], D d [:, i], α m ) Alg. D d [:, nbrs] ∪ D intd 13: end for 14: Output: D

Figure 4 . 10 :

 410 Figure 4.10: Comparison of the compactness for original and warped snapshots. Note that the curve labelled "dual cone" are the same curves from Figs. 2.8 and 2.13b, respectively. Curves labelled "dual cone dtw" show the compactness of the enriched basis.

Figure 4 .

 4 Figure 4.11: Projection error of nested level n + 1 snapshots on dictionary of level n. Note that solid dual orth and dual cone lines are similar to the ones in Fig.2.14. The dual orth points not seen in the figure drop to the order of 10 -14 (nearly machine precision)

Algorithm 4 . 3 2 : 6 : 7 :λ = D d λ 8 : 9 :

 4326789 DTW powered greedy active-set method 1: Input : Queried value of parameter µ Parameters : {α m } M m=1 , L 3: Given: Primal basis Φ and dual Dictionary D d Reduced operators Φ T KΦ and Φ T f 4: Initialize: l = 0, I = {•}, λ = 0 5: while λ not converge and l ≤ L do Compute u, λ, I using current dictionary D d Steps 4 to 12 of Alg. 3.2 Reconstruct D d ← adapt dictionary(D d , I, {α m }) Alg. 4.2 I = {•} Reset active set after adapting D d 10: l ← l + 1 11: end while 12: Reconstruct u = Φ u 13: Output: u, λ 4.3.1 Numerical results For the Hertz problem and the ironing problem training sets defined in Chapter 3, reconstruction errors with enriched dual dictionaries using Alg. 4.3 are shown in Fig. 4.12. For the ironing problem (Fig. 4.12b), the DTW-based interpolation improves the quality of reconstruction for smaller training sets. The improvement asymptotically dampens out as the training set size increases. This is natural as snapshots come quite close in large training sets and DTW-based interpolations do not generate any useful information.

Figure 4 . 12 :

 412 Figure 4.12: Comparison of reconstruction error for the Hertz and ironing problem using the original dual dictionary D d and DTW-adapted dual dictionary D adapt d

Titre:

  Approximations parcimonieuses et de rang faible pour la mécanique des contacts Mot clés : modèles d'ordre réduit, mécanique des contacts, approximations dictionnaires, réduction non-linéaire de dimensions Résumé : Les stratégies typiques de modélisation réduite de la mécanique des contacts utilisent des approximations de rang inférieur. L'hypothèse sous-jacente est l'existence d'un sousespace et d'un sous-cône de faible dimension pour les champs de déplacement et de pression de contact, respectivement. Cependant, la pression de contact présente un caractère local, car la zone de contact peut varier en fonction de paramètres tels que le chargement ou la géométrie. Dans cette thèse, la pertinence des approximations de rang faible pour la mécanique des contacts est étudiée et des voies alternatives fondées sur des techniques de régression parcimonieuses sont explorées. Il est montré que la nature locale de la pression de contact conduit à une perte de la séparabilité, limitant ainsi la pré-cision des méthodes de rang faible. Par la suite et pour atténuer les problèmes d'inséparabilité, les approximations utilisant des dictionnaires contenant un grand nombre d'atomes sont étudiées. Deux stratégies sont établies à l'aide de méthodes d'approximation parcimonieuses. L'une est fondée sur une approche d'ensemble actif où les éléments du dictionnaire sont sélectionnés de manière gloutone, et l'autre sur des approximations d'enveloppes convexes qui éliminent le besoin d'appliquer explicitement des contraintes de non-pénétration. Enfin, les méthodes de réduction non linéaires sont explorées. Elles permettent d'enrichir le dictionnaire en utilisant des interpolations nonlinéaires, réduisant ainsi le coût de calcul associé à la construction du dictionnaire dans la phase hors-ligne.

Table 2 . 1 :

 21 Truncation tolerances and ranks for both bases for the rope-obstacle problem. ε

			Primal basis	Dual basis
			(POD)	(CP-greedy)
	Full	Tolerance Rank	ε 11	ε 11
	Truncated	Tolerance Rank	ε 11	2 × 10 -1 7
	indicates numerical precision

  .3. However, an alternative definition given by(2.4) is used here. The new definition allows extending the concept of compactness to CPG reduced bases, where equivalents of singular values are undefined. Also, in case of an orthogonal basis, this definition is equivalent to the normalized singular values used in Fig.2.3. tr is a matrix containing (either primal or dual) snapshots in P tr arranged column-wise. Π Ψ[:m] (S tr ) is a projection operators that projects each column in S tr on the first m vectors of the basis Ψ (indicated in Python notation by Ψ[: m]). The snapshot matrix S tr , projection operator Π and the basis Ψ are generic symbols and are used in the following combinations, indicated by the following labels:

		C(m) :=	S tr -Π Ψ[:m] (S tr ) F ||S tr || F	(2.4)
	where S Label in figures Snapshots in S tr Projection operator Π	Basis Ψ
	primal orth	Primal	Orthogonal	POD primal basis Φ
	dual orth	Dual	Orthogonal	POD dual basis Υ
	dual cone	Dual	Cone	CPG dual basis Θ

Table 2 . 2 :

 22 Truncation tolerances and ranks for both bases for Hertz problem. ε indicates numerical precision

			Primal basis	Dual basis
			(POD)	(CP-greedy)
	Full	Tolerance Rank	ε 12	ε 12
	Truncated	Tolerance Rank	ε 12	5 × 10 -2 8

Table 3 . 1 :

 31 Characteristics of the two dictionary-based approximation methods for contact mechanics problems

		Greedy active-set	Convex hull approxima-tion
	Primal Basis type	Truncated low-rank	Dictionary
	Dual Basis type	Dictionary	Dictionary
	Monolithic dictionary	No	Yes
	Iteration method	Fixed point	Fixed point
			Convex	combinations
	Constraint enforcement	Active-set	of	non-penetrating
			snapshots
		Greedy activation and
	Sparsity induction	disactivation of dual	Non-negative FOCUSS
		dofs
	Additional assumptions	None	Convexity of the feasible region

Table 4 . 1 :

 41 Input variables values and the observed DTW alignments for the illustrative example

The terms hyperplane and hypersurface refer to a plane and surface in dimensions higher than 3 respectively.

column-cone is the subset of the column-space restricted to non-negative coefficients.

colsp indicates column-space

A convex combination is a linear combination with non-negative coefficients that sum up to unity

Convex hull of a set of points (or snapshots) is the set containing all convex combinations of the points. In other words, it is also the smallest convex set that envelopes the set of points.

usually the term "linear" is dropped in unambiguous cases

Though the DTW algorithm was originally intended for temporal data, here it will be used to warp spatial dimension.

Acknowledgement

List of Algorithms A | Spurious penetrations

In the greedy active-set method, the snapshots of the dual dictionary are selected based on the violations of the non-penetration constraints. The dictionary vector with the highest correlation with the current state of penetration is chosen. However, instead of zero penetration, a small value of penetration τ is allowed (see (3.5c) and (3.6)). This is done to avoid spurious selection of dictionary elements, which happens due to a truncated primal basis. Such a reconstruction example of Hertz problem for parameter value d = 0.25 is shown in Fig. A.1, where a spurious peak in the contact pressure is evident. In this case, a dual dictionary of size 30 and primal basis truncated at 10 -6 are used. The algorithm selected the appropriate snapshots from the dual dictionary, but also selected some unnecessary snapshots: the first (D d [START_REF] Lauzeral | A model order reduction approach to create patient-specific mechanical models of human liver in computational medicine applications[END_REF]) and twelfth (D d [START_REF] Yastrebov | Numerical Methods in Contact Mechanics[END_REF]) dictionary vectors. Importantly, the D d [START_REF] Lauzeral | A model order reduction approach to create patient-specific mechanical models of human liver in computational medicine applications[END_REF] vector contributes to large error in reconstruction as evident in the figure. The false selection of this snapshot is linked to the truncation of primal basis, and therefore, can be avoided by setting the value of τ same as δ. , where few dictionary snapshots were spuriously selected by the greedy active-set algorithm if non-penetration condition is applied "hardly". The falsely selected snapshots are indicated by their index in the dual dictionary. The snapshot which contributes to most of the error is D d [START_REF] Lauzeral | A model order reduction approach to create patient-specific mechanical models of human liver in computational medicine applications[END_REF], corresponds to loading parameter d = 0.01