Greedy maximin distance sampling based model order reduction of prestressed and parametrized abdominal aortic aneurysms

This work proposes a framework for projection-based model order reduction (MOR) of computational models aiming at a mechanical analysis of abdominal aortic aneurysms (AAAs). The underlying full-order model (FOM) is patient-specific, stationary and nonlinear. The quantities of interest are the von Mises stress and the von Mises strain field in the AAA wall, which result from loading the structure to the level of diastolic blood pressure at a fixed, imaged geometry (prestressing stage) and subsequent loading to the level of systolic blood pressure with associated deformation of the structure (deformation stage). Prestressing is performed with the modified updated Lagrangian formulation (MULF) approach. The proposed framework aims at a reduction of the computational cost in a many-query context resulting from model uncertainties in two material and one geometric parameter. We apply projection-based MOR to the MULF prestressing stage, which has not been presented to date. Additionally, we propose a reduced-order basis construction technique combining the concept of subspace angles and greedy maximin distance sampling. To further achieve computational speedup, the reduced-order model (ROM) is equipped with the energy-conserving mesh sampling and weighting hyper reduction method. Accuracy of the ROM is numerically tested in terms of the quantities of interest within given bounds of the parameter domain and performance of the proposed ROM in the many-query context is demonstrated by comparing ROM and FOM statistics built from Monte Carlo sampling for three different patient-specific AAAs.


Introduction
The potential of computational analysis to support clinical decision making is of great value for both physicians and patients. In particular the possibility to gain spatially and temporally resolved information on the patient-specific pathology at minimal intervention with the patient's body is driving this field of research. The human cardiovascular system is a specific example for the application of computational models [1] for risk assessment [2,3], planing of medical intervention and assessment of its effect [4][5][6] or general understanding of disease progression. More specifically, the pathology under consideration in this work is the abdominal aortic aneurysm (AAA).
An AAA corresponds to a dilatation of the aorta, which shows degraded mechanical properties in the widened segment [7] and is prone to rupture with highly-probable lethal outcome [8]. Given that aortic wall degradation and rupture is related to material failure, mechanical analysis of AAAs has been used for understanding and quantifying the development and progression of the disease [3,[9][10][11].
The AAA finite element models in this work are patient-specific, large-scale, stationary as well as materially and geometrically nonlinear. AAA geometries are extracted from medical screening images following the protocol in [12]. Given that imaged geometries are under blood pressure, an accurate computational model needs to impose a meaningful stress state, keeping the imaged configuration fixed. This is achieved in a modified updated Lagrangian formulation (MULF) [13,14] prestressing stage, wherein a physiological stress state is imprinted for an assumed diastolic blood pressure load. The vessel is subsequently deformed under further loading up to an assumed systolic blood pressure.
A common factor in most if not all works related to accurate state-of-the-art computational analysis of AAAs is a lack of knowledge on essential parameters related to mathematical modeling. This lack of knowledge results from the high inter-and intrapatient variety of AAA properties [15,16] and the limited accessibility to patient-specific data, given that the object of interest is located within the human body. From a computational perspective, this lack of knowledge typically results in the application of statistical methods, which attempt to propagate uncertainty through the computational model and involve sampling. Since the computational full-order models (FOMs) under consideration are nonlinear and large-scale, sampling with a high number of model evaluations quickly becomes too expensive to be practical in terms of computing power.
A well known approach to overcome the burden of impracticable requirements on computing power is projection-based model order reduction (MOR), which typically includes the following steps [17]. In a computationally expensive offline stage, the FOM is evaluated and a low-dimensional subspace is extracted from resulting solution snapshots in terms of the column span of an orthogonal matrix (the so-called reduced-order basis (ROB)). The ROB in turn is used to diminish the number of model degrees of freedom (DOFs) (also referred to as dimension or order in the current context). If constructed accurately, the reduced-order model (ROM) can replace the FOM in the given context of application.
The objective of the current work is to: 1. present a framework for the construction of a dimensionally reduced model (DROM) as well as a both dimensionally reduced and hyper reduced model (DHROM) for prestressed AAAs applying the Galerkin projection [18] for dimensional reduction and the energy-conserving mesh sampling and weighting (ECSW) method [19,20] for hyper reduction. The AAA models are parametrized in two material (low-strain range and high-strain range stiffness) and one geometric (AAA wall thickness) parameter. 2. demonstrate the applicability of both ROMs for assessment of the von Mises stress and the von Mises strain field in the aortic wall within bounds for the model parametrization. These quantities of interest are relevant in AAA rupture strati-fication and are therefore of essential importance for the progression of the disease [3,9,10]. 3. demonstrate the robustness of the presented framework by investigating three patient-specific computational examples which differ in geometry, parameter domain bounds as well as the number of DOFs.
Several techniques in the realm of the mechanical analysis of aneurysms have been proposed in the past to overcome the burden of limiting computing power. One example is the application of a computationally cheap intermediate mapping, which is utilized for sampling in place of the FOM. An example can be found in [3], wherein the authors fitted an inverse power-law function to represent the relation between aneurysm wall thickness and peak wall stress. In [21] a polynomial chaos expansion is built in order to investigate aneurysm wall stress assuming uncertainty in two material parameters, the wall thickness as well as the arterial pressure. A stochastic collocation method can be found in [22], wherein the authors interpolate the Navier-Stokes flow solution in order to evaluate the mean shear stress over the vessel wall.
Alternatively, the application of a cheap and possibly inaccurate model in terms of a multi-fidelity approach is presented in [10,23]. Therein, the cheap model is not supposed to replace the high-fidelity model, instead it rather serves as a means to decrease the number of high-fidelity model evaluations by providing additional information. A similar stochastic structure of the high-fidelity and the low-fidelity model is a prerequisite.
Also the applicability of projection-based MOR for computational feasibility of largescale aneurysm models has been demonstrated in the past. In [24], the authors address variable inflow angles and build a ROM for AAA hemodynamics. In [25], a ROM for the prediction of periodic regime hemodynamics of a cerebral aneurysm is derived.
We motivate the application of projection-based MOR for the following reasons. A surrogate model constructed by projection-based MOR will recover the FOM, if the degree of reduction is reversed. In this sense, projection-based MOR is consistent with FOM physics and contrasts the idea of an intermediate mapping as described above, given that such a mapping only exploits local FOM physics by sampling. The mentioned multi-fidelity approach incorporates the contribution of inaccurate information to specific quantities of interest. As opposed to projection-based MOR, a surrogate model producing highdimensional information and being able to serve as inexpensive FOM replacement is not created.
To the authors knowledge, no parametrized projection-based ROM has been presented for prestressed, large-scale, patient-specific and nonlinear solid mechanics AAA models to date. In particular, application of projection-based MOR to a MULF prestressing stage is a challenging task, which is investigated in this work. This involves a mathematical reformulation of our MULF prestressing stage, given that the original formulation accumulates an imprinted deformation gradient instead of computing displacement modes and therefore is not suitable for snapshot collection. Additionally, a sampling strategy combining greedy maximin distance sampling on parameter space subdomains and the concept of subspace angles is presented for snapshot collection and subsequent construction of the ROB.
The remainder of this paper is organized as follows. We first introduce the patientspecific AAA computational model. Special interest in view of projection-based MOR is devoted to the prestressing stage. Next, we present the DROM as well as the DHROM, assuming a given ROB and ECSW displacement modes and continue by describing the approach for construction of the ROB and ECSW displacement modes. For this purpose, a greedy maximin distance sampling approach and a stopping criterion based on subspace angles is applied. Finally, we present numerical experiments on three patient-specific AAA models and conclude.

Methods
The framework presented in this work includes a large-scale finite element model for AAA simulation, a projection-based MOR process and a sampling strategy for the construction of the ROM. These building blocks are described in the following.

Computational modeling of abdominal aortic aneurysms
In this section, we introduce the computational model in terms of its governing equations. We differentiate between the prestressing stage and the deformation stage, which, when combined, yield a mechanical state of the aortic segment under systolic blood pressure. Particular focus is placed on the prestressing stage, given that special treatment is required for the purpose of snapshot collection.

Patient-specific computational model
Our computational model consists of an aortic segment, which fully includes the AAA as well as short segments of the iliac arteries, see [12] for a detailed description of the workflow from imaging to finite element simulation. The aortic vessel is treated as an elastic solid consisting of an intraluminal thrombus (ILT) and the aortic wall. Pressure is exerted on the luminal (i.e. inner) surface of the ILT and the aneurysm is loaded to an assumed systolic blood pressure, which is the mechanical state of interest. The proximal and distal end surfaces of the model are constrained by a zero-displacement Dirichlet condition for vessel fixation. Figure 1 exhibits an example of a patient-specific computational domain.

Model equations
The governing equations read with The weak form of the governing equations is given by the principle of virtual work (PVW) δW, δW int and δW ext denote the total, internal and external virtual work, P denotes the first Piola-Kirchhoff stress tensor, N is the outward normal vector in the reference configuration and p,0 denotes the reference configuration pressure load surface (i.e. the luminal ILT surface). We emphasize that the traction boundary condition T depends on the displacement field u, see Eq. (4). Therein F (u) = I + ∂u ∂X is the deformation gradient with respect to the reference configuration, X ∈ 0 denotes reference configuration material coordinates, J (u) is the deformation gradient determinant and p is the pressure.
We make use of hyperelastic constitutive relations P = ∂ ∂F (6) introducing the strain-energy function and apply an isochoric-volumetric split for ILT as well as the vessel wall strain-energy wall (Ī 1 , J) = wall iso (Ī 1 ) + wall vol (J ), wherein are the first and second principal invariant of the modified right Cauchy Green tensor In more detail, we model the isochoric strain-energy contribution of the ILT as given in [12,26] ILT iso (Ī 1 , and the isochoric strain-energy contribution of the vessel wall as given in [9,12] wall iso (Ī 1 ) = α( The parameter c is a stiffness parameter of the ILT, while α (referred to as α-stiffness in the following) and β (referred to as β-stiffness in the following) can be interpreted as lowstrain range and high-strain range stiffness of the vessel wall, respectively. The volumetric parts wall vol , ILT vol of the strain-energies are chosen as given in [12,27] x vol (J ) = with x ∈ {ILT, wall} and κ wall , κ ILT being sufficiently large to reflect almost incompressible material behavior.

MULF prestressing
AAA geometries obtained from computed tomography imaging are exerted to blood pressure. From a continuum mechanics perspective, this corresponds to a non stressfree reference configuration [14,28,29]. Our simulations are therefore divided in two stages: The prestressing stage, which aims at imprinting a physiological stress-state into the imaged (i.e. fixed) geometric configuration at assumed diastolic blood pressure, is performed first. At second, the vessel is loaded to an assumed systolic blood pressure at evolving geometry in the deformation stage. We apply the Modified Updated Lagrangian Formulation (MULF) [14] prestressing approach in the first stage. MULF is an efficient prestressing method which especially was validated for the simulation of AAAs [10,12,13,30]. In the MULF prestressing approach an imprinted prestress deformation gradient F p is built up incrementally with boundary conditions evaluated at the imaged configuration.
Snapshot collection as required for data-driven construction of a ROB (cf. section "Construction of reduced-order model components") is not possible for MULF prestressing, given that displacement modes are not generated. To overcome this problem, we present a reformulation of MULF prestressing, shifting the wanted quantity from the prestress deformation gradient F p to a virtual prestress displacement field u p .
For consistency, we briefly review the original MULF prestressing formulation from a continuum mechanics perspective (details on implementation in the realm of the finite element method can be found in [14]) and state the mentioned reformulation in direct comparison with the original.
As starting point we recall the following kinematic relations. Given a virtual displacement fieldũ, from a virtual configuration˜ X to the current configuration x, a displacement field u from 0 X to , a deformation gradient F = ∂x ∂X and a virtual deformation gradientF = ∂X ∂X , we state As a result, the identical first Piola-Kirchhoff stress field P can be expressed as defining Applying the introduced notation into the PVW, we review the original MULF prestressing and subsequent deformation stage as In prestressing stage, find F p such that : In deformation stage, find u d (with given F p ) such that : Equation (22) implicitly defines the prestress deformation gradient F p , which is evaluated applying an assumed diastolic blood pressure load T (0, p dia ) at the known imaged geometry. Equation (23) utilizes the precomputed deformation gradient F p in order to evaluate the deformation stage displacement field u d applying an assumed systolic blood pressure load T (u d , p sys ) at the deformed geometry.
Recalling Eqs. (20) and (21), we can equivalently state the prestressing and deformation stage PVW as In prestressing stage, find u p such that : In deformation stage, find u d (with given u p ) such that : A comparison of (22), (23) with (24), (25) reveals the following. Instead of seeking a prestress deformation gradient F p , we solve for a virtual prestress displacement field u p fulfilling the PVW at a diastolic blood pressure load of the imaged geometry T (0, p dia ). u p is then used in the deformation stage to account for the stress in the imaged configuration at a systolic blood pressure load of the deformed configuration T (u d , p sys ).
We emphasize that the reformulation from (22), (23) to (24), (25) corresponds to a mathematical transformation of variables, physics remains unchanged. We also emphasize that both formulations are a well-posed approximation to the ill-posed inverse design problem as further detailed in [14]. From the perspective of projection-based MOR however, formulation (24), (25) enables a collection of virtual prestress displacement mode snapshots, an essential step in the data-driven construction of the ROB.

Finite element discretization
Applying the usual finite element discretization to the PVW for the MULF prestressing and deformation stage gives In deformation stage, find wherein u (e) = (e) d (e) , δu (e) = (e) δd (e) are the continuous element-wise displacement field and weighting function, which are interpolated by finite element shape functions contained in (e) and the element-wise displacement and weighting degree of freedom (DOF) vectors d (e) , δd (e) , respectively. Furthermore, we introduced the computational domain mesh element set E as well as the set F of elements loaded by the pressure load boundary condition. Given element-wise internal and external force vectors such that Eqs. (26) and (27) in assembled form read In prestressing stage, find d p such that : In deformation stage, find d d (with given d p ) such that : Thereby the global internal force vector Summarizing, we denote the high-fidelity finite element model residual as r : wherein the deformation stage only can be evaluated after the prestressing stage, which yields the virtual prestress displacement field d p as a solution. The nonlinear finite element system of equations in residual form reads and is solved applying Newton-Raphson iterations.

Reduction of the full-order model
In this section we briefly review the well known Galerkin projection, which yields a dimensionally reduced computational model. For nonlinear problems, the Galerkin projection is usually not sufficient to gain substantial computational speedup, given that the fullorder residual still needs to be assembled. For this reason, we additionally review the energy-conserving mesh sampling and weighting [19,20] hyper reduction method, which approximates the full-order residual with only a small subset of assembled mesh elements.

Galerkin projection on linear subspaces
The Galerkin projection has proven its applicability in structural mechanics problems [19,31,32]. Assuming a given orthogonal ROB V (its construction will be discussed in section "Construction of reduced-order model components") the dimensionally reduced model (DROM) retrieved from the Galerkin projection reads withd ∈ R n assuming n N . The argument of the residual is restricted to the column span of the ROB Vd ∈ span(V ), which corresponds to a reduction of the number of DOFs. Consistently, the number of equations is reduced by multiplication with the transposed ROB V T r(Vd) ∈ R n . As a result, application of the Newton-Raphson iteration scheme leads to wherein J r is the residual Jacobian with respect to the displacement field d. Equations (38), (39) reveal that only low-dimensional linear systems of equations have to be solved.

Hyper reduction of internal force contribution
The Galerkin Projection (37) leads to a dimensionally reduced model, however the FOM residual r(Vd) still needs to be evaluated together with its Jacobian J r (Vd) throughout Newton-Raphson iterations (38). Especially assembly of the internal force component of the residual (cf. Eq. (34)) is time consuming, given that every element of the computational mesh needs to be evaluated.
To reduce the cost of evaluation and assembly of the residual and its Jacobian, we apply the energy-conserving mesh sampling and weighting (ECSW) hyper reduction scheme [19,20] and give a brief review in the remainder of this section for completeness and adaption to the current context of application.
The idea is to replace the internal force vector f int ∈ R N with a surrogatef int ∈ R N , which will result in an accurate approximation after projection, that is f int is retrieved by a weighted assembly of a small mesh element subset and is derived from the requirement of an accurate approximation of the internal virtual work, which can be written as or for a dimensionally reduced model d, δd ∈ span(V ) Applying a sum over all element internal force contributions we can rewrite using L (e)T to extract DOFs of element e from a vector with global DOF numbering into a smaller vector with element DOF numbering. We now seek for an approximationW int (d, δd) of the internal virtual work such that with In contrast to (43), (45) only contains a summation over a reduced element setẼ. Additional non-negative element weights w (e) ∈ R + are introduced for approximation (44) to become feasible with a small cardinality of the reduced element set. A remaining question is the actual choice of elements inẼ as well as the element weights w (e) . For this reason, (44) is turned into an optimization problem by restriction to a finite set of displacement modesŜ being a set of known displacement modes (referred to as ECSW displacement modes here, the actual selection of modes will be discussed in section "Construction of reduced-order model components"). Note that S consists of m (even number) modes corresponding to virtual prestress displacement modes d p(i) and the sum of virtual prestress displacement and deformation stage displacement modes In its unassembled shape, the restriction of Eq. (44) to ∀d ∈Ŝ reads In order to keep the cardinality of the reduced element setẼ low, approximation (48) has to be accurate with a minimal number of non-zero weights. The corresponding optimization problem is The zero-norm (•) 0 counts the number of non-zero entries and is used as the objective function applied to the vector of element weights w, which is constrained to have non-negative values expressed by its minimum entry min(w) being non-negative. This constraint is required in order to ensure a positive semi-definite Jacobian of the internal force vector [19]. The other constraint is a fulfillment of Eq. (48) up to the relative tolerance ε h . Consequently, with vector valued entries and with vector valued entries wherein i ∈ {0, . . . , m − 1}, j ∈ {0, . . . , |E| − 1} andd i are vectors from the setŜ. Optimization problem (49) can be approximately solved with a sparse non-negative least-squares solver. Thereby sparse refers to the solution vector w, in the sense that the number of non-zero entries is kept minimal. For details on the iterative solver the reader is referred to [19].
An (approximate) solution to (49) returns the element weights as well as the reduced element set by an extraction of elements with non-zero weights. As a consequence, the hyper reduced internal force vector and its Jacobian read with the element stiffness J ∂d (e) . We denote the dimensionally reduced as well as hyper reduced model (DHROM) as V Tr (Vd) = 0 (56) and state the corresponding Newton-Raphson iterations whereinr is a residual approximation usingf int andJ r is the corresponding hyper reduced residual Jacobian.

Construction of reduced-order model components
In the given many-query context, the residual (34) depends on a modifiable set of model parameters. Introducing a parameter vector with [lb i ; ub i ] μ i being the lower and upper bounds for parameter μ i , we extend the notation of Eq. (35) to and attempt to find a ROB V and a set of ECSW modes S such that the resulting DROM as well as DHROM will accurately approximate the FOM solution for all μ ∈ P.
A prerequisite for an accurate ROM is that the FOM solution d(μ) can be accurately represented within the column span of the ROB This motivates a data-driven approach for construction of ROBs by a collection of FOM solution snapshots at different parametric configurations and subsequent orthogonalization with or without data compression [33].
In case of the presented stationary AAA computational model two snapshots per parametric configuration (virtual prestress displacement and deformation stage displacement) are retrieved and organized in the so-called snapshot matrix with S ∈ R N ×m . Two data-driven approaches for the construction of the ROB have gained special interest in projection-based MOR. Proper orthogonal decomposition (POD) can be used to orthogonalize S yielding a ROB V pod ∈ R N ×m such that [17] wherein V n pod corresponds to a selection of the first n columns of V pod with n ≤ m and (•) F is the Frobenius norm. Consequently, POD is used whenever solution snapshots can be accurately represented by a low-dimensional subspace.
The second approach for construction of the ROB are greedy methods [34]. The idea herein is to successively build the ROB by evaluating selected configurations within the parameter domain and enrich the span of the ROB by the span of the newly computed snapshots. Based on an (hopefully inexpensive and sharp) a posteriori error estimator, greedy methods attempt to find solution snapshots which are represented worst by the ROB constructed up to this point. However, even this local optimization problem quickly becomes too expensive. Several approaches [35][36][37] have been presented to date to overcome this computational bottleneck.
In order to avoid evaluating a posteriori error estimates in every greedy iteration, we apply a selection of solution snapshots from a greedy maximin distance sampling together with a stopping criterion based on subspace angles and an exclusion of subdomains. We dedicate section "Maximin distance design" and "Subspace angles for the comparison of subspaces" to the notion of maximin distance design and subspace angles, respectively. Section "A greedy maximin distance sampling approach for the construction of solution subspaces" introduces the actual sampling algorithm.

Maximin distance design
Space-filling designs is a topic from design of experiments. Maximin distance (MMD) is introduced in [38] as a criterion which can be used to rate the space-filling property of a design or to construct space-filling designs by optimization of that criterion. A greedy version with reduced computational complexity is presented in [39] under the name "Coffee-House Design" and a recent review on maximin distance sampling can be found in [40]. In contrast to a globally optimal MMD design, the greedy MMD design can be evaluated iteratively.
We apply the following terminology. A point is a specific instance of the parameter vector μ, also referred to as parametric configuration. Points are distributed by the sampling algorithm in the parameter domain. A sample corresponds to FOM solution snapshots at a given point.
Algorithm 1 depicts the steps for the selection of a greedy MMD point. Given an input grid i ⊂ P as subset of the parameter space and a set of previously chosen points c ⊂ P, the next point μ ∈ i is chosen such that the minimal distance to a neighboring point p ∈ c is maximized in a reference hypercube. Thereby χ and χ −1 map from physical domain to reference hypercube and vice versa, respectively.

Algorithm 1 MaxiMinPoint( i , c ) (select a greedy MMD point)
Input: input grid i ⊂ P, previously chosen points c ⊂ P Output: selected grid point μ transform grids to reference hypercube 2: μ = arg max q∈˜ i minp ∈˜ c q −p 2 get next point in reference hypercube 3: return χ −1 (μ) return point in physical domain The steps for a greedy MMD design on a training grid t ⊂ P are depicted in Algorithm 2. Figure 2 illustrates a greedy MMD design, wherein the first parametric configuration was chosen at random.
The idea of MMD sampling for the purpose of surrogate modeling in general is discussed broadly in literature [41][42][43][44]. Specific applications for instance can be found in [45], wherein the authors use the notion of MMD in their algorithm to sample cut lines and planes of the parameter domain in order to construct a radial-basis-function approximation surrogate. In [46], MMD sampling is used to distribute points in Voronoi cells for multifidelity radial-basis-function metamodeling.

Subspace angles for the comparison of subspaces
Subspace angles (or principal angles) are a concept from matrix computations [47]. Given two matrices Y ∈ R N ×n , Z ∈ R N ×m with n ≤ m, subspace angles can be defined recursively as the minimum value while the corresponding principal vectors follow from the minimization arguments The sets depend on the minimization arguments y j and z j of the previous k iterations. A maximum subspace angle of 0 • indicates that span(Y ) ⊆ span(Z), while a maximum subspace angle of 90 • indicates that there is at least one direction in span(Y ) which is orthogonal to span(Z). More general, the maximum subspace angle can be interpreted as a distance measure from span(Y ) to span(Z). In the following we will refer to the maximum subspace angle as the subspace angle distance (SAD). Figure 3 illustrates subspace angles for 2D subspaces embedded in a 3D space. Algorithm 3 [47] states the computation of subspace angles applying a singular value decomposition.

Algorithm 3 SSA(Y , Z) (computation of subspace angles)
Input: Y ∈ R N ×n , Z ∈ R N ×m with n ≤ m Output: subspace angles α In projection-based MOR, subspace angles have been used for the purpose of interpolation and sampling. In [48][49][50], the authors present and apply a subspace angle interpolation of ROBs for flow problems. In [51], subspace angle interpolation is performed with respect to the diffusion coefficient for a Diffusion-Convection-Reaction problem. The application of subspace angles as a stopping criterion for sampling has been presented in [52][53][54].

A greedy maximin distance sampling approach for the construction of solution subspaces
We expand the greedy MMD design with a stopping criterion based on the SAD and introduce adaptivity to the sampling by a division of the parameter domain into subdomains. Those subdomains are subsequently excluded from sampling and the algorithm stops, when all subdomains have been excluded. = ( sd,1 , . . . , sd,n sd −1 , sd,0 ) define subdomain tuple for iteration 6: while True do 7: for sd,i ∈ do iterate over subdomains 8: s(μ)) 11: α = max(α) 12: if α < α m then in case of small SAD We use the ROB for dimensional reduction by the Galerkin projection (37), while the snapshot matrix is used to compute the set S (47). Consequently, the number of ECSW displacement modes |S| coincides with the dimension of the subspace span(V ).
Note that the ROB enrichment strategy of Algorithm 4 (line 15:,16:) extends the span of the ROB by the span of the current snapshot matrix s(μ) in every greedy maximin iteration. As a consequence, we can state that span(V i ) ⊆ span(V j ) for j ≥ i, wherein V i and V j denotes the ROB after greedy maximin iteration i and j, respectively.
Following the terminology in [43], we classify the proposed approach as global, sequential / adaptive and fine-grained. Additionally, the MMD criterion as well as initial iterations over subdomains introduce the property of domain exploration, while subsequent sampling of a subset of subdomains amounts to local exploitation.
In more detail, by introducing subdomains the sampling algorithm can evaluate more samples in specific parameter domain regions as compared to others, if this need is identified by the subspace angle criterion. As a result, the purely explorative greedy maximin sampling receives a feedback from the parameter domain and sampling is refined. Refer to section "Application of greedy maximin distance sampling" for a demonstration in the given context of prestressed and parametrized AAAs.
A stopping criterion for sampling based on subspace angles already has been presented in [52][53][54]. In [52,53] the authors present adaptive sampling of a linear time-invariant state-space system, while in [54] adaptive selection of linearization points in trajectory piecewise linear approximation is in the focus of interest. In contrast to [52][53][54] the approach in this work aims at the construction of a global ROB instead of interpolation between parametric configurations. Additionally, the notion of MMD yields to finely granular sampling applicable to parameter domains with multiple dimensions.

Results and discussion
The proposed framework is applied to three patient-specific computational examples of AAAs. The ROB is constructed by greedy subdomain MMD sampling with 8 subdomains following Algorithm 4. The accuracy of the resulting DROMs and DHROMs is evaluated in terms of the quantities of interest (von Mises stress field and von Mises strain field in the aortic vessel wall) and wall clock timings are reported. The choice of von Mises type quantities of interest is based on preceding numerical studies on AAAs with emphasis on solid mechanics and rupture risk (e.g. as presented in [9,55,56]). Nonetheless, other quantities of interest could have been selected, given that AAA pathological progression is still subject to research. Finally, accuracy of the proposed MOR framework in a statistical sense is demonstrated by comparing maximum von Mises stress and von Mises strain probability distributions gained from FOM and DHROM sampling. The ILT is discretized using linear tetrahedral and pyramid elements, wherein pyramids are introduced to connect the ILT to the aortic wall, which is discretized using linear hexahedral elements with F-bar element technology [57]. Table 1 depicts information on the model discretization.

Patient-specific computational models
Referring to section "Computational modeling of abdominal aortic aneurysms", we quantify boundary conditions by a diastolic blood pressure of p dia = 87 mmHg (11.6kPa)    and a systolic blood pressure of p sys = 121 mmHg (16.1 kPa). The ILT stiffness c is interpolated linearly from a luminal stiffness of c = 2.62 kPa to a medial stiffness of c = 1.98 kPa and from the medial stiffness to an abluminal stiffness of c = 1.73 kPa [26]. Together with the aortic wall thickness t the model parametrization is given as with the subscripts l and u denoting the lower and upper bound. Table 2 exhibits parameter domain lower and upper bounds for the three models. The bounds were computed from patient-specific Log-normal probability distributions from [58] for each entry of the parameter vector μ. In more detail, the parameter domain bounds are chosen as (γ l , γ u ) = (Q log (0.025; μ γ , σ γ ), Q log (0.975; μ γ , σ γ )) for γ ∈ {α, β, t} with being the p-percentile value for a Log-normal distribution with expectation μ γ and standard deviation σ γ . erf denotes the error function. Consequently, the range within the chosen parameter domain bounds covers 95% of realizations of μ.
We perform 15 equally spaced load steps for the prestressing stage and 10 equally spaced load steps for the deformation stage. Multiple thousands of simulations were performed and postprocessed for the results presented in the following sections. Individual unconverged simulations were dropped from analysis.
For linear systems of equations arising in FOM simulations, we use an iterative, parallel GMRES solver with algebraic multigrid preconditioning implemented in Trilinos [59]. For the ROM linear systems of equations we apply a direct solver [60], given that arising linear systems have less than 100 unknowns.

Application of greedy maximin distance sampling
In our first numerical experiment, we create a one-shot (i. e. no adaptation) design distributing 200 points in the parameter domain (Algorithm 4 with sd = { sd,0 }, α m = 0.0 and stopping at | c | = 200). If a simulation fails to converge, a neighboring point is taken in the set of selected points c instead and the FOM is recomputed. The parameter domain grid is created from all combinations of 100 equidistantly placed points in each direction of the parameter domain axes. The initial point is chosen as the "minimum-value" point μ = [α l , β l , t l ] T (see Table 2) for each patient-specific example.
As a result, the greedy MMD design returns identical points (except for few individually shifted points due to convergence failure) in the reference cube for all three computational examples. The corresponding MMD in the reference cube is depicted in Fig. 7d. Simultaneously, the SAD (Algorithm 4 (line 11:)) is depicted in Fig. 7a-c, wherein we highlight SADs corresponding to (− − −)-octant configurations of the parameter domain (i.e. α < α m , β < β m , t < t m ) in blue and SADs corresponding to (+ + +)-octant configurations of the parameter domain (i.e. α > α m , β > β m , t > t m ) in dark red (given α m = 1 2 (α l + α u ), β m = 1 2 (β l + β u ) and t m = 1 2 (t l + t u ) as the axes mid values). The resulting distribution of SADs is affected by two contributions. The first contribution is the MMD for a newly set point. As one can observe from Fig. 7d, the MMD strongly decreases initially, while a stagnation occurs with an increasing number of samples. This behavior also reflects in the SAD, which shows a pronounced decay in the beginning and increased scattering with ongoing stagnation of the MMD. The second contribution are different sensitivities of FOM snapshots with respect to the parameter domain. For instance, (+ + +)-octant value parametrizations (dark red points) yield lower subspace angles than (− − −)-octant value parametrizations (blue points), such that the (+ + +)octant of the parameter domain can be said to show lower sensitivity in solution snapshots. For practical reasons (see numerical examples presented next), we are only interested in the region of pronounced decay of the SAD. As a consequence, distance in the reference parameter space is a suitable and efficient sampling criterion for the problem at hand.
We include adaptivity to the greedy MMD design by introducing subdomains as presented in section "A greedy maximin distance sampling approach for the construction of solution subspaces". In more detail, we create eight equally shaped subdomains by splitting each parameter domain axis in 2 intervals and run Algorithm 4 with α m = 0.1, sd = { sd,0 . . . , sd,7 } and the initial configuration μ = [α l , β l , t l ] T . Figure 8 depicts the decay of SADs for each patient-specific computational model. As one can observe, sampling runs until the last sample in each subdomain yields a SAD below α m .  Table 3 depicts the number of distributed points in the individual subdomains. Note the differences in the point distributions, especially prominent for patient 3. Figure 9 depicts the selected parametric configurations in the 3D parameter domain. As one can see, the most sensitive subdomain 0 (compare Table 3) corresponds to the low-stiffness and thin vessel wall range of the parameter domain. This is plausible from a physical perspective, given that soft and thin-walled tissue will deform more than stiff and thick-walled tissue.

Patient-specific reduced-order models
We compute ROBs from greedy subdomain MMD sampling with parametric configurations as depicted in Fig. 9. A Galerkin projection (37) yields the patient-specific DROMs. Hyper reduction is achieved via ECSW (cf. section "Hyper reduction of internal force contribution") parallelized on 4 processors with global tolerance set to ε h = 10 −4 , see Eq. (49). The parallelization is implemented in terms of a domain decomposition as described in [20].     Figure 10 illustrates selected mesh elements, while Table 4 shows the number of DOFs and selected mesh elements. Inspecting Fig. 10, note the increasingly accurate sampling in the neighborhood of vessel fixation (proximal and distal vessel ending) and in regions with increased curvature of the vessel wall (compare with Figs. 5 and 6).

Accuracy of the reduced-order model
We evaluate accuracy in terms of the relative error in the quantities of interest (von Mises stress field σ vM and von Mises strain field e vM in the aortic wall). The relative error is given as wherein x ∈ {σ vM , e vM } corresponds to FOM quantities andx ∈ {σ vM ,ẽ vM } corresponds to ROM approximations. A validation grid with 1000 points in the parameter domain is used. The grid results from all combinations of 10 equidistantly placed points in each direction of the parameter domain axes. Note that the resulting grid corresponds to a full factorial design being created independently of the points used for the construction of the ROB by greedy subdomain MMD sampling. Figures 11 and 12 depict the corresponding errors. The majority (> 98%) of relative errors are below 1%, while individual runs show a relative error above 1%. We conclude, that DROM as well as DHROM are accurate models for the von Mises stress and von Mises strain field in the aortic wall in a statistical sense. Individual simulations might show increased relative errors, caution is required when applying the ROMs for the prediction of point estimates. We investigate the influence of α m and ε h , i.e. the maximum subspace angle defining the stopping criterion in Algorithm 4 and the relative tolerance for the ECSW algorithm introduced in Eq. (49). On the introduced validation grid with 1000 points, we evaluate as the mean value of all relative errors given the number of performed simulations n sim , x ∈ {σ vM , e vM } denoting von Mises stress or von Mises strain as the quantity of interest and RE i (x, x) being the relative error (71) of sample i.  Fig. 13b corresponds to DHROM results with ECSW meshes at different tolerances and an unchanged ROB at α m = 0.1. As can be observed from Fig. 13a, lower values for α m yield larger ROBs, while at the same time the mean relative error decreases indicating a more accurate DROM as expected. Figure 13b illustrates that stricter ECSW tolerances ε h lead to an increased number of selected mesh elements with decreasing mean relative error indicating a more accurate DHROM. We conclude, that both α m and ε h strongly influence ROM accuracy.

Monte Carlo sampling on the reduced-order model
To demonstrate applicability of the constructed ROMs for approximation of probability distributions in the quantities of interest, we compare the 99.9 percentile aortic wall von Mises stress (referred to as maximum von Mises stress in the following) and the 99.9 percentile aortic wall von Mises strain (referred to as maximum von Mises strain in the following) probability distributions retrieved from Monte Carlo sampling of the FOM and Monte Carlo sampling of the DHROM. Both models are evaluated on 10000 identical (per patient) parametric configurations drawn from the corresponding patient-specific Log-normal probability distributions. Table 5 depicts mean and standard deviation of the quantities of interest. As one can see, FOM and DHROM results are very close, relative errors are < 1%. Figure 14 depicts kernel-density-estimated probability distributions gained from FOM and DHROM samples. We apply the Gaussian kernel-density-estimator scipy.stats.gaussian_kde available in the SciPy [61] (version 1.3.0) ecosystem of the Python programming language. The plots show negligible differences between probability distributions gained from FOM and DHROM sampling.

Timing
We report wall clock timings of the patient-specific computational models as well as corresponding speedups in Table 6. All simulations in this section were performed on a workstation with Intel Xeon W-2133 (3.60GHz) processors.
The values in Table 6 are mean values corresponding to seven simulations (per patient) evaluated at face mid-points as well as the mid-point of the patient-specific parametric domains, compare with Table 2 for domain lower and upper bounds.
As the reader can observe, only slight speedup can be achieved with DROM models, given that the full-order residual as well as its Jacobian need to be evaluated. A rather Table 6 Timing for patient-specific models. All computations were performed on 4 cores (Intel Xeon W-2133 (3.60GHz)). The reported timings are mean values of seven simulations (per patient). Speedup is computed using FOM timing as a reference  Table 7 Timing for offline stage steps of patient-specific models performed on 4 cores (Intel Xeon W-2133 (3.60GHz)). Timing is given as multiple of a single FOM evaluation time, which in turn is estimated from the mean of seven simulations per patient substantial speedup can be achieved by DHROM models, recalling that only a small portion of the computational mesh is evaluated and assembled. Table 7 depicts offline stage timings, subdivided into the ROB construction stage by greedy subdomain MMD sampling (as described in section "Application of greedy maximin distance sampling") and the hyper-reduction by ECSW (as described in "Patientspecific reduced-order models"). For details on theory please refer to section "A greedy maximin distance sampling approach for the construction of solution subspaces" and section "Hyper reduction of internal force contribution".

Conclusions
We presented a framework for projection-based MOR of patient-specific AAA models. A dimensionally reduced model was built by a Galerkin projection on a low-dimensional subspace and a dimensionally reduced as well as hyper reduced model was built by the Galerkin projection and energy-conserving mesh sampling and weighting. Specific attention was dedicated to the MULF prestressing stage, given that the original MULF aims at a calculation of an imprinted deformation gradient and therefore is not suited for snapshot collection. A sampling algorithm relying on the maximin distance criterion was presented for the construction of low-dimensional solution subspaces. Therein, a stopping criterion based on subspace angles and an exclusion of subdomains was applied. Finally, three patient-specific computational examples with different complexities were demonstrated.
The proposed sampling algorithm led to comparable results in terms of reduced-order basis size as well as number of sampled mesh elements for all three patient-specific computational models. Subsequent experiments on ROM accuracy revealed relative von Mises stress and von Mises strain field errors below 1% for more than 98% of all simulations on a validation grid. We conclude that the proposed MOR framework is robust across patient-specific AAA geometries and parameter domains.
Direct Monte Carlo sampling on the dimensionally reduced as well as hyper reduced model was performed calculating the maximum von Mises stress and maximum von Mises strain in the vessel wall. Comparison with the corresponding FOM reference solution revealed a very good match between FOM and ROM kernel-density-estimated probability distributions for the maximum von Mises stress and the maximum von Mises strain in the aortic wall.
The proposed sampling algorithm led to appropriate dimensionally reduced (and hyper reduced) models as a conclusion from numerical experiments. However, a certification in terms of an upper bound for the error in the quantities of interest was not presented in this work. Furthermore, motivated by practical reasoning, this paper investigated numerical experiments on parametrizations in terms of two material and one geometric parameter. Deviations from this setup need further validation. We did not include calcification and did not distinguish between healthy and aneurysmatic sections of the simulated vessel in terms of material behavior. These are two example features that would yield a more realistic full-order model.