Nonintrusive reduced order model for parametric solutions of inertia relief problems

The Inertia Relief (IR) technique is widely used by industry and produces equilibrated loads allowing to analyze unconstrained systems without resorting to the more expensive full dynamic analysis. The main goal of this work is to develop a computational framework for the solution of unconstrained parametric structural problems with IR and the Proper Generalized Decomposition (PGD) method. First, the IR method is formulated in a parametric setting for both material and geometric parameters. A reduced order model using the encapsulated PGD suite is then developed to solve the parametric IR problem, circumventing the so-called curse of dimensionality. With just one offline computation, the proposed PGD-IR scheme provides a computational vademecum that contains all the possible solutions for a pre-defined range of the parameters. The proposed approach is nonintrusive and it is therefore possible to be integrated with commercial FE packages. The applicability and potential of the developed technique is shown using a three dimensional test case and a more complex industrial test case. The first example is used to highlight the numerical properties of the scheme, whereas the second example demonstrates the potential in a more complex setting and it shows the possibility to integrate the proposed framework within a commercial FE package. In addition, the last example shows the possibility to use the generalized solution in a multi-objective optimization setting.


Introduction
Unconstrained structures are widespread in the automotive, aerospace and naval industry. As is well known, due to the singularity of the stiffness matrix, conventional static analyses cannot be performed if the system undergoes rigid body motions. At the same time, imposing dummy constraints in order to make a free-body system statically determinate leads to unrealistic reaction forces and, as a consequence, an unrealistic distribution of the internal stresses. The inertia relief (IR) method represents an attractive alternative for solving unconstrained structural problems without resorting to the more expensive full dynamic analysis [1]. The main idea is to counteract the unbalanced applied loads by a set of rigid body accelerations, the latter providing body forces which are distributed over the structure in such a way that the applied forces are equilibrated and the static analysis can be performed. The technique is available into most of the commercial finite element packages and it has been widely used by the industry in different fields [2,3,4,5,6,7,8,9,10].
The static global stiffness analysis of a body in white (BIW) is a common example that involves the computational simulation of an unconstrained configuration using the IR method. The BIW global stiffness plays a significant role in the design process of a car. An important challenge of this problem is the number of parameters (e.g. geometry, material) to be considered during the analysis of a BIW. As any change in the material or geometrical characteristics of the car components might have considerable effects on the global behaviour of the structure, the number of simulations that are required to account for the whole range of the involved parameters becomes prohibitive when classical numerical approaches are employed. As a consequence, the possibility to perform parametric studies, shape optimization or inverse identification in the context of a BIW remains a challenge.
A way to circumvent this issue and to reduce the computational complexity of parametric problems is to employ a reduced order model (ROM). In the last two decades, researchers from the most diverse areas of science and engineering have developed different ROM techniques, with the common goal of finding low-order models, described by a reduced order basis, which are able to capture the essential behaviour of a complex system. Well known computational approaches based on this idea are Krylov-based methods [11], the reduced basis method [12] and the proper orthogonal decomposition (POD) [13,14,15]. These techniques, known as a posteriori methods, first solve the full-order problem for a suitably chosen set of parameters, providing a set of snapshots of the solution. This step, usually referred to as the offline stage, is used to extract the most relevant characteristics of the solution. Then, during the online phase, the solution for any new parameter value can be expressed as a linear combination of the previously computed basis functions.
A valuable alternative is represented by the proper generalized decomposition (PGD) method [16,17,18,19], which is an a priori approach. The main idea behind the PGD method is to consider the parameters as extra coordinates of the problem, increasing the dimensionality of the problem at hand, and to assume that the solution of the highdimensional problem can be approximated by a separable function. During the offline stage, usually performed by employing high performance computing resources, the PGD algorithm computes on-the-fly a set of basis functions, usually called modes. The PGD approximation depends explicitly on the parameters, so it represents a computational vademecum containing the solution for every possible combination of the parameters. In the online stage, the solution can be particularized for any set of the parameters in real-time. The method has been tested in the most diverse fields, such as flow problems [20,21,22,23,24], thermal problems [25,26,27], solid mechanics [28,29], fracture mechanics [30,31], elastic metamaterials and coupled magneto-mechanical problems [32,33].
Despite the wide range of problems where the PGD has shown its potential, the application to geometrically parametrized problems remains particularly challenging, due to the difficulty to obtain a separable expression of the discrete problem. Previous works that dealt with parametric shapes are usually limited to simple geometric dependence [17,34,35,36,37]. Other authors proposed a technique based on the idea that a parent domain can be associated to the parametric domain in order to introduce the parametric dependency on the geometry in the governing equations [34,38]. More recently, another approach was proposed [39] in which the control points characterising the NURBS curves or surfaces used in CAD representation are defined as the geometric parameters of the problem.
One of the main drawbacks of the original PGD approach, when compared to other a posteriori approaches, is the intrusivity of its implementation, which precludes its wider application in an industrial context, where commercial software are usually employed. This limitation has motivated the development of nonintrusive implementations of the PGD rationale for solid [40] and fluid [41] mechanics problems. Also motivated by the goal of achieving the nonintrusive applications of the PGD, the Encapsulated PGD Toolbox has been recently developed by Díez et al. [42,43]. The toolbox consists of a collection of PGD-based routines that are able to perform algebraic operations for multidimensional tensors. One relevant advantage of the toolbox is that each routine is encapsulated and can be used as a black box, enabling the nonintrusive coupling with commercial FE packages. This feature is of major importance for the application of the PGD method in an industrial setting.
This work proposes a nonintrusive PGD-IR method for the solution of an unconstrained structure characterized by material and/or geometric parameters. The proposed parametric IR approach involves a "cascade" application of the PGD method in order to solve three sequential parametric problems, where the parametric solution of one problem is taken as the input of the next parametric problem. In order to automate the process, an ad-hoc solver was implemented that makes use of the Encapsulated PGD Toolbox. In the present work, a nonintrusive interaction between the external commercial finite element (FE) software MSC-Nastran and an in-house code implemented in Matlab is considered.
The structure of the remainder of the paper is as follows. Sec. 2 presents the problem statement in terms of an elastodynamic boundary value problem. Sec. 3 briefly reviews the idea behind the IR method for a non-parametric problem. The proposed PGD-IR approach is presented in 4, where the "cascade" application of the PGD approach is proposed to solve the PGD-IR problem. Also, the algebraic approach to deal with geometric parameters is detailed. In Sec. 5 two numerical examples are used to show the potential of the proposed method. In the first example, a simple linear elastic 3D structure with one material and one geometric parameter are considered to underline the main properties of the developed ROM. The second example uses a more realistic case of a dummy car to demonstrate the nonintrusive interaction with the commercial FE software MSC-Nastran. A multiobjective optimization study is shown which proves that the method can be employed as a fast and reliable tool to guide the designers in the intricate decision-making procedure. Finally, Sec. 6 summarizes the conclusions of the work that has been presented. The strong form of the elastodynamic problem using the classical Voigt notation [44] can be written as where u is the displacement field,ü denotes the acceleration, σ is a vector containing the extensional and shear stress components of the Cauchy stress tensor, b is a external body force vector, T is the final time of interest, E is a matrix accounting for the normal direction to the boundary, u D and t N are the prescribed displacement and traction vectors on the Dirichlet and Neumann boundaries respectively and u 0 and v 0 are the initial position and velocity respectively. In three dimensions, the matrix operator ∇ S and the matrix E are given by with n being the outward unit normal vector to ∂Ω. For a linear elastic material, the generalized Hooke's law expresses a linear relation between the stress vector, σ, and the strain vector, ε, namely where ε := ∇ S u and D is a symmetric positive definite matrix depending upon the Young modulus, E, and the Poisson ratio, ν. In three dimensions The weak formulation of the strong form of Eq. (1) reads as follows: given u D on Γ D and

Spatial discretization
A partition of the domain Ω in a set of n el disjoint elements Ω e is considered. Following the classical isoparametric framework, the approximation of the displacement field is defined in a reference element, Ω, with local coordinates ξ, as where U j are nodal values, N j are polynomial shape functions of order p in the reference element and n en is the number of nodes per element. The so-called isoparametric mapping, given by is employed to establish the relation between the reference element, Ω, and a generic physical element, Ω e , with nodes {x j } j=1,...,nen . Employing the isoparametric mapping of Eq. (7), the element and boundary integrals are mapped to the reference space. By using the approximation of the displacement field given by Eq. (6) and selecting the space of weighting functions to be equal to the space spanned by the interpolation functions, the following system of ordinary differential equations is obtained As usual in a finite element context, the global mass matrix M, the global stiffness matrix K and the global forcing vector F are obtained by assembling the elemental contributions given by In the above expressions B e := (J e ϕ ) −1 ∇ S N is the strain-displacement matrix, J e ϕ is the Jacobian of the isoparametric mapping, J e Γ is the Jacobian of the restriction of the isoparametric mapping to an element face and the matrix N , in three dimensions, is given by

The inertia relief method
In this section, a short review of the IR method is presented. As already mentioned, the IR method [1] is available in many commercial FE packages and it is widely used in industry to solve unconstrained structural problems without resorting to the more expensive full dynamic analysis.
When constant unbalanced external loads are applied to an unconstrained structure, the whole system undergoes a steady-state rigid body acceleration in each free direction and, due to the mass of the system, inertial forces are generated that deform elastically the body. The trajectory as rigid body of the system (as if it was infinitely stiff) is described by a displacement field U r (t) such that KU r (t) = 0, and therefore also KU r (t) = KÜ r (t) = 0. The global displacement U has to be complemented with the elastic deformation, namely The elastic deformation field U e is important to analyze the internal stresses created by the motion and also to assess other quantities of interest like the torsional stiffness, which is one of the aims of this paper.
The key idea of the IR method is to compute U e solving a static problem where the forces F eq are equilibrated, that is the resultant forces and moments are zero. Despite matrix K is singular, the fact that F eq is equilibrated guaranties that system (12) is solvable with a family of infinite solutions, all equivalent up to a rigid body motion.
Isostatic constrains (as many as rigid body modes, 3 in 2D and 6 in 3D) have to be set to compute one of these solutions (they all produce the same strains and stresses).
The idea of the inertia relief is to compute the equilibrated forces as noting thatÜ r is the rigid body mode (recall that KÜ r (t) = 0) such that F eq is equilibrated. Thus, equation (12) is derived from (8) assuming thatÜ e = 0 (which stands under a constant load, and therefore constant acceleration).
The rigid body acceleration vectorÜ r can be expressed as a linear combination of the (6 in 3D) rigid body modes, namelyÜ where each column rigid body transformation matrix Φ corresponds to one of the n r rigid body modes (n r = 6 in three dimensions) and the coefficient vector α is seen as containing the acceleration of each of the rigid body modes. Introducing the expression of Eq. (14) in Eqs. (13) and (12), and pre-multiplying by Φ T , the following equation is obtained which guaranties that the right-hand side term in Eq. 12 is an equilibrated system of forces (sum of forces and sum of moments equal to zero). It is worth noting that Φ T K = 0 because the eigenmodes are mutually orthogonal and the eigenvalues (frequencies) associated to the rigid body modes are zero. The vector of unknown accelerations α providing the equilibrated forces F eq is therefore computed by solving the system where Φ T MΦ and Φ T F are a reduced 6 × 6 mass matrix and a reduced 6 × 1 load vector, respectively. To completely define the rigid body acceleration vectorÜ r in Eq. (14) it is only necessary to compute the rigid body modes of the structure, that is the 6 columns of matrix Φ. They correspond to the kernel of the global stiffness matrix, so they are computed as the solution of KΦ = 0. To this end, the set of indices corresponding to the n d = d × n n degrees of freedom, with n n being the number of mesh nodes, is partitioned into the reference set s and the remaining set l. To simplify the notation, and without loss of generality, the set s is assumed to correspond to the last n r degrees of freedom. The system of equations to obtain the rigid body modes is then written as As K ll is symmetric and positive definite, the degrees of freedom of the rigid body modes corresponding to the l set can be expressed in terms of the degrees of freedom of the rigid body modes corresponding to the reference set, namely A natural assumption consists of choosing Φ s = I nr , where I nr denotes the identity matrix of dimension n r × n r , so that each column of the matrix Φ s represents a unit translation or rotation in the direction of the corresponding reference degrees of freedom. With all these premises, the relative elastic displacement U e in Eq. (12) is computed. It is worth noting that, in the IR framework, displacements are measured relative to the moving reference set of degrees of freedom s, which is subjected to a constant acceleration and undergoes infinite displacements. As a consequence, the rigid body displacement U r is not of interest and can be eliminated from the solving equation. Finally, the system to be solved to compute the relative elastic displacement, which in the remainder is simply referred to as U, reads Imposing a zero displacement in the degrees of freedom of the s set, U s = 0, ensures that the following system is solvable and provides the required relative elastic displacement at the degrees of freedom of the l set, where The IR method can be summarized in three steps as depicted in Fig. 1.

The parametric inertia relief method 4.1 Problem definition
The IR problem is now extended to the case of an unconstrained structure characterised by parametric properties. Let us introduce a set of n p material or geometric parameters denoted by µ = [µ 1 , µ 2 , . . . , µ np ] T ∈ M ⊂ R np . The set of parametric domains M is defined as the Cartesian product of a predefined interval for each one of the parameters, namely M := M 1 × M 2 × · · · × M np , with µ j ∈ M j for j = 1, . . . , n p . The semi-discrete system of Eq. (8) for a parametric problem can be written as In order to solve Eq. (21) with the IR method, three parametric steps have to be performed, following the rationale of the IR method described in Sec. 3 for the non-parametric case. The first step consists of computing the rigid body modes as Second, the vector of accelerations is given by Finally, the relative elastic displacement is computed as In Eqs. (22) to (24), µ is treated as a set of additional independent variables (or parametric coordinates), instead of problem parameters. As a consequence, the generalized solution of the three equations depends explicitly on the parameters and takes values in the multidimensional domain D = Ω × M. Standard numerical methods (e.g. finite elements, finite volumes, finite differences) would require the solution of each step of the IR method in the high dimensional space D, which is not feasible in practical problems. In this work, the PGD is proposed as a ROM able to circumvent the so-called curse of dimensionality and to provide the generalized solution of the parametric IR problem.

Cascade application of the encapsulated PGD approach
The goal of this section is to solve the parametric IR problem by means of the PGD technique. Following the standard PGD rationale, let us assume that the solution U(µ) of Eq. (21) can be approximated by a linear combination of an a priori unknown number N U of terms (or modes), namely Each PGD mode i is given by the product of a spatial term, U i , defined on the discretized space Ω and a set of parametric functions u i j (µ j ) depending, in a separated form, on each parameter µ j , for j = 1, 2, . . . , n p . The spatial term, U i , is a vector of the size of the finite element displacement vector, whereas each parametric dimension µ j is discretized with n j points with coordinates µ p j j , where p j = 1, 2, . . . , n j . The spatial and parametric modes are usually normalized and the amplitude of each mode, β i U , indicates the relevance of the i-th mode to the final separated solution.
In order to compute the terms of the summation in Eq. (25), the PGD solver typically employs a greedy approach. Assuming that the previous n−1 modes are known, the greedy algorithm computes sequentially the n-th term given by the spatial mode U n and the parametric terms u n j (µ j ) for j = 1, 2, . . . , n p . The enrichment process automatically stops when a user-defined level of accuracy is reached, that is when the amplitude β n U of the last term is smaller than a user defined tolerance. Since the unknown spatial and parametric terms U n and u n j (µ j ) are multiplying, the problem of computing the n-th term in Eq. (26) is nonlinear. More precisely, it is a nonlinear least-squares problem defined to find the best rank-one approximation (meant as the product of sectional functions) of the unknown term U n np j=1 u n j (µ j ). As commonly done in a PGD context, an alternated direction scheme is applied, which consists in solving the problem separately for each unknown function, assuming that all the others are known, until a stationary solution is reached. It is worth emphasizing that despite a non-linear problem needs to be solved to obtain each PGD mode, the computational cost of the problem increases linearly with the number of introduced parameters, making the solution of high-dimensional problems affordable. Recently, Díez et al. [43] developed the Encapsulated PGD Toolbox, which is a collection of PGD-based algorithms able to perform algebraic operations (e.g. product, division, storage, compression, solving linear system of equations, etc.) for multidimensional data represented in a discretized tensorial separated format. The main advantage of the library (freely available at https://git.lacan.upc.edu/zlotnik/algebraicPGDtools.git) is that each routine is encapsulated, meaning that it can be used as a black box. This is particularly attractive for the end user and it facilitates the interaction with commercial software.
To illustrate the idea behind the encapsulated PGD Toolbox, Fig. 2 describes the structure of the encapsulated-PGD routine that solves parametric linear systems of equations.
In a straightforward way, the same structure can be extended to other arithmetic operators. Shortly, given an algebraic linear system of equations A(µ) x(µ) = b(µ) depending on the set of parameters µ, the toolbox is able to return an explicit description of x(µ), also called computational vademecum, containing the solution for every possible combination of the parameters. The only requirement to employ the encapsulated PGD approach is to pre-process the input quantities, i.e. the parametric matrix A(µ) and vector b(µ), such that they are expressed in a PGD separated form. Given the input data, the user only needs to employ the encapsulated PGD in the offline stage, to obtain the PGD approximation by means of the above mentioned greedy algorithm and alternate direction scheme. The output consists of the sought computational vademecum and, during an online stage, the solution can be evaluated in real time for any set of parameters at a negligible computational cost. The PGD-IR approach proposed in this work makes use of the encapsulated PGD toolbox. The three parametric IR Eqs. (22) to (24) are solved sequentially, in the sense that the solution of each equation is directly needed in the next one. Consequently, as depicted in Fig. 3, a cascade PGD scheme can be employed, in which the output of each step, obtained in a separated format by simply calling the encapsulated PGD linear solver, can be directly used as an input of the next one, until the final solution of the global problem is computed. It is worth noting that, in order to use the toolbox, the user has to provide a separable representation of the input data. In particular, the stiffness K(µ) and mass M(µ) matrices must be written as where the spatial terms are K i ∈ R n d ×n d and M i ∈ R n d ×n d , whereas the parametric terms are k i j (µ j ) ∈ R n j and m i j (µ j ) ∈ R n j , for j = 1, 2, . . . , n p . Similarly, the input nodal force vector F(µ) must be written as with F i ∈ R n d and f i j (µ j ) ∈ R n j . In the above expressions N K , N M and N F are the number of modes required to produce a separable approximation of K(µ), M(µ) and F(µ) respectively.
It is important to underline that it is not always trivial to find a separated representation of the input data, as given by equations (27) and (29), especially when geometric parameters are considered in the problem. This issue will be addressed in the next section.
Finally, in order to solve the parametric equations of the PGD-IR approach, steps 2 and 3 depicted in Fig. 3 require some extra operations between the parametric objects, such as products, additions or compression. These operations can be easily performed by the Encapsulated PGD toolbox. Remark 1. The first two steps of the PGD-IR shown in Fig. 3 are parametric problems only if geometric parameters are considered, because, by definition, the rigid body modes of a structure do not depend on the material properties.

Geometric parameters: a nonintrusive algebraic approach to separate input quantities
The extension of the proposed nonintrusive PGD framework to geometrically parametrized problems represents a challenging task. This is due to the fact that, if geometric parameters are introduced in the problem, it is not trivial to find separable representation of the input quantities.
If a closed form separated expression of the stiffness and mass matrices is sought, the weak form of the problem must be modified to account for the parametric geometry. A common approach consists of formulating the problem in a reference domain, leading to several limitations that are briefly discussed in Appendix A. The most important limitation in the context of the current work is that the implementation based on a reference domain requires access to the code, precluding its application in an industrial framework, where commercial codes are typically employed.
In this section a nonintrusive algebraic approach is proposed, which is able to deal with general geometric parametrizations. The main idea is to perform a sampling of the parametric matrices and to express them in a separated format. The approach requires the computation of the parametric matrices for different values of the geometric parameters, whilst maintaining the connectivity matrix of the FE mesh. To this end a mesh morphing approach is adopted in this work. Every time a sampling of the parametric matrices is required, an initial mesh is deformed according to the geometric parameters and the global stiffness and mass matrices are computed. It is worth noting that this approach can be easily integrated in commercial packages that are equipped with a mesh morphing capability. Alternatively, the user can define the preferred mesh morphing approach and produce a set of meshes to be imported in the preferred FE software. It is also worth mentioning that the sampling does not require the solution of the FE system of equations as only the global stiffness and mass matrices are of interest for the proposed PGD-IR approach. Once the set of global stiffness and mass matrices is available, they are expressed in a separated format using the encapsulated PGD toolbox.
To illustrate the proposed nonintrusive approach, let us consider the stiffness matrix K ∈ R n d ×n d , depending on n p parameters µ = [µ 1 , µ 2 , . . . , µ np ] T ∈ M ⊂ R np . The parametric dimension µ j ∈ M j , for j = 1, 2, . . . , n p , is discretized using n j points with coordinates µ p j j , where p j = 1, 2, . . . , n j . The full-order sampling of the parametric matrix consists of evaluating K(µ) in the set of n tot points used to discretize the parametric domain M = M 1 × M 2 × · · · × M np , where n tot = np j=1 n j . Each point is characterized by its sectional indices (p 1 , p 2 , . . . , p np ), which are duly sorted by using a linear array index i such that Note that the association between the multi index (p 1 , p 2 , . . . , p np ) and the index i is also obtained by updating i = i + 1 inside n p nested loops, with no need to use explicitly expression (30). Employing the association between the multi-index (p 1 , p 2 , . . . , p np ) and the linear index i, the parametric stiffness matrix K(µ) can be written as where F p 1 ,p 2 ,...,pn p is such that F p 1 ,p 2 ,...,pn p (µ p 1 1 , µ p 2 2 , . . . , µ pn p np ) = 1 and it is equal to zero for any other values of the discrete indices p j . Using the linear indexing i introduced in Eq. (30), Eq. (31) becomes where K i = K(µ p 1 1 , µ p 2 2 , . . . , µ pn p np ), and k i j (µ p l j ) = δ p l ,p j for any p l = 1, 2, . . . , n j , while p j is given by i as defined in Eq. (30). Finally, Eq. 32 represents the desired separated representation of the stiffness matrix.
Depending on the number of parameters, n p , and the number of nodes chosen to discretize the parametric domains, n j , the separated expression of the parametric stiffness matrix might involve a large number of terms, n tot . It is possible to reduce the computational cost of the following calculations by employing the PGD-compression, available in the encapsulated PGD-toolbox [43]. The idea is to perform an L 2 projection of the expression of Eq. (32) to reduce the number of terms in the summation while maintaining an accurate representation of K(µ). In a similar fashion, a separated representation of the mass matrix can be also obtained. As it will be shown by means of numerical examples, the main advantage of the proposed algebraic technique is its flexibility which in general allows to add an arbitrary number of geometric parameters as variables of the problem. In addition, the nonintrusive character of the proposed ROM, makes the approach proposed in this work suitable for industrial applications.

Numerical examples
In this section two numerical examples are presented in order to show the properties of the proposed method. The first example is used to illustrate the numerical properties of the proposed PGD-IR method when both material and geometric parameters are considered. In the second example, the method is applied to a more realistic industrial case involving three parameters. Furthermore, a multiobjective optimization study is performed, which proves the potential of the PGD-IR method in the context of design optimization problems.

Parametric inertia relief with material and geometric parameters
A pure torsion test case is considered for an unconstrained linear elastic 3D structure characterized by one material and one geometric parameter, that are treated as additional coordinates of the problem. For a better readability, the two variables are denoted here with different symbols, that is µ ∈ M µ and θ ∈ M θ for the material and geometric parameters respectively.
As depicted in Fig. 4, the reference domainΩ consists of a block with dimensions where L x = 6, L y = 12 and L z = 1. The torsional load is given by two parallel forces of constant magnitude F = 10 acting on the positive and negative z direction respectively and applied at the points P = (2, 4, 1/2) and Q = (−2, 4, 1/2). Fig. 4 also shows the spatial discretization employed, consisting on a regular mesh with 236 nodes and 742 linear tetrahedral elements.
The physical domain Ω(θ) depends upon the geometric parameter and it is split into two non-overlapping subdomains Ω A (θ) and Ω B (θ). The parametric Young's modulus E is defined as where the Young modulus E A (µ) is considered varying in the range M µ = [10,410], and M µ is discretized with a uniform distribution of n µ = 41 points. The Poisson's ratio and the density are assumed constant in the whole domain and taken as ν = 0.3 and ρ = 1 respectively.
The geometrically parametrized domain Ω(θ) is described with the Cartesian coordinates x, and it is defined as the image of the reference domainΩ, with reference coordinateŝ x, via a geometric mapping Ψ(x, θ), namely The parameter θ is taken to be in the interval M θ = [0, 0.5], and M θ is discretized with a uniform distribution of n θ = 21 points.  The objective of this numerical test is to employ the proposed PGD-IR approach to obtain a computational vademecum able to describe the variation of the solution with respect to the material and geometric parameters.
Following the proposed PGD-IR framework, the first step consists of choosing a reference set of six degrees of freedom able to counteract the rigid body motions of the structure. Next, in order to employ the encapsulated PGD toolbox, it is necessary to define the input data (i.e. stiffness matrix, mass matrix, force vector) in a separated format. By using the linear dependence of the stiffness matrix on the Young's modulus, an analytical separable representation of the stiffness matrix with respect to µ can be easily constructed. For the geometric parameter θ, the algebraic PGD toolbox is employed, as discussed in detail in Sec. 4.3. For every nodal value of the geometric parameter θ p = [θ 1 , θ 2 , . . . , θ n θ ] T , the geometrically deformed mesh is generated according to the mapping of Eq. (34), and two stiffness-like matrices K A (θ p ) and K B (θ p ) are computed. The quantity K A (θ p ) is calculated by imposing the Young's modulus (E A , E B ) = (1, 0), thus accounting for the contribution of the finite elements belonging to the subdomain Ω A (θ p ) to the global stiffness matrix. Analogously, K B θ p corresponds to the choice (E A , E B ) = (0, 1) and accounts for the contribution of the finite elements belonging to the subdomain Ω B (θ p ). Once these matrices are sampled in the parametric nodes n θ , a separated form of the parametric global stiffness matrix is readily available, namely with k i (θ p ) = δ p,i , for every p = 1, 2, . . . , n θ . In this example, a PGD-compression was performed, which is always advisable when the number of PGD-terms is large, and an accurate approximation of the stiffness matrix was obtained in the known PGD format In this case, after performing compression with a tolerance tol = 10 −5 , the number of PGD terms was reduced to N K = 10. Following the same procedure, the PGD approximation of the parametric mass matrix is also obtained, namely Please note that the mass matrix is actually independent on the Young modulus, that is m i (µ) = 1. However, the general expression of Eq. (37) is used to maintain a consistent notation for all the inputs of the PGD-IR approach. Finally, the global forcing vector is also written in the general separated form where, again, it is worth emphasizing that F is the standard FE forcing vector and f (µ) = f (θ) = 1, because the right hand side is not dependent on the material parameter and, for the given set of forces applied to the structure is also independent on the geometric parameter.
The computation of the separated form of the stiffness and mass matrices and the forcing vector completes the pre-process required to apply the proposed PGD-IR approach. Next, the three steps of the PGD-IR approach can be sequentially completed. As detailed in Remark 1, the three steps involve a parametric problem because not only material parameters are considered but also geometric parameters, leading to a generalized solution that can be written as It is worthy to mention that the proposed PGD-IR approach was implemented in a Matlab routine which acts as a black-box, following the philosophy of the encapsulated PGD toolbox. In fact, the routine only requires to receive the input quantities in a separated form in order to return the output in the same separated form. Fig. 6 plots the evolution of the amplitude β U of each PGD mode. It can be observed that the amplitude rapidly decreases as the number of modes is increased. With 15 computed modes the amplitude of the last mode is almost four orders of magnitude lower than the amplitude of the first mode. In addition, the results show that the first four modes capture the most relevant information of the generalized solution as the fifth and subsequent modes have an amplitude at least two orders of magnitude lower than the amplitude of the first mode.
The first four normalized spatial modes are shown in Fig. 7, whereas the first four parametric modes are displayed in Fig. 8. The spatial modes provide an illustration of  the deformation induced by the four most relevant modes of the generalized solution. The parametric modes corresponding to the material illustrate that the four modes have the maximum contribution to the generalized solution for µ = 10. As the material property approaches the maximum value of µ = 410, the third and fourth mode have less influence on the solution. Finally, the modes corresponding to the geometric parameter have a more global character, proving the extra difficulty in solving geometrically parametrized problems. In order to get a particularized solution for a chosen set of the parameters (μ,θ), the correspondent function values u i µ (μ) and u i θ (θ) are evaluated for each PGD-mode i and then multiplied by the correspondent spatial mode and amplitude. Fig. 9 shows the particularized solutions in terms of deformed configuration and equivalent von Mises stress field for nine specific sets of parameters. The dominant character of the first spatial mode of Fig. 7 can be clearly observed, whereas the magnitude of the stress highly depends on the parametric choice. Please remember that these particularized solutions were obtained in real-time during an online post-process step.
In order to validate the PGD results, the accuracy with respect to the full-order FE computations is measured as the relative error between the PGD and FE solutions in the It is worth noting that to compute this error measure, the problem is solved by means of the standard FE method for each possible combination of the parameters, that is n µ ×n θ = 21 × 41 = 861 FE simulations. Fig. 10 shows the evolution of the relative error with respect to the number of PGD modes. As expected, the level of accuracy increases as the number of modes increases, up to a user-defined tolerance, which in this case was chosen equal to 10 −3 . Note that the PGD solution converges to the desired tolerance with only nine PGD modes. An interesting   advantage of the PGD method with respect to the standard FE method concerns the storage memory. In fact, the obtained PGD computational vademecum needs˜74 KB of storage memory versus the˜6650 KB needed to store all the 861 full-order FE solutions. Computational time is not particularly significant in the PGD context. In fact, the main goal is to provide a method which is able to explore an arbitrary large parametric space with only one offline computation. Nevertheless, an interesting comparison is shown in Table 1 where the number of iterations needed by the alternating direction scheme for the computations of each PGD mode is provided. As the cost of each iteration corresponds to the cost of a full-order FE simulation, the results in Table 1 show that the cost of the proposed PGD-IR is equivalent to 161 full-order solutions, compared to the 861 full-order computations required by the standard FE approach.
Finally, a major advantage of computing a PGD computational vademecum is the possibility to explore the design space and check, in real time, the effects of the design parameters on a predefined quantity of interest (QoI). As an example, the relative displacement in the z direction, ∆U PQ (z), of the points P and Q (see Fig. 4) is selected a QoI. The variation of the chosen QoI in the parametric space is depicted in Fig. 11.

Industrial application: dummy car test
The PGD-IR method is now employed to solve a more realistic problem, which is the static global torsional stiffness analysis of the BIW structure of a generic car. The geometry of the BIW is shown in Fig. 12. Two couples of parallel and opposite forces are applied at the front and rear shock towers, such that two opposite torsional moments of magnitude 1 Nm are generated. The FE model, is discretized with isoparametric quadrilateral shell elements. The material is linear elastic and it is characterized by a Young's modulus E = 207 GPa, Poisson's ratio ν = 0.29 and density ρ = 7.82 kg/m 3 . In this example, the thickness of three car components highlighted in Fig. 13, that usually play a role in the characterization of the global stiffness of the car, are introduced as extra coordinates of the problem. The three parameters are denoted by µ = [µ 1 , µ 2 , µ 3 ] T and they vary in the  intervals M j = [0.7, 1.5] mm, for j = 1, 2, 3. The three parametric domains are discretized with n 1 = n 2 = n 3 = 9 equidistant nodes. The goal of this example is to demonstrate potential of the proposed PGD-IR approach, able to produce a generalized solution that enables a designed to check how the overall static stiffness of the vehicle is affected by any change of the introduced parameters. This is done by computing the equivalent torsional stiffness (ETS), which is defined as a function of the front and back twisting rotations of the car body when a torsion load is applied (see Fig. 14), namely where the two angles α AB and α CD are defined as In the above expressions, U z (P ) denotes the displacement in the z direction at point P and L PQ denotes the distance between the points P and Q.
The proposed PGD-IR approach is employed following the same procedure described in the previous example. In the preprocess stage, the commercial FE package MSC-Nastran is now used to sample the parametric input matrices. A script was prepared to automatically produce a new Nastran input file (.bdf and .dat files) for each possible combination of the parameters. The generated files were then read by the Nastran solver, where the matrices were assembled (without solving the problem) and stored in a plain text format. Afterwards, the matrices were read by the developed Matlab routine to be expressed in the required separated form, namely with n tot = n 1 ×n 2 ×n 3 = 729, while K i and k i j (µ j ) being defined as described in Section 4.3. The mass matrix was obtained in a similar fashion and the separated force vector was computed. Finally, the separated expression of both the parametric stiffness and mass is compressed to minimize the number of terms in the separated form.
As already shown in the previous example, once the input data K PGD (µ), M PGD (µ) and F PGD (µ) are pre-processed and expressed in a separated form, the cascade scheme of the encapsulated PGD is used to sequentially solve the three steps involved in the developed PGD-IR approach, such that the final solution is obtained as The amplitude of the PGD modes is shown in Fig. 15. In this example, with a more complex geometry and a larger number of parameters, it can be observed that more modes are required to produce an accurate description of the multi-dimensional solution. With 21 modes the amplitude of the modes is approximately two orders of magnitude lower than the amplitude of the first mode.
The first four spatial modes are depicted in Fig. 16 amplified by a factor of˜1000 and Fig. 17 shows the first four normalized parametric functions.
To illustrate the full potential of the proposed PGD-IR approach, the application of the PGD-IR in a multi-objective optimization process is considered. The goal is to find the combination of parameters that maximize the ETS, while minimizing the mass of the three car components considered in this example. To this end, two objective functions are considered, namely where g 1 (µ) represents the mass of the material needed to manufacture the three car components, equal to the product of the material density ρ and the parametric volume. The latter is given by the sum of the products between the car components areas (A 1 , A 2 , A 3 ) times their variable thicknesses (µ 1 , µ 2 , µ 3 ). Clearly, this quantity is strictly related to the production cost. The objective function g 2 (µ) represents the parametric ETS defined in Eq. (41).
With the computed generalized solution U PGD (µ), an evaluation of the objective functions within a multi-objective optimization process only requires the particularization of the solution for a given set of the parameters. With the proposed PGD-IR approach, this evaluation can be performed in real time, making the overall cost of the optimization stage almost negligible. This is in contrast with a traditional approach where each evaluation of the objective functions require the assembly and solution of a new FE system of equations.
In this example, the optimization problem was solved by means of the gamultiobj function available in the Global optimization Toolbox released by Matlab. The function is able to find the Pareto front of multiple objective functions using a genetic algorithm. A Pareto front is a set of optimal points in the parametric space that represent a tradeoff between the objective functions. More specifically, a point is considered optimal if no objective can be improved without sacrificing at least one other objective.   18 shows the Pareto front in terms of the objective functions as well as the whole range of configurations that result from the parametric space M 1 × M 2 × M 3 . It is important to note that the optimization study allows to significantly reduce the range of solutions to be considered by a designer in the decision-making process. In fact, the variability in the two quantities of interest (mass and ETS), induced by the three parameters and calculated by PGD for all the possible combinations (PGD points in the left Fig. 18), is much larger than the number of points belonging to the Pareto front. Fig. 18 (right) plots the Pareto points in the parametric space (µ 1 , µ 2 , µ 3 ), where the correspondence to the Pareto front is described by same colors. In this example, the Pareto front was computed by assigning the same weight to the objective functions. Nevertheless, it is straightforward to obtain other fronts if the user wants to put more emphasis on one of the objective functions.

Conclusions
A nonintrusive algebraic PGD approach combined with the IR method for the solution of unconstrained problems being characterized by material and geometric parameters has been presented. The developed solver makes use of the Encapsulated PGD Toolbox developed by Díez et al. [43], which enables to perform algebraic operations for multidimensional data and allows to solve sequentially the three parametric problems required by the IR method.
An algebraic approach has been proposed to deal with geometric parameters by morphing a mesh generated for a reference configuration. The proposed method acts as a black-box, such that a nonintrusive interaction with commercial FE packages is possible.
Two numerical examples are used to underline the main properties of the method. The first example considers an academic test case with one material and one geometric parameter. The ability to compute a computational vademecum is shown and the accuracy of the generalized solution is measured by comparing the PGD solution to a set of stadard FE full-order solutions. It it shown that the proposed PGD-IR approach is able save almost the 99% of storage memory, requiring only the 20% of computational time needed by the FE method to solve the problem for every possible set of parameters. The second problem involves an industrial application for the static global stiffness analysis of a BIW structure of a car characterized by three parameters. This example shows the potential of the proposed PGD-IR approach and its ability to be integrated with a commercial FE package, such as MSC-Nastran. Finally, a multi-objective optimization was performed in order to show how the proposed approach can represent an important support to designers during the decision-making process. With the developed technique it is to possible to produce a computational vademecum displayed in a portable device to support the design engineers in the decision-making by evaluating in real time the impact of certain parameters on the global response of the structure.