 Research article
 Open Access
 Published:
Greedy maximin distance sampling based model order reduction of prestressed and parametrized abdominal aortic aneurysms
Advanced Modeling and Simulation in Engineering Sciences volume 8, Article number: 18 (2021)
Abstract
This work proposes a framework for projectionbased model order reduction (MOR) of computational models aiming at a mechanical analysis of abdominal aortic aneurysms (AAAs). The underlying fullorder model (FOM) is patientspecific, 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 manyquery context resulting from model uncertainties in two material and one geometric parameter. We apply projectionbased MOR to the MULF prestressing stage, which has not been presented to date. Additionally, we propose a reducedorder basis construction technique combining the concept of subspace angles and greedy maximin distance sampling. To further achieve computational speedup, the reducedorder model (ROM) is equipped with the energyconserving 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 manyquery context is demonstrated by comparing ROM and FOM statistics built from Monte Carlo sampling for three different patientspecific 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 patientspecific 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 highlyprobable 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 patientspecific, largescale, 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 stateoftheart 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 patientspecific 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 fullorder models (FOMs) under consideration are nonlinear and largescale, 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 projectionbased model order reduction (MOR), which typically includes the following steps [17]. In a computationally expensive offline stage, the FOM is evaluated and a lowdimensional subspace is extracted from resulting solution snapshots in terms of the column span of an orthogonal matrix (the socalled reducedorder 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 reducedorder 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 energyconserving mesh sampling and weighting (ECSW) method [19, 20] for hyper reduction. The AAA models are parametrized in two material (lowstrain range and highstrain 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 stratification 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 patientspecific 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 powerlaw 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 NavierStokes 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 multifidelity approach is presented in [10, 23]. Therein, the cheap model is not supposed to replace the highfidelity model, instead it rather serves as a means to decrease the number of highfidelity model evaluations by providing additional information. A similar stochastic structure of the highfidelity and the lowfidelity model is a prerequisite.
Also the applicability of projectionbased 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 projectionbased MOR for the following reasons. A surrogate model constructed by projectionbased MOR will recover the FOM, if the degree of reduction is reversed. In this sense, projectionbased 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 multifidelity approach incorporates the contribution of inaccurate information to specific quantities of interest. As opposed to projectionbased 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 projectionbased ROM has been presented for prestressed, largescale, patientspecific and nonlinear solid mechanics AAA models to date. In particular, application of projectionbased 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 projectionbased 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 patientspecific AAA models and conclude.
Methods
The framework presented in this work includes a largescale finite element model for AAA simulation, a projectionbased 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.
Patientspecific 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 zerodisplacement Dirichlet condition for vessel fixation. Figure 1 exhibits an example of a patientspecific 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)
\(\delta W, \delta W_{\mathrm {int}}\) and \(\delta W_{\mathrm {ext}}\) denote the total, internal and external virtual work, \(\varvec{P}\) denotes the first PiolaKirchhoff stress tensor, \(\varvec{N}\) is the outward normal vector in the reference configuration and \(\Gamma _{p,0}\) denotes the reference configuration pressure load surface (i.e. the luminal ILT surface). We emphasize that the traction boundary condition \({\varvec{T}}\) depends on the displacement field \({\varvec{u}}\), see Eq. (4). Therein \(\varvec{F}({\varvec{u}}) = \varvec{I} + \frac{\partial {\varvec{u}}}{\partial {\varvec{X}}}\) is the deformation gradient with respect to the reference configuration, \({\varvec{X}} \in \Omega _0\) denotes reference configuration material coordinates, \(J({\varvec{u}})\) is the deformation gradient determinant and p is the pressure.
We make use of hyperelastic constitutive relations
introducing the strainenergy function \(\Psi \) and apply an isochoricvolumetric split for ILT as well as the vessel wall strainenergy
wherein
are the first and second principal invariant of the modified right Cauchy Green tensor \({\bar{\varvec{C}}} = {\varvec{F}}^T_{\mathrm {iso}} {\varvec{F}}_{\mathrm {iso}}\) with \({\varvec{F}}_{\mathrm {iso}} = J^{\frac{1}{3}} \varvec{F}\). In more detail, we model the isochoric strainenergy contribution of the ILT as given in [12, 26]
and the isochoric strainenergy contribution of the vessel wall as given in [9, 12]
The parameter c is a stiffness parameter of the ILT, while \(\alpha \) (referred to as \(\alpha \)stiffness in the following) and \(\beta \) (referred to as \(\beta \)stiffness in the following) can be interpreted as lowstrain range and highstrain range stiffness of the vessel wall, respectively. The volumetric parts \(\Psi _{\mathrm {vol}}^{\mathrm {wall}}, \Psi _{\mathrm {vol}}^{\mathrm {ILT}}\) of the strainenergies are chosen as given in [12, 27]
with \(\mathrm {x} \in \{\mathrm {ILT}, \mathrm {wall}\}\) and \(\kappa ^{\mathrm {wall}}\), \(\kappa ^{\mathrm {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 stressstate 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 \(\varvec{F}_p\) is built up incrementally with boundary conditions evaluated at the imaged configuration.
Snapshot collection as required for datadriven construction of a ROB (cf. section “Construction of reducedorder 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 \(\varvec{F}_p\) to a virtual prestress displacement field \(\varvec{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 \({\tilde{\varvec{u}}}\), from a virtual configuration \({\tilde{\Omega }} \ni {\tilde{\varvec{X}}}\) to the current configuration \(\Omega \ni \varvec{x}\), a displacement field \(\varvec{u}\) from \(\Omega _0 \ni \varvec{X}\) to \(\Omega \), a deformation gradient \(\varvec{F} = \frac{\partial {\varvec{x}}}{\partial {\varvec{X}}}\) and a virtual deformation gradient \({\tilde{\varvec{F}}} = \frac{\partial {\tilde{\varvec{X}}}}{\partial {\varvec{X}}}\), we state
As a result, the identical first PiolaKirchhoff stress field \(\varvec{P}\) can be expressed as
defining
Applying the introduced notation into the PVW, we review the original MULF prestressing and subsequent deformation stage as
Equation (22) implicitly defines the prestress deformation gradient \(\varvec{F}_p\), which is evaluated applying an assumed diastolic blood pressure load \(\varvec{T}({\varvec{0}}, p_{\mathrm {dia}})\) at the known imaged geometry. Equation (23) utilizes the precomputed deformation gradient \(\varvec{F}_p\) in order to evaluate the deformation stage displacement field \(\varvec{u}_d\) applying an assumed systolic blood pressure load \(\varvec{T}({\varvec{u}}_d, p_{\mathrm {sys}})\) at the deformed geometry.
Recalling Eqs. (20) and (21), we can equivalently state the prestressing and deformation stage PVW as
A comparison of (22), (23) with (24), (25) reveals the following. Instead of seeking a prestress deformation gradient \(\varvec{F}_p\), we solve for a virtual prestress displacement field \(\varvec{u}_p\) fulfilling the PVW at a diastolic blood pressure load of the imaged geometry \(\varvec{T}({\varvec{0}}, p_{\mathrm {dia}})\). \(\varvec{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 \(\varvec{T}(\varvec{u}_d, p_{\mathrm {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 wellposed approximation to the illposed inverse design problem as further detailed in [14]. From the perspective of projectionbased MOR however, formulation (24), (25) enables a collection of virtual prestress displacement mode snapshots, an essential step in the datadriven construction of the ROB.
Finite element discretization
Applying the usual finite element discretization to the PVW for the MULF prestressing and deformation stage gives
wherein \(\varvec{u}^{(e)} = \varvec{\Phi }^{(e)} \varvec{d}^{(e)}, \delta \varvec{u}^{(e)} = \varvec{\Phi }^{(e)} \delta \varvec{d}^{(e)}\) are the continuous elementwise displacement field and weighting function, which are interpolated by finite element shape functions contained in \({\varvec{\Phi }}^{(e)}\) and the elementwise displacement and weighting degree of freedom (DOF) vectors \(\varvec{d}^{(e)}, \delta \varvec{d}^{(e)}\), respectively. Furthermore, we introduced the computational domain mesh element set \(\mathcal {E}\) as well as the set \(\mathcal {F}\) of elements loaded by the pressure load boundary condition.
Given elementwise internal and external force vectors such that
Eqs. (26) and (27) in assembled form read
Thereby the global internal force vector \(\varvec{f}_{\mathrm {int}} = \sum _{e \in \mathcal {E}} \varvec{L}^{(e)} \varvec{f}_{\mathrm {int}}^{(e)} \), global external force vector \(\varvec{f}_{\mathrm {ext}} = \sum _{e \in \mathcal {F}} \varvec{L}^{(e)} \varvec{f}_{\mathrm {ext}}^{(e)} \), global displacement DOF vector \(\varvec{d}= \sum _{e \in \mathcal {E}} \varvec{L}^{(e)} \varvec{d}^{(e)}\) as well as global weighting DOF vector \(\delta \varvec{d}= \sum _{e \in \mathcal {E}} \varvec{L}^{(e)} \delta \varvec{d}^{(e)}\) result from an assembly of the corresponding elementwise vectors, while \(\varvec{L}^{(e)}\) is the usual finite element assembly operator towards the global system.
Summarizing, we denote the highfidelity finite element model residual as
wherein the deformation stage only can be evaluated after the prestressing stage, which yields the virtual prestress displacement field \(\varvec{d}_p\) as a solution. The nonlinear finite element system of equations in residual form reads
and is solved applying NewtonRaphson iterations.
Reduction of the fullorder 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 energyconserving mesh sampling and weighting [19, 20] hyper reduction method, which approximates the fullorder 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 \(\varvec{V}\) (its construction will be discussed in section “Construction of reducedorder model components”)
the dimensionally reduced model (DROM) retrieved from the Galerkin projection reads
with \({\hat{\varvec{d}}} \in \mathbb {R}^n\) assuming \(n \ll N\). The argument of the residual is restricted to the column span of the ROB \(\varvec{V} \hat{\varvec{d}} \in \mathrm {span}(\varvec{V})\), which corresponds to a reduction of the number of DOFs. Consistently, the number of equations is reduced by multiplication with the transposed ROB \(\varvec{V}^T {\varvec{r}}(\varvec{V} {\hat{\varvec{d}}}) \in \mathbb {R}^n\). As a result, application of the NewtonRaphson iteration scheme leads to
wherein \(\varvec{J}_r\) is the residual Jacobian with respect to the displacement field \(\varvec{d}\). Equations (38), (39) reveal that only lowdimensional 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 \(\varvec{r}({\varvec{V}} {\hat{\varvec{d}}})\) still needs to be evaluated together with its Jacobian \(\varvec{J}_r(\varvec{V}{\hat{\varvec{d}}})\) throughout NewtonRaphson 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 energyconserving 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 \({\varvec{f}}_{\mathrm {int}} \in {\mathbb {R}}^{N}\) with a surrogate \({\tilde{\varvec{f}}}_{\mathrm {int}} \in {\mathbb {R}}^N\), which will result in an accurate approximation after projection, that is
\({\tilde{\varvec{f}}}_{\mathrm {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 \(\varvec{d}, \delta \varvec{d}\in \mathrm {span}(\varvec{V})\)
Applying a sum over all element internal force contributions we can rewrite
using \(\varvec{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 approximation \(\tilde{W}_{\mathrm {int}}(\varvec{d}, \delta \varvec{d})\) of the internal virtual work such that
with
In contrast to (43), (45) only contains a summation over a reduced element set \(\tilde{\mathcal {E}}\). Additional nonnegative element weights \(w^{(e)} \in {\mathbb {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 \(\tilde{\mathcal {E}}\) 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 \(\mathcal {\hat{S}}\)
with
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 reducedorder model components”). Note that \(\mathcal {S}\) consists of m (even number) modes corresponding to virtual prestress displacement modes \(\varvec{d}_{p(i)}\) and the sum of virtual prestress displacement and deformation stage displacement modes \(\varvec{d}_{p(i)} + \varvec{d}_{d(i)}\) with \(i \in \{0, \ldots , \frac{m}{2}  1\}\).
In its unassembled shape, the restriction of Eq. (44) to \(\forall {\hat{\varvec{d}}} \in \hat{\mathcal {S}}\) reads
In order to keep the cardinality of the reduced element set \(\tilde{\mathcal {E}}\) low, approximation (48) has to be accurate with a minimal number of nonzero weights. The corresponding optimization problem is
The zeronorm \(\left\Vert (\bullet ) \right\Vert _0\) counts the number of nonzero entries and is used as the objective function applied to the vector of element weights \(\varvec{w}\), which is constrained to have nonnegative values expressed by its minimum entry \(\min (\varvec{w})\) being nonnegative. This constraint is required in order to ensure a positive semidefinite Jacobian of the internal force vector [19]. The other constraint is a fulfillment of Eq. (48) up to the relative tolerance \(\varepsilon _h\). Consequently,
with vector valued entries
and
with vector valued entries
wherein \(i \in \{0, \ldots , m1\}, j \in \{0, \ldots , \mathcal {E}1 \}\) and \({\hat{\varvec{d}}}_i\) are vectors from the set \(\hat{\mathcal {S}}\).
Optimization problem (49) can be approximately solved with a sparse nonnegative leastsquares solver. Thereby sparse refers to the solution vector \(\varvec{w}\), in the sense that the number of nonzero 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 nonzero weights. As a consequence, the hyper reduced internal force vector and its Jacobian read
with the element stiffness \(\varvec{J}^{(e)}_{\mathrm {int}} = \frac{\partial \varvec{f}^{(e)}_{\mathrm {int}}}{\partial \varvec{d}^{(e)}}\).
We denote the dimensionally reduced as well as hyper reduced model (DHROM) as
and state the corresponding NewtonRaphson iterations
wherein \({\tilde{\varvec{r}}}\) is a residual approximation using \({\tilde{\varvec{f}}}_{\mathrm {int}}\) and \({\tilde{\varvec{J}}}_r\) is the corresponding hyper reduced residual Jacobian.
Construction of reducedorder model components
In the given manyquery context, the residual (34) depends on a modifiable set of model parameters. Introducing a parameter vector
with \([\mathrm {lb}_i; \mathrm {ub}_i] \ni \mu _i\) being the lower and upper bounds for parameter \(\mu _i\), we extend the notation of Eq. (35) to
and attempt to find a ROB \(\varvec{V}\) and a set of ECSW modes \(\mathcal {S}\) such that the resulting DROM as well as DHROM will accurately approximate the FOM solution for all \({\varvec{\mu }} \in \mathcal {P}\).
A prerequisite for an accurate ROM is that the FOM solution \({\varvec{d}}({\varvec{\mu }})\) can be accurately represented within the column span of the ROB
This motivates a datadriven 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 socalled snapshot matrix
with \(\varvec{S} \in {\mathbb {R}}^{N \times m}\).
Two datadriven approaches for the construction of the ROB have gained special interest in projectionbased MOR. Proper orthogonal decomposition (POD) can be used to orthogonalize \(\varvec{S}\) yielding a ROB \(\varvec{V}_{\mathrm {pod}} \in \mathbb {R}^{N \times m}\) such that [17]
wherein \(\varvec{V}^{n}_{\mathrm {pod}}\) corresponds to a selection of the first n columns of \(\varvec{V}_{\mathrm {pod}}\) with \(n \le m\) and \(\left\Vert (\bullet ) \right\Vert _F\) is the Frobenius norm. Consequently, POD is used whenever solution snapshots can be accurately represented by a lowdimensional 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
Spacefilling 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 spacefilling property of a design or to construct spacefilling designs by optimization of that criterion. A greedy version with reduced computational complexity is presented in [39] under the name “CoffeeHouse 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 \({\varvec{\mu }}\), 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 \(\Sigma _i \subset \mathcal {P}\) as subset of the parameter space and a set of previously chosen points \(\Sigma _c \subset \mathcal {P}\), the next point \({\varvec{\mu }} \in \Sigma _i\) is chosen such that the minimal distance to a neighboring point \({\varvec{p}} \in \Sigma _c\) is maximized in a reference hypercube. Thereby \(\chi \) and \(\chi ^{1}\) map from physical domain to reference hypercube and vice versa, respectively.
The steps for a greedy MMD design on a training grid \(\Sigma _t \subset \mathcal {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 radialbasisfunction approximation surrogate. In [46], MMD sampling is used to distribute points in Voronoi cells for multifidelity radialbasisfunction metamodeling.
Subspace angles for the comparison of subspaces
Subspace angles (or principal angles) are a concept from matrix computations [47]. Given two matrices \(\varvec{Y} \in \mathbb {R}^{N\times n}, \varvec{Z} \in \mathbb {R}^{N\times m}\) with \(n \le m\), subspace angles can be defined recursively as the minimum value
while the corresponding principal vectors follow from the minimization arguments
The sets
and
depend on the minimization arguments \({\varvec{y}}_j\) and \({\varvec{z}}_j\) of the previous k iterations.
A maximum subspace angle of \(0^{\circ }\) indicates that \(\mathrm {span}(\varvec{Y})\) \(\subseteq \mathrm {span}(\varvec{Z})\), while a maximum subspace angle of \(90^{\circ }\) indicates that there is at least one direction in \(\mathrm {span}(\varvec{Y})\) which is orthogonal to \(\mathrm {span}(\varvec{Z})\). More general, the maximum subspace angle can be interpreted as a distance measure from \(\mathrm {span}(\varvec{Y})\) to \(\mathrm {span}(\varvec{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.
In projectionbased 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 DiffusionConvectionReaction 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.
Algorithm 4 exposes the individual steps. After having calculated an initial ROB (line 2:) from the local snapshot matrix \(\varvec{s}\) at the initial parametric configuration \({\varvec{\mu }}\) (line 1:), the algorithm iterates over a (predefined) tuple \(\Sigma \) (line 5:) of subdomains \(\Sigma _{sd,i} \subset \Sigma _t\) for \(i \in \{0, \ldots n_{sd}1\}\) (\(n_{sd}\) consequently is the number of subdomains), wherein \(\Sigma _{sd,i} \cap _{i \ne j} \Sigma _{sd,j} = \emptyset \ \mathrm {for} \ i,j \in \{0, \ldots , n_{sd}1\}\) and \(\cup _{i=0}^{n_{sd}1} \Sigma _{sd,i} = \Sigma _t\). In every iteration, a parametric configuration is chosen within \(\Sigma _{sd,i}\) by a greedy MMD step (line 8:). Note that in line 8: distances to all previously selected points \(\Sigma _c\) are taken into account, although the new point is selected exclusively from \(\Sigma _{sd,i}\). Subdomains are excluded from sampling (line 13:) depending on the threshold \(\alpha _m\) (line 12:). The algorithm stops, if all subdomains have been excluded (line 19:,20:). On output, Algorithm 4 returns a set of selected grid points \(\Sigma _c\), a global orthogonal ROB \(\varvec{V}\) as well as a globally collected snapshot matrix \(\varvec{S}\). We use the ROB for dimensional reduction by the Galerkin projection (37), while the snapshot matrix is used to compute the set \(\mathcal {S}\) (47). Consequently, the number of ECSW displacement modes \(\mathcal {S}\) coincides with the dimension of the subspace \(\mathrm {span}(\varvec{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 \(\varvec{s}(\varvec{\mu })\) in every greedy maximin iteration. As a consequence, we can state that \(\mathrm {span}(\varvec{V}_i) \subseteq \mathrm {span}(\varvec{V}_j) \ \mathrm {for} \ j \ge i\), wherein \(\varvec{V}_{i}\) and \(\varvec{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 finegrained. 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 timeinvariant statespace 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 patientspecific 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.
Patientspecific computational models
Figures 4, 5 and 6 visualize the computational mesh, a cut through the computational domain depicting a separation between the ILT and the aortic wall and exemplary von Mises stress distributions for patient 1, patient 2 and patient 3, respectively.
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 Fbar element technology [57]. Table 1 depicts information on the model discretization.
Referring to section “Computational modeling of abdominal aortic aneurysms”, we quantify boundary conditions by a diastolic blood pressure of \(p_{\mathrm {dia}} = 87 \ \mathrm {mmHg}\) (\(11.6 \mathrm {kPa}\)) and a systolic blood pressure of \(p_{\mathrm {sys}} = 121 \ \mathrm {mmHg}\) (\(16.1 \ \mathrm {kPa}\)). The ILT stiffness c is interpolated linearly from a luminal stiffness of \(c = 2.62 \ \mathrm {kPa}\) to a medial stiffness of \(c = 1.98 \ \mathrm {kPa}\) and from the medial stiffness to an abluminal stiffness of \(c = 1.73 \ \mathrm {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 patientspecific Lognormal probability distributions from [58] for each entry of the parameter vector \({\varvec{\mu }}\). In more detail, the parameter domain bounds are chosen as
with
being the ppercentile value for a Lognormal distribution with expectation \(\mu _{\gamma }\) and standard deviation \(\sigma _{\gamma }\). \(\mathrm {erf}\) denotes the error function. Consequently, the range within the chosen parameter domain bounds covers \(95\%\) of realizations of \({\varvec{\mu }}\).
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 oneshot (i. e. no adaptation) design distributing 200 points in the parameter domain (Algorithm 4 with \(\Sigma _{sd} = \{\Sigma _{sd,0}\}, \alpha _m = 0.0\) and stopping at \(\Sigma _c = 200\)). If a simulation fails to converge, a neighboring point is taken in the set of selected points \(\Sigma _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 “minimumvalue” point \({\varvec{\mu }} = [\alpha _l, \beta _l, t_l]^T\) (see Table 2) for each patientspecific 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. \(\alpha< \alpha _m, \beta< \beta _m, t < t_m\)) in blue and SADs corresponding to \((+++)\)octant configurations of the parameter domain (i.e. \(\alpha> \alpha _m, \beta> \beta _m, t > t_m\)) in dark red (given \(\alpha _m = \frac{1}{2}(\alpha _l + \alpha _u)\), \(\beta _m = \frac{1}{2}(\beta _l + \beta _u)\) and \(t_m = \frac{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 \(\alpha _m = 0.1\), \(\Sigma _{sd} = \{\Sigma _{sd,0} \ldots , \Sigma _{sd,7}\}\) and the initial configuration \({\varvec{\mu }} = [\alpha _l, \beta _l, t_l]^T\). Figure 8 depicts the decay of SADs for each patientspecific computational model. As one can observe, sampling runs until the last sample in each subdomain yields a SAD below \(\alpha _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 lowstiffness and thin vessel wall range of the parameter domain. This is plausible from a physical perspective, given that soft and thinwalled tissue will deform more than stiff and thickwalled tissue.
Patientspecific reducedorder models
We compute ROBs from greedy subdomain MMD sampling with parametric configurations as depicted in Fig. 9. A Galerkin projection (37) yields the patientspecific DROMs. Hyper reduction is achieved via ECSW (cf. section “Hyper reduction of internal force contribution”) parallelized on 4 processors with global tolerance set to \(\varepsilon _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 reducedorder model
We evaluate accuracy in terms of the relative error in the quantities of interest (von Mises stress field \({\varvec{\sigma }}_{vM}\) and von Mises strain field \({\varvec{e}}_{vM}\) in the aortic wall). The relative error is given as
wherein \({\varvec{x}} \in \{{\varvec{\sigma }}_{vM}, {\varvec{e}}_{vM}\}\) corresponds to FOM quantities and \({\tilde{\varvec{x}}} \in \{{\tilde{\varvec{\sigma }}}_{vM}, {\tilde{\varvec{e}}}_{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 \(\alpha _m\) and \(\varepsilon _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 \in \{\sigma _{vM}, e_{vM}\}\) denoting von Mises stress or von Mises strain as the quantity of interest and \(RE_i({\tilde{\varvec{x}}}, {\varvec{x}})\) being the relative error (71) of sample i.
Figure 13 depicts results for patient 1 as exemplary model. In more detail, Fig. 13a corresponds to DROM results with ROBs created by Algorithm 4 at different SADs, while Fig. 13b corresponds to DHROM results with ECSW meshes at different tolerances and an unchanged ROB at \(\alpha _m = 0.1\). As can be observed from Fig. 13a, lower values for \(\alpha _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 \(\varepsilon _h\) lead to an increased number of selected mesh elements with decreasing mean relative error indicating a more accurate DHROM. We conclude, that both \(\alpha _m\) and \(\varepsilon _h\) strongly influence ROM accuracy.
Monte Carlo sampling on the reducedorder 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 patientspecific Lognormal 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 kerneldensityestimated probability distributions gained from FOM and DHROM samples. We apply the Gaussian kerneldensityestimator 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 patientspecific computational models as well as corresponding speedups in Table 6. All simulations in this section were performed on a workstation with Intel Xeon W2133 (3.60GHz) processors.
The values in Table 6 are mean values corresponding to seven simulations (per patient) evaluated at face midpoints as well as the midpoint of the patientspecific 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 fullorder residual as well as its Jacobian need to be evaluated. A rather 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 hyperreduction by ECSW (as described in “Patientspecific reducedorder 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 projectionbased MOR of patientspecific AAA models. A dimensionally reduced model was built by a Galerkin projection on a lowdimensional subspace and a dimensionally reduced as well as hyper reduced model was built by the Galerkin projection and energyconserving 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 lowdimensional solution subspaces. Therein, a stopping criterion based on subspace angles and an exclusion of subdomains was applied. Finally, three patientspecific computational examples with different complexities were demonstrated.
The proposed sampling algorithm led to comparable results in terms of reducedorder basis size as well as number of sampled mesh elements for all three patientspecific 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 patientspecific 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 kerneldensityestimated 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 fullorder model.
Availability of data and materials
All data referring to the provided examples is contained in the manuscript.
Abbreviations
 MOR:

Model order reduction
 AAA:

Abdominal aortic aneurysm
 FOM:

Fulloder model
 ROM:

Reducedorder model
 MULF:

Modified updated Lagrangian formulation
 ROB:

Reducedorder basis
 DOF:

Degree of freedom
 DROM:

Dimensionally reduced model
 DHROM:

Dimensionally reduced and hyper reduced model
 ECSW:

Energyconserving mesh sampling and weighting
 ILT:

Intraluminal thrombus
 PVW:

Principle of virtual work
 POD:

Proper orthogonal decomposition
 MMD:

Maximin distance
 SAD:

Subspace angle distance
References
 1.
Formaggia L, Quarteroni A, Veneziani A. Cardiovascular Mathematics: Modeling and simulation of the circulatory system, vol. 1. : Springer Science & Business Media; 2010.
 2.
Leach JR, Kao E, Zhu C, Saloner D, Hope MD. On the Relative Impact of Intraluminal Thrombus Heterogeneity on Abdominal Aortic Aneurysm Mechanics. Journal of Biomechanical Engineering. 2019;141(11).
 3.
Polzer S, Gasser TC. Biomechanical rupture risk assessment of abdominal aortic aneurysms based on a novel probabilistic rupture risk index. Journal of The Royal Society Interface. 2015;12(113):20150852.
 4.
Hemmler A, Lutz B, Kalender G, Reeps C, Gee MW. Patientspecific in silico endovascular repair of abdominal aortic aneurysms: application and validation. Biomechanics and modeling in mechanobiology. 2019;18(4):983–1004.
 5.
Perrin D, Badel P, Orgeas L, Geindreau C, rolland du Roscoat S, Albertini JN, et al. Patientspecific simulation of endovascular repair surgery with tortuous aneurysms requiring flexible stentgrafts. Journal of the mechanical behavior of biomedical materials. 2016;63:86–99.
 6.
Lee LC, Ge L, Zhang Z, Pease M, Nikolic SD, Mishra R, et al. Patientspecific finite element modeling of the Cardiokinetix Parachute® device: effects on left ventricular wall stress and function. Medical & biological engineering & computing. 2014;52(6):557–66.
 7.
Niestrawska JA, Viertler C, Regitnig P, Cohnert TU, Sommer G, Holzapfel GA. Microstructure and mechanics of healthy and aneurysmatic abdominal aortas: experimental analysis and modelling. Journal of The Royal Society Interface. 2016;13(124):20160620.
 8.
Ockert S, Boeckler D, Allenberg J, Schumacher H. Rupturiertes abdominelles aortenaneurysma. Gefaesschirurgie. 2007;12(5):379–91.
 9.
Raghavan M, Vorp DA. Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysm: identification of a finite strain constitutive model and evaluation of its applicability. Journal of biomechanics. 2000;33(4):475–82.
 10.
Biehler J, Gee MW, Wall WA. Towards efficient uncertainty quantification in complex and largescale biomechanical problems based on a Bayesian multifidelity scheme. Biomechanics and modeling in mechanobiology. 2015;14(3):489–513.
 11.
Vorp DA. Biomechanics of abdominal aortic aneurysm. Journal of biomechanics. 2007;40(9):1887–902.
 12.
Maier A. Computational modeling of rupture risk in abdominal aortic aneurysms. Technische Universität München. 2012.
 13.
Gee M, Reeps C, Eckstein H, Wall W. Prestressing in finite deformation abdominal aortic aneurysm simulation. Journal of biomechanics. 2009;42(11):1732–9.
 14.
Gee MW, Förster C, Wall W. A computational strategy for prestressing patientspecific biomechanical problems under finite deformation. International Journal for Numerical Methods in Biomedical Engineering. 2010;26(1):52–72.
 15.
Reeps C, Maier A, Pelisek J, Härtl F, GrabherMeier V, Wall W, et al. Measuring and modeling patientspecific distributions of material properties in abdominal aortic aneurysm wall. Biomechanics and modeling in mechanobiology. 2013;12(4):717–33.
 16.
Raghavan ML, Hanaoka MM, Kratzberg JA, de Lourdes Higuchi M, Da Silva ES. Biomechanical failure properties and microstructural content of ruptured and unruptured abdominal aortic aneurysms. Journal of biomechanics. 2011;44(13):2501–7.
 17.
Quarteroni A, Manzoni A, Negri F. Reduced basis methods for partial differential equations: an introduction, vol. 92. : Springer; 2015.
 18.
Carlberg K, Barone M, Antil H. Galerkin v. leastsquares PetrovGalerkin projection in nonlinear model reduction. Journal of Computational Physics. 2017;330:693–734.
 19.
Farhat C, Avery P, Chapman T, Cortial J. Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energybased mesh sampling and weighting for computational efficiency. International Journal for Numerical Methods in Engineering. 2014;98(9):625–62.
 20.
Farhat C, Chapman T, Avery P. Stability and accuracy properties of the energyconserving sampling and weighting (ECSW) method for the hyper reduction of nonlinear finite element dynamic models. Int J Numer Methods Eng. 2015;102:1077–110.
 21.
Quicken S, Donders WP, van Disseldorp EM, Gashi K, Mees BM, van de Vosse FN, et al. Application of an adaptive polynomial chaos expansion on computationally expensive threedimensional cardiovascular models for uncertainty quantification and sensitivity analysis. Journal of biomechanical engineering. 2016;138(12):121010.
 22.
Sankaran S, Marsden AL. A stochastic collocation method for uncertainty quantification and propagation in cardiovascular simulations. Journal of biomechanical engineering. 2011;133(3):031001.
 23.
Biehler J, Wall W. The impact of personalized probabilistic wall thickness models on peak wall stress in abdominal aortic aneurysms. International journal for numerical methods in biomedical engineering. 2018;34(2):e2922.
 24.
Chang GH, Schirmer CM, ModarresSadeghi Y. A reducedorder model for wall shear stress in abdominal aortic aneurysms by proper orthogonal decomposition. Journal of biomechanics. 2017;54:33–43.
 25.
Negri F. Efficient Reduction Techniques for the Simulation and Optimization of Parametrized Systems. Ecole Polytechnique Fédérale de Lausanne. 2015.
 26.
Gasser TC, Görgülü G, Folkesson M, Swedenborg J. Failure properties of intraluminal thrombus in abdominal aortic aneurysm under static and pulsating mechanical loads. Journal of vascular surgery. 2008;48(1):179–88.
 27.
Doll S, Schweizerhof K. On the development of volumetric strain energy functions. J Appl Mech. 2000;67(1):17–21.
 28.
De Putter S, Wolters B, Rutten M, Breeuwer M, Gerritsen F, Van de Vosse F. Patientspecific initial wall stress in abdominal aortic aneurysms with a backward incremental method. Journal of biomechanics. 2007;40(5):1081–90.
 29.
Lu J, Zhou X, Raghavan ML. Inverse elastostatic stress analysis in predeformed biological structures: demonstration using abdominal aortic aneurysms. Journal of biomechanics. 2007;40(3):693–6.
 30.
Hemmler A, Lutz B, Reeps C, Kalender G, Gee MW. A methodology for in silico endovascular repair of abdominal aortic aneurysms. Biomechanics and modeling in mechanobiology. 2018;17(4):1139–64.
 31.
Rutzmoser J. Model Order Reduction for Nonlinear Structural Dynamics. Technische Universität München. 2018.
 32.
Carlberg K, Tuminaro R, Boggs P. Preserving Lagrangian structure in nonlinear model reduction with application to structural dynamics. SIAM Journal on Scientific Computing. 2015;37(2):B153–84.
 33.
Haasdonk B, Ohlberger M. Reduced basis method for finite volume approximations of parametrized linear evolution equations. ESAIM: Mathematical Modelling and Numerical Analysis. 2008;42(2):277–302.
 34.
Quarteroni A, Rozza G, Manzoni A. Certified reduced basis approximation for parametrized partial differential equations and applications. Journal of Mathematics in Industry. 2011;1(1):3.
 35.
Hesthaven JS, Stamm B, Zhang S. Efficient greedy algorithms for highdimensional parameter spaces with applications to empirical interpolation and reduced basis methods. ESAIM: Mathematical Modelling and Numerical Analysis. 2014;48(1):259–83.
 36.
Maday Y, Stamm B. Locally adaptive greedy approximations for anisotropic parameter reduced basis spaces. SIAM Journal on Scientific Computing. 2013;35(6):A2417–41.
 37.
Jiang J, Chen Y, Narayan A. Offlineenhanced reduced basis method through adaptive construction of the surrogate training set. Journal of Scientific Computing. 2017;73(2–3):853–75.
 38.
Johnson ME, Moore LM, Ylvisaker D. Minimax and maximin distance designs. Journal of statistical planning and inference. 1990;26(2):131–48.
 39.
Müller WG. Coffeehouse designs. In: Atkinson A, Bogacka B, Zhigljavsky AA, editors. Optimum design 2000. : Springer; 2001. p. 241–8.
 40.
Pronzato L. Minimax and maximin spacefilling designs: some properties and methods for construction. Journal de la Societe Francaise de Statistique. 2017;158(1):7–36.
 41.
Yondo R, Bobrowski K, Andrés E, Valero E. A review of surrogate modeling techniques for aerodynamic analysis and optimization: current limitations and future challenges in industry. In: Advances in Evolutionary and Deterministic Methods for Design, Optimization and Control in Engineering and Sciences. Springer; 2019. p. 19–33.
 42.
Van Der Herten J, Van Steenkiste T, Couckuyt I, Dhaene T. Surrogate Modelling with Sequential Design for Expensive Simulation Applications. Computer Simulation. 2017;p.;173.
 43.
Crombecq K, Laermans E, Dhaene T. Efficient spacefilling and noncollapsing sequential design strategies for simulationbased modeling. European Journal of Operational Research. 2011;214(3):683–96.
 44.
Garud SS, Karimi IA, Kraft M. Design of computer experiments: A review. Computers & Chemical Engineering. 2017;106:71–95.
 45.
Liu H, Hervas JR, Ong YS, Cai J, Wang Y. An adaptive RBFHDMR modeling approach under limited computational budget. Structural and Multidisciplinary Optimization. 2018;57(3):1233–50.
 46.
Cai X, Qiu H, Gao L, Wei L, Shao X. Adaptive radialbasisfunctionbased multifidelity metamodeling for expensive blackbox problems. AIAA Journal. 2017;p. 2424–2436.
 47.
Golub GH, Van Loan CF. Matrix computations. 3rd ed. : Johns Hopkins University Press; 1996.
 48.
Lieu T, Lesoinne M. Parameter adaptation of reduced order models for threedimensional flutter analysis. In: 42nd AIAA Aerospace Sciences Meeting and Exhibit;2004.
 49.
Lieu T, Farhat C, Lesoinne M. Reducedorder fluid/structure modeling of a complete aircraft configuration. Computer methods in applied mechanics and engineering. 2006;195(41–43):5730–42.
 50.
Lieu T, Farhat C. Adaptation of aeroelastic reducedorder models and application to an F16 configuration. AIAA journal. 2007;45(6):1244–57.
 51.
Akman T. Local improvements to reducedorder approximations of optimal control problems governed by diffusionconvectionreaction equation. Computers & Mathematics with Applications. 2015;70(2):104–31.
 52.
Bazaz MA, Nahve S, Nabi M, Janardhanan S, Rehman M. Adaptive parameter space sampling in matrix interpolatory pMOR. In: 2015 International Conference on Recent Developments in Control, Automation and Power Engineering (RDCAPE). IEEE; 2015. p. 83–9.
 53.
Varona MC, Lohmann B, Nabi M. Automatic adaptive sampling in parametric model order reduction by matrix interpolation. In: 2017 IEEE International Conference on Advanced Intelligent Mechatronics (AIM). IEEE; 2017. p. 472–7.
 54.
Kalra S, Nabi M. TPWL Simulation of Large Nonlinear Circuits Using Subspace Angle Based Adaptive Sampling. IEEE Transactions on Circuits and Systems II: Express Briefs. 2019;67(3):575–9.
 55.
Martufi G, Lindquist Liljeqvist M, Sakalihasan N, Panuccio G, Hultgren R, Roy J, et al. Local diameter, wall stress, and thrombus thickness influence the local growth of abdominal aortic aneurysms. Journal of Endovascular Therapy. 2016;23(6):957–66.
 56.
Bruder L, Pelisek J, Eckstein HH, Gee MW. Biomechanical rupture risk assessment of abdominal aortic aneurysms using clinical data: A patientspecific, probabilistic framework and comparative casecontrol study. PloS one. 2020;15(11):e0242097.
 57.
de Souza Neto E, Perić D, Dutko M, Owen D. Design of simple low order finite elements for large strain analysis of nearly incompressible solids. International Journal of Solids and Structures. 1996;33(20–22):3277–96.
 58.
Biehler J, Kehl S, Gee MW, Schmies F, Pelisek J, Maier A, et al. Probabilistic noninvasive prediction of wall properties of abdominal aortic aneurysms using Bayesian regression. Biomechanics and modeling in mechanobiology. 2017;16(1):45–61.
 59.
Heroux MA, Willenbring JM. Trilinos users guide. Sandia National Laboratories. 2003.
 60.
Davis TA. Algorithm 832: UMFPACK V4. 3–an unsymmetricpattern multifrontal method. ACM Transactions on Mathematical Software (TOMS). 2004;30(2):196–9.
 61.
Virtanen P, Gommers R, Oliphant TE, et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods. 2020;17:261–72.
Acknowledgements
The authors gratefully acknowledge support and funding by the Leibniz Rechenzentrum München (LRZ) of the Bavarian Academy of Sciences under the contract pn34cu.
Funding
Open Access funding enabled and organized by Projekt DEAL. Funding was provided by Bayerische Akademie der Wissenschaften (DE), LRZ (Grant No. pn34cu).
Author information
Affiliations
Contributions
AS developed the idea, performed the software implementation, conducted numerical experiments and wrote the manuscript. MWG finetuned the research idea, suggested numerical experiments and revised the paper. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Schein, A., Gee, M.W. Greedy maximin distance sampling based model order reduction of prestressed and parametrized abdominal aortic aneurysms. Adv. Model. and Simul. in Eng. Sci. 8, 18 (2021). https://doi.org/10.1186/s40323021002037
Received:
Accepted:
Published:
Keywords
 Abdominal aortic aneurysm
 Nonlinear model order reduction
 Prestressing
 Finite element method