Fast parametric analysis of trimmed multi-patch isogeometric Kirchhoff-Love shells using a local reduced basis method

This contribution presents a model order reduction framework for real-time efficient solution of trimmed, multi-patch isogeometric Kirchhoff-Love shells. In several scenarios, such as design and shape optimization, multiple simulations need to be performed for a given set of physical or geometrical parameters. This step can be computationally expensive in particular for real world, practical applications. We are interested in geometrical parameters and take advantage of the flexibility of splines in representing complex geometries. In this case, the operators are geometry-dependent and generally depend on the parameters in a non-affine way. Moreover, the solutions obtained from trimmed domains may vary highly with respect to different values of the parameters. Therefore, we employ a local reduced basis method based on clustering techniques and the Discrete Empirical Interpolation Method to construct affine approximations and efficient reduced order models. In addition, we discuss the application of the reduction strategy to parametric shape optimization. Finally, we demonstrate the performance of the proposed framework to parameterized Kirchhoff-Love shells through benchmark tests on trimmed, multi-patch meshes including a complex geometry. The proposed approach is accurate and achieves a significant reduction of the online computational cost in comparison to the standard reduced basis method.


Introduction
In the past decades, the integration of geometric design and finite element analysis has attracted a lot of attention within the computational engineering community.The introduction of Isogeometric Analysis (IGA) [1] paved the way for an integrated framework from Computer Aided Design (CAD) to numerical simulation followed by a wide range of successful applications in several fields.Adopting the isogeometric paradigm, the effort involved in meshing and geometry clean-up can be circumvented by employing the same representation for both geometric description and numerical analysis.The prevailing technology for geometric representation in CAD is based on splines, namely B-splines and in particular non-uniform rational B-splines (NURBS).The reader is further referred to [1,2] for a detailed overview of the method.
In the context of shell analysis, IGA offers clear advantages.In particular, Kirchhoff-Love shell formulations are based on fourth order partial differential equations (PDEs) and require C 1 -continuity that hinders the use of classical C 0continuous finite elements.Due to the higher continuity of splines, isogeometric discretizations are well suited for the modeling of higher-order PDEs.An overview of isogeometric methods for Kirchhoff-Love shell analysis is given in [3].
Nevertheless, real-world applications of complex geometries require the treatment of further issues.Typically, complex shapes are represented in CAD using Boolean operations such as intersection and union.This results in trimmed meshes, where the underlying discretization is unfitted with respect to the physical object.First, tackling trimmed surfaces needs to be properly addressed in the numerical simulation.A detailed review on trimming and related challenges such as numerical integration, conditioning, and others, is given in [4].We also refer to [5,6] in the context of trimmed shell analysis.Besides trimming, complex shapes are commonly represented by multiple, non-conforming patches in CAD and require suitable coupling strategies.In particular for shells, C 1 -continuity between adjacent patches is required for the analysis.There are several methods in the literature to achieve displacement and rotational continuity in a weak sense, such as the penalty, Nitsche's, and mortar methods.A comprehensive overview of these methods in the context of IGA is given in [7].We also refer to previous works on mortar methods [8][9][10][11], penalty approaches [12][13][14], and Nitsche's method [15][16][17][18] including their application to the analysis of trimmed shells.For the case of conforming patches, C 1 -continuity across interfaces can be even imposed strongly, as shown recently in [19].
On the other hand, there exist several applications where the analysis needs to be performed rapidly for many different parametric configurations.This is, for example, the case in design and shape optimization, where the geometry is updated on the fly within the optimization loop.The reader is further referred to [20][21][22], for an overview on isogeometric shape optimization.To achieve a computational speedup, research efforts have been devoted to the development of efficient reduced order models (ROMs) for the solution of parameterized problems.In particular, we refer to previous works in the context of combining IGA with ROMs [23][24][25].Most of these works employ the reduced basis method to construct projection-based ROMs in combination with hyper-reduction techniques such as the Empirical Interpolation Method [26] to handle geometrical parameterizations for non-affine problems.A detailed overview on reduced basis methods is given in [27,28].We also refer to [29,30] in the context of isogeometric ROMs for multi-patch geometries.Nevertheless, the study of efficient ROMs in the context of isogeometric shell analysis is a subject that still lacks thorough investigation.The application of reduced basis methods to isogeometric Kirchhoff-Love shells has been explored in [31] for simple geometries.However, a general ROM framework suitable for trimmed, multi-patch geometries of industrial relevance is still missing.
This contribution aims to investigate the use of efficient ROM strategies for fast isogeometric Kirchhoff-Love shell analysis formulated on parameterized trimmed and multi-patch geometries.Based on our previous work [32], we employ a local reduced basis method to construct efficient ROMs for problems formulated on parameterized unfitted geometries obtained by trimming operations.Local ROMs have been a subject of previous research works within the model reduction community [33][34][35][36].The method we use in this work is based on extension and parameter-based clustering of snapshots to apply reduction techniques such as the Proper Orthogonal Decomposition and Discrete Empirical Interpolation Method [37,38].
We recast the Kirchhoff-Love shell formulation into parameterized trimmed geometries within a multi-patch setting.The proposed reduction framework enables an efficient offline/online decomposition and is agnostic to the underlying coupling method, therefore different strategies can be in principle considered to enforce displacement and rotational continuity in a weak sense.In addition, we discuss the application of the reduction framework to parametric shape optimization problems.We also refer to previous works on ROM-based parametric optimization [39][40][41].
The manuscript is structured as follows: Section 2 provides a brief overview of basic concepts related to B-splines and trimming formulated on parameterized domains.In Section 3 we present the parameterized formulation for the Kirchhoff-Love shell, while Section 4 provides the necessary definitions related to the coupling of multipatches.We discuss the local reduced basis method in Section 5 and its application to parametric shape optimization in Section 6.In Section 7 we study several numerical experiments including non-conforming discretizations and a complex geometry to assess the performance of the proposed framework.The main conclusions that can be drawn from this study are finally summarized in Section 8.

Parameterized trimmed geometries in isogeometric analysis
In this section, we provide a brief review of some basic concepts related to B-splines and trimming in isogeometric analysis.In particular, we formulate the trimming operation in the context of parameterized geometries for single-patch domains.For a detailed review of splines in the context of isogeometric analysis, the reader is referred to [1,2,42].Moreover, we refer to [4] for a comprehensive review of trimming in isogeometric analysis.

B-spline basis functions
To illustrate the basic concept of B-splines, we introduce a knot vector in the parametric space [0, 1] as a non-decreasing sequence of real values denoted by Ξ = {ξ 1 , . . ., ξ n+p+1 }.Here, the integer n denotes the number of basis functions and p the degree.We further introduce a univariate B-spline basis function b j ij ,pj where p j denotes the degree and i j the index of the function in the j-th parametric direction.The B-spline functions B i,p (ξ) can be easily defined in multiple dimensions exploiting the tensor product of univariate B-splines: where d is the dimension of the parametric space.Moreover, the vector i = (i 1 , . . ., i d) is a multi-index denoting the position in the tensor-product structure and p = (p 1 , . . ., p d) the polynomial degree corresponding to the parametric coordinate ξ = (ξ 1 , . . ., ξ d).For simplicity, we will assume from now on that the vector p is identical in all parametric directions and can be replaced by a scalar value p.The B-spline basis is C p−k -continuous at every knot, where k is the multiplicity of the knot, and is C ∞ elsewhere.The concept of B-splines can be easily extended to rational B-splines and the interested reader is further referred to [2] for a detailed exposition.

Parameterized trimming
In the following, we briefly present the mathematical foundation of trimming in isogeometric analysis and recast its formulation into the context of parameterized geometries.
For this purpose, we adopt the notation previously introduced in [43,44].
Let us first define the parameterized physical domain Ω(µ) ⊂ R d described by the geometrical parameters µ ∈ P ⊂ R M , where d is the dimension of the physical space for our problem, P is the space of parameters, and M is the number of parameters.We will show later how to obtain Ω(µ) through trimming operations.Now we consider the non-trimmed, single-patch physical domain Ω 0 ⊂ R d and its counterpart in the parameter space Ω0 = [0, 1] d.In the following, we assume that the non-trimmed domain Ω 0 is parameter-independent without loss of generality.In principle, the extension to the parameter dependent case is straightforward.The multi-patch setting will be discussed below in Section 4.
Given a control point mesh P i , a spline geometric map F : Ω0 → Ω 0 can be defined as As a result, by an abuse of notation, the non-trimmed physical domain is obtained by In what follows we are interested in parameterizing the trimmed regions.For this purpose, we introduce K parameter-dependent trimmed regions Ω 1 (µ), . . . ,Ω K (µ) ⊂ R d that are cut away by the trimming operation.The obtained physical domain can then be expressed as: The concept is illustrated in Fig. 1.Following this definition, the boundary of the domain consists of a trimmed part ∂Ω(µ)\∂Ω 0 and a part that corresponds to the original domain ∂Ω(µ) ∩ ∂Ω 0 .For the sake of simplicity, we will assume from now on that Dirichlet boundary conditions are not applied on the trimming boundary.
The interested reader may refer to [45,46] regarding the weak imposition of Dirichlet boundary conditions on the trimming boundary in combination with stabilization techniques.
Furthermore, it should be noted that the elements and basis functions are defined on the non-trimmed domain Ω 0 and the original domain remains unchanged by the trimming operation.Thus, we define the B-spline space S p 0 , which is constructed upon Ω 0 and is independent of the parameters µ, as In fact, let us now rewrite the trimming operation in the parametric domain Ω0 .
The counterpart of the parameter-dependent trimmed regions Ω i (µ), i = 1, . . ., K (see also (3)) in the parametric domain are defined as Ωi ⊂ Ω0 , while it holds that Ω i (µ) = F( Ωi ), i = 1, . . ., K. Thus, Equation (3) can be reinterpreted as Let us now introduce the B-spline space of degree p restricted to the trimmed domain Ω(µ): where supp(B i,p ) is the support of the non-trimmed basis functions.We remark that for numerical integration, a re-parameterization of cut elements is performed in the trimmed parametric domain Ω(µ).In this work we follow the high-order reparameterization procedure discussed in [47].The trimming operation can be understood as restricting the map F to an active region of the original domain, which is the visible part after the trimmed regions are cut away.Indeed, only the basis functions whose support intersects Ω(µ) are active.The remaining basis functions are inactive and do not contribute to the solution discretization of the problem.Moreover, the active basis functions and the dimension of the space may change for different values of the parameters µ.Therefore, to apply reduce order modeling techniques we will rely on the definition of the non-trimmed domain Ω 0 and associated B-spline space S p 0 following [32].This aspect will be further discussed in Section 5.
Hereinafter, for the sake of simplicity of exposition, and without loss of generality, we consider that no normal rotations are prescribed on the boundary, i.e., Γ D,θ (µ) = ∅, and that only homogeneous Dirichlet boundary conditions are applied on Γ D (µ) = Γ D,u (µ).
We further recall the spline geometric map F in (2) and construct a covariant basis, where the basis vectors are defined as Here F ,α denotes the partial derivative of the spline geometric mapping with respect to the α-th curvilinear coordinate.The midsurface normal vector a 3 can be constructed as a normalized cross-product of the two in-plane vectors a α , that is Then, we can construct the contravariant basis vectors that satisfy the Kronecker relationship as a α • a β = δ β α .Note that it holds a 3 = a 3 .The reader is further referred to [18] for an elaborate discussion on fundamentals of differential geometry that are relevant to the Kirchhoff-Love shell formulation.
Let us now introduce a discrete space V h (µ) to approximate our problem: where the spline space S p (Ω(µ)) 3 has at least C 1 -continuity, as typically required by isogeometric Kirchoff-Love formulations (see, e.g., [3]).Higher continuity requirements will be discussed in Section 4.2 for the case of Nitsche's interface coupling for multipatch domains.Now let us define the discrete weak formulation of the parameterized Kirchhoff-Love shell problem as: find The parameterized bilinear form a(•, •; µ) is given as: and the parameterized linear functional f (•; µ) reads: Note that θ n (v h ) = −n • ∇u ⊤ h a 3 is the normal rotation, n is the unit outward normal vector to the boundary Γ(µ), and the membrane and bending strain tensors are defined, respectively, as: where sym(•) denotes the symmetric part of a tensor, ∇ * is the surface gradient, and P = I − a 3 ⊗ a 3 is the in-plane projector, with I the identity tensor.
Let us recall that in Kirchhoff-Love shell kinematics we assume that the transverse shear strains vanish, i.e., the normal vectors remain straight and normal during deformation.The rotations are therefore constrained as follows: Now, let us define the energetically conjugate stresses for both the membrane and bending strains.For this purpose, we assume a linear elastic constitutive model and define the following fourth-order elasticity tensor: where repeated indices imply summation from 1 to 3, and E and ν are the elasticity modulus and Poisson's ratio, respectively.Assuming a constant thickness t and performing through-thickness integration we obtain the membrane and bending stresses as: For a rigorous derivation of the weak formulation the reader is further referred to [18].

Multi-patch geometries
Let us now consider that the parameterized domain Ω(µ) is split into non-overlapping subdomains, i.e., patches, such that where Ω k (µ k ) ∩ Ω l (µ l ) = ∅ for l ̸ = k.Here the patches Ω k (µ k ) are in principle trimmed and their definition follows the setting introduced in Section 2.2.For ease of exposition, we assume that the parameters associated to the k-th patch coincide with the parameters describing the global computational domain such that µ k = µ for k = 1, . . ., N p , although in principle different choices are possible (see [30] for more details).Furthermore, we define the common interface γ j (µ), between two adjacent patches such that where N Γ is the number of interfaces.We note that the interface can be both trimmed or non-trimmed.The B-spline space (4) can be now rewritten as where B k i,p and F k are the B-spline functions and the geometric map associated to the k-th patch, respectively.The dimension of the associated space is denoted as N k,0 = dim(S p k,0 ).And, analogously to (6), the B-spline space of degree p corresponding to the trimmed patch Ω k (µ) now reads where Ωk (µ) is the k-th parametric trimmed domain in (5).Following the choice ( 9), the multi-patch discrete space can be then written as As in the single-patch case, and unless stated otherwise, we assume the spaces V h,k (µ) to be C 1 -continuous within each patch.Nevertheless, the multi-patch space V h (µ) is discontinuous across the interfaces γ j (µ).The continuity required by problem (10) will be imposed in a weak sense.Thus, let us now define the coupling conditions for each interface γ j (µ).First, we denote the displacement fields restricted to ∂Ω k (µ) and ∂Ω l (µ) as u k (µ) and u l (µ), respectively.Then the coupling conditions can be expressed as that, by means of the standard jump operators we can rewrite as There exist several coupling strategies for the imposition of the above continuity constraints in a weak sense.In the numerical experiments of Section 7, we use both a super-penalty [13] and Nitsche's [18] methods to achieve displacement and rotational continuity, both described below.In principle, the reduced order modeling techniques to be applied (Section 5) are also suitable for other coupling strategies, such as the mortar method where Lagrange multipliers can be eliminated in the context of the reduced basis method [48].
The weak formulation (10) of the parameterized problem can be extended to the multi-patch case by enriching the bilinear form with additional terms that weakly enforce the coupling conditions.Thus, the shell multi-patch problem reads: find where a k and f k are the bilinear and linear parameterized forms ( 11) and ( 12), respectively, restricted to the patch k, while the coupling terms a Γ j will be further discussed in Sections 4.1 and 4.2.The discretization yields the following parameterized linear system of dimension where is the force vector, and N h (µ) is the global number of degrees of freedom.The problem ( 24)-( 25) is the high-fidelity or full order model (FOM) upon which we aim to construct a reduced order model (ROM) for fast parametric simulations.
In what follows, we describe the super-penalty [13] and Nitsche's [18] coupling methods (Sections 4.1 and 4.2, respectively).While the first approach is simpler and more efficient, its application is limited to the case in which each interface can be associated to a parametric face for at least one of the two contiguous patches at the interface.On the other hand, Nitsche's method overcomes such limitation, but is more convoluted as it involves further terms and third-order derivatives, and, consequently, requires C 2 -discretization spaces.

Projected super-penalty approach
In what follows, we will briefly present the projected super-penalty approach that we will use to impose the continuity constraints ( 23) for some of the numerical experiments in Section 7. The reader is further referred to [13,49] for a more detailed overview of the coupling approach.
At each interface γ j (µ), for j = 1, . . ., N Γ , let us first choose arbitrarily one of the neighboring patches as active.This active patch must be chosen such that γ j (µ) is fully contained in one of its four patch parametric faces.We extract the patch knot vector associated to that face and remove the first and last knots, denoting the resulting vector as Ξ j .Then, using Ξ j we construct a one-dimensional lower degree spline space S p−2 (γ j (µ)).Finally, we denote as Π j the L 2 -projection, related to the interface γ j (µ), onto the reduced vector space S p−2 (γ j (µ)) d for displacements, and onto S p−2 (γ j (µ)) for normal rotations.The discretized bilinear form is enriched with penalty terms that weakly enforce the coupling conditions (23) and exploit the properties of the L 2 -projection.Thus, these coupling terms read The parameters of the penalty method are chosen as: where h j denotes the interface mesh size and the measure |γ j (µ)| is the length of the coupling interface γ j (µ).The exponent factor c exp is chosen as c exp = p − 1 in the numerical experiments of Section 7, which yields optimal convergence of the method in the H 2 norm.A detailed discussion on the choice of c exp in relation to optimal convergence and conditioning of the underlying system of equations is given in [13].
We remark that this coupling approach is locking-free at interfaces and the choice of penalty parameters can be done automatically based on the problem setup.Nevertheless, the computational cost grows for higher order splines (p > 3) and the condition number related to the chosen penalty parameters may affect the accuracy of the solution.In addition, the coupling at interfaces where both sides are trimmed is still an open issue.To this end, we will consider the Nitsche's method for the coupling of patches in more general cases of complex geometries.Remark 1.To overcome such problem, in [13] the internal knots of the trimming interfaces are neglected for the computation of the intersection mesh at the interface, which can yield sub-optimal results in particular for non-smooth interfaces.

Nitsche's method
Let us now introduce the Nitsche's method for recovering C 1 -continuity at the multipatch interfaces.A detailed overview of the method is given in [18].The coupling conditions are imposed in a weak sense by augmenting the discretized bilinear form with penalty, consistency, and symmetry terms.Thus, the coupling terms in (24) read Here, ⟨•⟩ δ denotes the average operator, B nn = n • B(u h )n is the bending moment and T is the ersatz force, defined in [18] as where b = −∇a 3 is the curvature tensor, t is the counter-clockwise positively-oriented unit tangent vector to Γ(µ), and is the twisting moment.The average operator is defined as: where a is an arbitrary function defined over Ω(µ), and a| Ω k (µ) and a| Ω l (µ) its restrictions to the patches Ω k (µ) and Ω l (µ) at the interface γ j (µ).In the numerical experiments of Section 7, the penalty constants are chosen as We remark that this coupling approach is variationally consistent and stable (see Remark 3 below).In addition, its discretized weak formulation results in a wellconditioned system of linear system and is more robust with respect to the chosen parameters compared to penalty approaches.Nevertheless, it is easy to realize that the ersatz force T requires third-order derivatives (terms ∇ • B(u h ) and ∂B nt (u h )/∂t in Equation ( 29)), what implies the necessity of C 2 -continuous discretization spaces.Furthermore, the number of terms and the order of the involved derivatives make this method's implementation more convoluted than the super-penalty approach (recall Section 4.1).Remark 2. The ersatz force T in Equation (29) involves third-order derivatives (terms ∇ • B(u h ) and ∂B nt (u h )/∂t), what implies the necessity of C 2 -continuous discretization spaces.In particular, and as detailed in [18], the patch approximation spaces where u 1 , u 2 , u 3 are contravariant components of the displacement, i.e., u h = u 1 a 1 + u 2 a 2 + u 3 a 3 .Consequently, in the numerical examples included in Section 7, whenever Nitsche's method is applied, C 2 -continuous spline spaces are considered.Remark 3. As discussed, for instance, in [50], in the case of trimmed domains that present small elongated cut elements at the interface, instability effects may arise in the evaluation of the normal fluxes in the formulation (28).This problem can be easily overcome in the case in which one of the patches is not trimmed at the interface, by just selecting δ = 1 in (30), i.e., computing the flux on the non-trimmed side only.However, this is not possible when both patches are trimmed and both present small elongated cut elements at the interface.
A possible alternative is the use of stabilization techniques that have been recently proposed for the Nitsche method in the context of isogeometric discretizations [45,50,51].However, and to the best of our knowledge, no mathematically sound stabilization method has been proposed yet for isogeometric Kirchhoff-Love shells.It is our belief that the minimal stabilization technique proposed in our previous work [50] may be handy in stabilizing the problem.Nevertheless, this study is out of the scope of this paper and has not been addressed yet.In the numerical experiments included in Section 7, no numerical instabilities were found.

Local Reduced Basis method
We are interested in problems where a large amount of solution evaluations are required for problem (25) for different values of the parameters vector µ.The use of reduced order modeling techniques, such as the reduced basis method, can speedup the computation of parameterized problems.The main idea is based on an offline/online split: In the offline phase, the first step is to compute snapshots of the FOM and extract a linear combination of reduced basis functions using techniques such as the greedy algorithm or the Proper Orthogonal Decomposition method.Then, a ROM is constructed by orthogonal projection into the subspace spanned by the reduced basis.In the online phase, a reduced problem is solved to obtain the solution for any given parameter.The reader is further referred to [27,28] for more details on the reduced basis method.Nevertheless, the application of standard reduced order modeling techniques on parameterized trimmed domains entails several challenges: • The spline space V h (µ) and its dimension N h (µ) depend on the geometric parameters µ.In particular, the set of active basis functions may change for different snapshots depending on the value of the parameters µ.This impedes the construction of snapshots matrices to extract a reduced basis, since snapshots may be vectors of different length, and requires suitable snapshots extension [52].
• The efficiency of the offline/online decomposition is based on the assumption that operators depend affinely on the parameters.In case of geometrically parameterized problems, this assumption does not always hold and has to be recovered by constructing affine approximations with efficient hyper-reduction techniques [26,38].• The solution manifold obtained by extended snapshots is highly nonlinear with respect to the parameters µ.The same holds also for the affine approximations of the parameter-dependent operators.The approximation of a nonlinear manifold with a single, linear reduced basis space may yield a very high dimension of the basis.This requires tailored strategies to ensure the efficiency of the method.
In the following we will briefly review a local reduced basis method to construct efficient ROMs on parameterized trimmed domains based on clustering strategies and the Discrete Empirical Interpolation Method (DEIM).The reader is further referred to our previous work [32] for a more detailed exposition.The main steps of the local method that will be discussed in the following sections are summarized as follows: 1. perform a trivial extension of snapshots and define the extended FOM on a common, background mesh, 2. cluster the parameters and associated snapshots in order to construct DEIM and reduced basis approximations, 3. in the offline phase, train local DEIM approximations and reduced bases for each cluster combination and perform the projection, 4. choose the cluster with the smallest distance to a given parameter during the online phase and solve a reduced problem of sufficiently small dimension.

Snapshots extension
Let us first discuss the construction of snapshots.The domain at hand Ω(µ) depends on µ, and thus the support of the B-spline basis functions is also parameter-dependent.Since both the spline space V h (µ) and its dimension N h (µ) depend on the geometric parameters µ, the dimension of each solution snapshot may differ from one parameter to the other.Therefore, the first step is to perform a suitable extension of the solution vector u(µ) such that all solution vectors have the same length.As the original nontrimmed domain Ω 0 is independent of the parameters, the natural choice is to use it as a common background domain where snapshots are extended.In this work, we consider a zero extension although other choices are also possible [52].This extension is performed over the non-trimmed, multi-patch space With these definitions at hand, the extended full order problem in (25) becomes: where K ∈ R N h,0 ×N h,0 is the extended stiffness matrix, u(µ) ∈ R N h,0 is the extended solution vector and f (µ) ∈ R N h,0 is the extended right-hand side vector.Therefore, the size of the above extended problem is µ-independent and is suitable for the computation of snapshots.

Clustering strategy
In what follows we will briefly review the clustering strategy to construct local ROMs.The main idea is to build multiple approximations based on smaller subspaces instead of one global space.Then in the online phase, the closest cluster is selected for a given parameters vector µ and a local reduced problem is solved.
For geometrically parameterized problems on trimmed domains we opt for a parameter-based clustering [32].We seek for a partition of the parameter space P in N c subspaces, such that In this work, we will use the k-means clustering algorithm [53] to partition the parameter space although other partitioning strategies are also possible [33,54,55].The main idea is to assign a given parameter vector µ to the cluster that minimizes the maximum distance D(µ, k) between boundaries of the trimmed regions Ω i (µ , where ∃C > 0 and μk is the centroid of the k-th cluster1 .Then, the parameter space P is partitioned into N c subspaces as Thus, the k-means clustering minimizes the distance between each parameter vector and the cluster's centroid with respect to the Euclidean norm ∥•∥ 2 .This strategy partitions trimmed discretizations that have similar active and inactive regions. For performing such partition, we work with a discrete counterpart of the continuous space P. Thus, we create a sufficiently fine and properly selected training sample set P s = {µ 1 , . . ., µ Ns } ⊂ P of dimension N s = dim(P s ), for which the N c centroids are sought and updated iteratively until the algorithm converges.We refer the interested reader to [34,Algorithm 5] for a detailed overview of the k-means algorithm.Once the parameter space partition (34) is created, the training set P s can be clustered accordingly as Note that a suitable number of clusters N c should be chosen in advance to perform the k-means algorithm.The suitability of this choice can be evaluated a posteriori by considering the k-means variance as In particular, the k-means variance is expected to decrease with increasing number of clusters and the smallest integer can be chosen as N c at the transition between a steep slope and a plateau [56].Note that in the numerical experiments of Section 7, this transition occurs at N c ≤ 10.
Once the clustering is performed, local ROMs are constructed in the offline phase.The construction of the localized ROM will be discussed in Sections 5.4 and 5.5, while a detailed overview of the offline phase of the algorithm is given in [32, Algorithms 1-2].Afterwards, in the online phase, for a given parameters vector µ, we determine the corresponding cluster P k according to (36) and select its associated local ROM for solving a reduced problem.This online phase of the algorithm is also presented in more detail in [32, Algorithm 3].

Proper Orthogonal Decomposition
In the following we will briefly recall the Proper Orthogonal Decomposition (POD) that we will use for the construction of the ROM later on.The POD is based on the singular value decomposition algorithm (SVD) and aims to extract a set of orthonormal basis functions [57].The SVD of a matrix S ∈ R m×n reads where the orthogonal matrices have columns containing the left and right singular vectors of S, respectively, Σ ∈ R m×n is a rectangular diagonal matrix that contains the singular values σ 1 ≥ σ 2 ≥ . . .σ r , and r ≤ min(m, n) is the rank of S. The POD basis of dimension N is then defined as the set of the first N left singular vectors of S, i.e., the N largest singular values.We can choose the dimension of the basis N such that the error in the POD basis is smaller than a prescribed tolerance ε POD [28], namely N is the smallest integer such that Note that the POD basis is orthonormal by construction and the basis functions can be understood as modes that retain most of the energy of the original system.The dimension of the latter is then reduced such that the energy captured by the neglected modes in smaller than or equal to ε POD in (40).

Local Reduced Basis problem
In order to construct an efficient ROM, we seek for local reduced bases V k ∈ R N h,0 ×N k for every cluster k, where N k is the local reduced space dimension that is ideally of sufficiently small dimension, i.e., N k ≪ N h,0 .The reduced basis is constructed in the offline phase separately for each cluster based on the strategy discussed in Section 5.2.
In this work we will consider the POD to construct each reduced basis V k , while its construction was briefly discussed in Section 5.3.Note that other techniques can in principle be also chosen for the construction, as, e.g., the greedy algorithm [30].
To construct a POD basis, let us again consider a fine training sample set P s = {µ 1 , . . ., µ Ns } ⊂ P with dimension N s = dim(P s ) introduced in Section 5.2.Then we form the snapshots matrix S ∈ R N h,0 ×Ns as where the vectors u j ∈ R N h,0 denote the extended solutions u(µ j ) for j = 1, . . ., N s .These snapshots are also partitioned into N c submatrices as {S 1 u , . . ., S Nc u }, according to P s = ∪ Nc k=1 P s k in (37).The local reduced basis V k is then extracted from each cluster S k u , separately applying the POD as discussed in Section 5.3.Thus, the basis reads: Let us now derive the local reduced basis problem.For any µ ∈ P k , the solution u(µ) can be approximated using the local reduced basis V k , as where u N (µ) ∈ R N k is the solution vector of the reduced problem.A projection-based ROM can be obtained from ( 33) by enforcing the residual to be orthogonal to the subspace spanned by V k , such that Thus, the local reduced basis problem reads: while the reduced matrix and vector are defined as The size of the reduced problem ( 45) is N k ≪ N h,0 , which makes it suitable for fast online computation given many different parameters µ ∈ P k .However, the solution of the reduced problem requires the assembly of the parameter-dependent operators K(µ) and f (µ).Therefore, an important assumption for the efficiency of the reduced basis method in general is that the operators depend affinely on the parameters µ.This assumption is not always fulfilled in the presence of geometrical parameters.Therefore, we will build affine approximations to recover the affine dependence.Since we aim to approximate a manifold of extended operators that is nonlinear with respect to µ, the dimension of the approximation space may be high.Therefore, the clustering strategy of Section 5.2 will be also considered to construct local affine approximations.Note that from now on we assume for ease of exposition that the clustering is performed only once for constructing both the reduced bases and affine approximations, although in principle this could be chosen differently [32].We now introduce the following local affine approximation for any µ ∈ P k : where θ k a,q : P k → R, for q = 1, . . ., Q k a , and θ k f,q : P k → R, for q = 1, . . ., Q k f , are µ-dependent functions, whereas K k q ∈ R N h,0 ×N h,0 and f k q ∈ R N h,0 are µ-independent forms 2 .Since the latter forms do not depend on the parameters µ, they can be precomputed and stored in the offline phase.Then, the online assembly requires only the evaluation of θ k a,q , θ k f,q , which is inexpensive assuming that To obtain the affine approximation in the form of ( 47), we will employ the Discrete Empirical Interpolation Method in combination with Radial Basis Functions Interpolation.This hyper-reduction strategy will be further discussed in Section 5.5.Once the affine approximation is recovered, inserting (47) into (46) yields, for any given parameter µ ∈ P k : where K N (µ) ∈ R N k ×N k and f N (µ) ∈ R N k are the reduced matrix and right-hand side vector, respectively.We remark that in (48) only the coefficients θ k a,q , θ k f,q depend on the parameters µ and are evaluated online, while all other quantities are assembled and stored in the offline phase.Finally, during the online phase, for any given parameter µ the respective cluster is selected as in (36) and the local reduced problem in (45) is solved considering the approximation assembly in (48).Finally, the high-fidelity approximation of the solution can be recovered through (43).We remark that the efficiency of the overall method depends on the size of the local reduced problem N k , on the number of local affine terms Q k a , Q k f , as well as the efficient online evaluation of the coefficients θ k a,q , θ k f,q .The latter aspects will be further elaborated in Section 5.5.

Local Discrete Empirical Interpolation Method
In this section we will briefly present the hyper-reduction strategy based on the Discrete Empirical Interpolation Method (DEIM) for matrices and vectors.The reader is further referred to [38] for a detailed presentation of the method.
As discussed before, the first crucial step for the efficiency of the ROM involves constructing local affine approximations in the form of Equation (47).These are constructed separately for each cluster during the offline phase.For ease of exposition, 2 Henceforward, and for the sake of clarity, whenever it is clear from the context, the ranges 1, . . ., Q k a and 1, . . ., Q k f of the index q referred to the terms of the local affine approximation (47) will be omitted.
we consider the same training sample set P s = {µ 1 , . . ., µ Ns } ⊂ P of dimension N s as the one in (41), although other choices are also possible [32].Then we form the snapshots matrices S a ∈ R N 2 h,0 ×Ns and S f ∈ R N h,0 ×Ns where the vectors , with i = 1, . . ., N s , denote the vectorization of the extended stiffness matrix and the extended right-hand side vector, respectively.Following the training sample set partitioning P s = ∪ Nc k=1 P s k in (37), the snapshots matrices S a and S f are also partitioned into N c submatrices {S 1 a , . . ., S Nc a } and {S 1 f , . . ., S Nc f }, accordingly.Then, we apply the POD to each submatrix to obtain the matrices K k q and vectors f k q in (47).Here, the number of affine terms Q k a and Q k f can be determined by prescribing a tolerance ε POD as in (40).It should be remarked that the latter should in general be lower than the tolerance used to construct the local reduced basis, so that the accuracy of the DEIM approximation does not impede the overall accuracy of the ROM [28].Now, let us discuss how to efficiently compute the parameter-dependent coefficients θ k a,q (µ) and θ k f,q (µ) in (47) for each cluster.For this purpose, we will use the known as magic points [58] according to the empirical interpolation procedure [26].For each local affine approximation k, a collection of Q k a , Q k f entries is selected based on a greedy algorithm that minimizes the interpolation error over the snapshots [28].In what follows, the selected magic points are denoted as J k a , J k f for the stiffness matrix and right-hand side vector, respectively.These entries fulfill exactly the following interpolation constraints for the stiffness matrix and right-hand side vector for each µ ∈ P k : We remark that the two right-hand sides in the equations above require the online assembly of a collection of Q k a /Q k f FOM matrix/vector entries for a given µ ∈ P k , which can be costly and intrusive.Therefore, we will opt for interpolating with radial basis functions (RBFs) following our previous work [32].Similar to the construction of local affine approximations, local RBF-interpolants are constructed separately for each cluster k.The main idea is to compute offline the values of θ k a,q (µ) and θ k f,q (µ) as in (50), for each training sample µ ∈ P s k , and train a fast interpolant using these computations.Then, in the online phase the local interpolants can be evaluated rapidly for any given µ ∈ P k , being the coefficients θ k a,q (µ) in ( 47) approximated as where N k s = dim(P s k ) is the number of training samples associated to the k-th cluster, ϕ k q,j is the radial basis function associated to the j-th center parameter point µ j and ∥•∥ 2 denotes the Euclidean norm.For the numerical experiments in Section 7 we will use cubic RBFs, although other types of functions can be also chosen [59].During the offline phase, the unknown weights ω k a,q,j are computed separately for each cluster such that they fulfill the interpolation constraint exactly for The procedure is identical for the coefficients θ k f,q associated to the right-hand side vector and is therefore omitted here.
So far we have presented two approximations with respect to the FOM in ( 25), namely the local reduced problem as well as affine approximations of the extended stiffness matrix K(µ) and right-hand side vector f (µ).The construction of localized reduced bases and local affine approximations via DEIM allows to confine the dimension of the bases and number of affine terms, respectively.Moreover, the RBFinterpolation of the coefficients in ( 51) enables a rapid evaluation in the online phase.It should be noted that the localized method requires additional offline effort to construct and store multiple bases, nevertheless, the main advantage is the reduction of the online computational cost.

ROM-based parametric shape optimization
In this section, we will briefly review the aplication of the discussed reduction strategies to optimization problems, following closely the notation in [57].Such problems require several evaluations of the solution and the objective function to be minimized, which can be expensive.Therefore, these can benefit from reduced basis approximations.
First, let us assume that µ controls the shape of the computational domain at hand Ω(µ) for the optimization process.As a first step, a reduced model is constructed offline for the design variables µ following the approach presented in Section 5.Then, the optimization is performed online.In what follows we will directly use the discrete reduced approximation u N = u N (µ) of Equation ( 43) for our exposition.
The parametric optimization problem that we will consider in this paper is the compliance minimization under a given volume constraint, which is a common choice in structural optimization.The minimization of the compliance implies that the structure deforms less, i.e., it becomes stiffer.Let us now formulate the optimization problem at hand as: where V (µ) is the volume of the domain Ω(µ), V 0 a prescribed maximum volume, and µ min , µ max the lower and upper bounds for the design variable µ, respectively.Here, the cost functional J N (u N , µ) is obtained by evaluating the reduced basis approximation of the problem solution.For the compliance case, the reduced objective function reads Remark 4. In the case of the extended FOM problem, the compliance functional (54) can be written as that, after introducing the approximations (43) and ( 46), becomes where the orthogonality of the basis V k was considered.
There are several ways to solve the optimization problem ( 53)-( 54).In the following we will consider a gradient-based approach, where the parameters µ are updated in an iterative fashion depending on the gradient of the cost functional J N .The gradient can be evaluated either analytically or based on a suitable approximation, e.g., with a finite difference scheme.The latter yields a black-box optimization approach that simply requires the reduced solution and the evaluation of the objective function.Nevertheless, the reduced model is perfectly suitable for computing parametric sensitivities due to its differentiability and affine parametric dependence [28,Proposition 5.3].Let us now reformulate the parameterized reduced problem as With these definitions at hand and following [28,Proposition 11.3], the gradient of the objective function JN (µ) = J N (u N , µ) reads: The evaluation of these gradient requires the solution of M sensitivity equations where D denotes the Fréchet derivative.Instead of directly solving the above equations at each step of the optimization process, we define an additional adjoint problem and denote its solution as It should be remarked that this approach requires in general the construction of a reduced model for the adjoint problem, which implies additional offline cost.This can be constructed following the approach in Section 5. Given the reduced adjoint approximation u N , the gradient in Equation ( 58) can be reformulated as: Note that the evaluation of the above derivatives with respect to µ is inexpensive assuming that the operators depend affinely on the parameters.Considering the compliance case, the solution of the adjoint problem can be directly obtained from the reduced solution u N .Thus, by inserting Equations ( 54) and (57) into Equation ( 60), the adjoint problem for the compliance case reads: Finally, inserting Equation ( 62) into (61), and expanding the derivative D µ G N (u N , µ), the gradient reads: Under the assumption of affine parametric dependence in (47) and considering ( 34) and ( 43), we rewrite the gradient for µ ∈ P k , as The derivatives ∂θ k a,q (µ)/∂µ are simple and inexpensive to evaluate by just differentiating the expression (51).The procedure for ∂θ k f,q (µ)/∂µ is identical.

Numerical results
In this section we present some numerical experiments for the Kirchhoff-Love shell problem to assess the capabilities of the presented ROM framework for parameterized trimmed and multi-patch geometries.In what follows we apply the two presented strategies to enforce interface coupling conditions in a weak sense, namely the projected super-penalty (Section 4.1) and Nitsche's method (Section 4.2).The numerical experiments are carried out using the open-source Octave/Matlab isogeometric package GeoPDEs [60] in combination with the re-parameterization tool for integration of trimmed geometries presented in our previous works [43,47], while the ROM construction exploits the open-source library redbKIT [61].We remark that the stiffness matrix is preconditioned with a diagonal scaling to avoid large conditions numbers and loss of accuracy due to small trimmed elements as discussed in [62].

Scordelis-Lo roof with holes
The first numerical example is a single-patch, trimmed variant of the Scordelis-Lo roof.The geometry and material properties are adopted from the well known benchmark, see, e.g., [63] for more details.Thus, the Young's modulus is E = 4.32 • 10 8 Pa, the Poisson's ratio ν = 0, and the thickness t = 0.25 m.The shell structure is subjected to self-weight with a vertical loading of f z = −90 N/m 2 as depicted in Fig. 2. Rigid diaphragm boundary conditions are imposed on the curved ends of the shell, that is, we fix the displacements in the xz-plane.In this example, we cut out two circular holes in the parametric domain as shown in Fig. 2. The radius of the holes in the parametric domain is fixed as r = 0.2.Their location is parameterized, where µ ∈ P = [0, 0.1] is a geometric parameter representing the location of the center of each circle that moves along the diagonal of the unit square.The geometry of the shell structure in the physical space is then obtained by an additional mapping as mentioned in Section 2.
The shell is discretized with cubic C 2 -continuous B-splines and the dimension of the non-trimmed space is N h,0 = 1083.Let us now construct a ROM with the local reduced basis method presented in Section 5.For this purpose, we use a training sample of dimension N s = 500 obtained by Latin Hypercube sampling [64] for the construction of the reduced basis.Note that this sampling is also suitable for more general problems with higher dimension of the parameter space [28].Fig. 3a depicts the decay of the singular values of the POD for different numbers of clusters.We remark that the POD basis is constructed such that it minimizes the squared projection error with respect to the matrix norm, that represents the algebraic counterpart of the H 2 norm.Moreover, the number of clusters for the DEIM approximations is fixed to 16 for all computations.It can be observed that the decay is more rapid with localized ROMs.Now we employ a test sample of dimension N t = 30 with a uniform random distribution to perform the error analysis whose results are presented in Fig. 3b.We observe that increasing the number of clusters further reduces the dimension of the reduced basis, while an accuracy of 10 −5 is achieved in the H 2 norm.The optimal number of clusters is selected as N c = 8 based on the k-means variance in Fig. 4. Note that the computation of the optimal number of clusters is performed offline and requires 0.21 s for the problem at hand.Fig. 5 depicts the vertical displacement solutions, where the local ROM with N c = 8 clusters is compared to the FOM.The online CPU time is for the ROM is 0.12 s including the closest cluster iterations, whereas the assembly and solution with the FOM requires 2.32 s.This results in a speedup of 19×.Note that the trimming operation for computing the FOM or vizualizing the results requires additionally 0.854 s.We remark that the offline time to compute the snapshots for the POD basis is 3 minutes exploiting the affine approximation (47).The offline time to compute the snapshot matrices and vectors for the DEIM approximations is 1.75 hours given a training sample of 2000 snapshots.We remark that this includes the trimming operation for each snapshot and can be further optimized with the use of parallelization.Moreover, we perform optimization using both the local ROM and the FOM.The upper and lower bounds for the design variable are defined by the parameterization at hand, that is we seek the optimal shape within the bounds 0 ≤ µ ≤ 0.1 during the optimization.For this problem we do not consider any volume constraints.First, we solve the optimization problem with the FOM using a finite difference scheme to compute sensitivities.The solution requires 7 iterations and 23 function evaluations in total.Then, we employ the local ROM with N c = 8 clusters for the optimization.The ROM with approximate sensitivities requires 17 function evaluations versus the ROM with exact sensitivities with 15 function evaluations.Note that a forward finite difference scheme is employed for the computation of the approximate sensitivities during the optimization.The solution is slightly faster for the ROM with exact sensitivities and comprises 76.1 ms versus 79.5 ms for the solution with finite difference approximation of the gradient.In all cases the optimal shape is obtained for µ = 0.0319 after 7 iterations.Fig. 6 depicts the optimization results.The optimization history is shown in Fig. 6a by depicting the evolution of the relative compliance during the optimization.For the optimal solution, the compliance is reduced 8% respect to the initial configuration (i.e., for µ = 0).Finally, Fig. 6b shows the vertical displacement solution for the final optimal shape, comparing ROM and FOM solutions.

Multi-patch simple geometries
In this section we will assess the capabilities of the ROM framework for multi-patch geometries.For this purpose we employ two simple geometries, namely the multipatch Scordelis-Lo roof and two non-conforming planar patches.The coupling of patches is achieved in both cases using the projected super-penalty method discussed in Section 4.1.

Non-trimmed multi-patch Scordelis-Lo roof
This example is intended to test the ROM framework on a non-trimmed, multi-patch setting.For this purpose, we employ again the Scordelis-Lo roof and split the geometry into two subdomains.The geometric setup is depicted in Fig. 7.We employ the projected super-penalty method to enforce interface coupling conditions.The material parameters, loading, and boundary conditions are the same as in the previous example.The shell structure is modeled with two conforming subdomains and the parameterization of the multi-patch design is shown in Figs.7b-7c for a coarse geometry.The common interface is depicted in red color.Note that the geometry is parameterized by moving the depicted control points in the vertical direction, where the geometric parameter µ ∈ P = [0, 10] defines their position in the y-direction prescribing the curvature of the shell structure.We remark the µ = 0 corresponds to the original coordinates of the Scordelis-Lo roof benchmark.In Figs.7b-7c we denote for simplicity P (µ) = P (x, y + µ, z).The parameters are first defined on a coarse mesh and then the geometry is further refined for the analysis as discussed in [30].The geometry is discretized with quadratic C 1 -continuous B-splines for the analysis resulting in N h = 1944 degrees of freedom.Now let us construct a ROM for the multi-patch geometry.For this purpose, we employ a training sample of dimension N s = 500 obtained by Latin Hypercube sampling.First, the affine approximations are constructed with the DEIM.Fig. 8 depicts the error of the DEIM approximations in the L ∞ norm for both the right-hand side vector and the stiffness matrix.It can be observed that the error already decays rapidly with one global approximation.The error reaches an accuracy of 10 −7 with Q a = 8 and Q f = 5 affine terms for the stiffness matrix and the right-hand side vector, respectively.As a further step, a reduced basis is constructed with the POD.Fig. 9 shows the decay of the singular values and the relative error of the reduced solution in the H 2 norm.The error analysis is performed using a test sample of dimensions N t = 30 and a uniform random distribution.Similarly to the DEIM approximations, the decay is already rapid by constructing only one global reduced basis space and the error reaches an accuracy of 10 −7 for a reduced basis space of dimension N = 7.In Fig. 10 we compare the vertical displacement solution of the ROM with the solution of the full order model for three different parameter values.
Moreover, we perform optimization considering a volume constraint.The initial volume is set as V 0 = 1.9 • 10 3 and we seek the optimal shape for the design variable bounded by 0 ≤ µ ≤ 10.The optimization is completed after 4 iterations and 14 function evaluations.The ROM-based optimization is computed in 67 ms, while the optimization with the high fidelity solution is solved in 518 ms.This implies a speedup of 7.76×.In both cases the exact gradients are computed at every optimization step.At the final iteration, the compliance is reduced by 16.45% compared to the initial configuration.Fig. 11 depicts the vertical displacement solutions for the optimal shape obtained for µ = 5.3127 with both the ROM and full order solution, while Table 1 summarizes the main results and obtained computation times.We conclude that this verifies the suitability of the ROM framework for multi-patch geometries coupled along non-trimmed interfaces.

Trimmed non-conforming planar patches
In this example we aim to assess the capabilities of the ROM for trimmed multi-patch geometries with both conforming and non-conforming discretizations.In particular, we investigate the behavior of the local ROM for patches that are coupled along parameterized trimming interfaces and are expected to behave poorly with a standard global reduced basis.For this purpose we consider a planar setting, where the computational domain is a unit square and is subdivided into two patches coupled along a curved 100 .Note that since both patches are trimmed, the internal knots of the coupling curve are neglected for the construction of the interface knot vector Ξ j following the discussion in Remark 1.The geometry is discretized with quadratic C 1 -continuous B-splines and the dimension of the non-trimmed space is N h,0 = 1944.The applied boundary conditions and loading are adopted from a manufactured smooth solution given in [13] such that (u x , u y , u z ) = (0, 0, sin(πx) sin(πy)).Now let us construct local ROMs for both the conforming and non-conforming case.For this purpose, we employ a training sample of dimension N s = 500 that we obtain by Latin Hypercube sampling.First, we investigate the k-means variance in Equation (38) in order to choose the number of clusters.As it can be seen in Fig. 13, the k-means variance does not decrease significantly after N c = 10 clusters.To perform the error analysis, we consider a testing sample of dimension N t = 30, which is obtained by a uniform random distribution.The relative error in H 2 norm is depicted in Fig. 14 for N c = 8, 10 and the results of both discretizations are compared to each other.As expected, the size of the reduced basis decreases for increasing number of clusters.We observe that the non-conformity slightly affects the maximum number of reduced basis functions N , while the error is of the same magnitude in all cases.It can be concluded that the ROM framework is suitable for multi-patch geometries coupled at trimmed interfaces with both conforming and non-conforming discretizations.

Joint of intersecting tubes
The last example aims to demonstrate the capabilities of the presented framework for complex geometries.For this purpose we consider the geometry of three intersecting tubes (see Fig. 15) that represent a generic configuration for, e.g., joints in steel support truss structures (see also [14,17] for a similar variant).We remark that optimizing such large-scale structures of industrial relevance is a challenging task, where the shape parameters of each joint element may differ and the ROM can be reused for several online evaluations.The optimization of such problems requires further measures for efficiency (e.g.domain decomposition strategies), which is out of the scope of the present work.For this particular geometry, both sides are trimmed at all interfaces, which hinders the use of the projected super-penalty approach.To this end, we employ the Nitsche's method to enforce interface coupling conditions.The material parameters are the Young's modulus E = 3 • 10 6 , the Poisson's ratio ν = 0.3, and the thickness t = 0.2.The radius of the main tube is R 1 = 1.Homogeneous Dirichlet boundary conditions are applied on the top and bottom sections of the main tube, while periodic boundary conditions are applied at the cylinders' closure.A vertical load f z = 10 is applied on all tubes.The geometry is discretized with cubic C 2continuous B-splines for the analysis and results in N h,0 = 2673 degrees of freedom.In the following we will demonstrate the capabilities of the ROM for two test cases of geometrical parameterization.

Parameterization of radius
First, we consider a geometrical parameter µ ∈ P = [0.6,0.8] that represents the radius of the connecting skewed tubes.The geometry setup and parameterization is depicted in Fig. 15.Note that the angles of the connecting tubes with respect to the horizontal plane are fixed to 45 • and −30 • for the top and bottom tubes, respectively (see Fig. 16).In Fig. 17, the trimming configuration, including the trimmed parametric domains and quadrature points are shown for a particular case of µ.To construct the ROM, we employ a training sample of dimension N s = 1000 that we obtain by Latin Hypercube sampling.Note that this refers to the global dimension before clustering, while the snapshots are computed in parallel to speedup the offline phase.The k-means variance is depicted in Fig. 18a for increasing number of clusters.We observe that the variance does not change significantly after 10 clusters, thus N c = 10 is chosen in what follows.The error analysis is performed using a test sample of dimension N t = 30 obtained by a uniform random distribution.Fig. 18b shows the error convergence for N c = 10 while in Fig. 19 the displacement solution obtained with the ROM is compared to the FOM for three parameter values from the test sample.The main results and computation times are summarized in Table 2.

Parameterization of angle
Let us now consider a geometrical parameter that represents the angle of the connecting tubes with respect to the horizontal plane.In this case, the radius of the intersecting tubes is fixed as R 2 = 0.8 and the absolute value of the angle varies as µ ∈ P = [25.7 • , 36 • ].Note that here the angles of both tubes are symmetric with respect to the horizontal plane.The geometrical variation of the angle is illustrated in     clusters is chosen as N c = 10 based on the k-means variance in Fig. 20b.Fig. 21 shows the absolute and relative error in H 2 norm with respect to the maximum number of reduced functions over all clusters, while Fig. 22 compares the displacement solutions using the ROM with the FOM for three parameter values of the test sample.Moreover, we employ the ROM to solve an optimization problem that minimizes the compliance within the parameter space P without volume constraints.Here, we consider a displacement constraint such that max(u(µ)) ≤ 2.5 • 10 −3 .That is, the compliance is minimized such that the maximum displacement does not exceed the prescribed value.
We remark that we also employ the ROM for the computation of the displacement constraint to speedup the optimization.Note that the initial shape at the beginning of the optimization corresponds to the minimum angle with µ = 25.7 • .The optimization using exact sensitivities requires 8 iterations and 29 function evaluations.Fig. 23a depicts the evolution of the relative compliance during the optimization.At the final iteration, the compliance is decreased by 17.9% compared to the initial configuration.Note that the optimization with approximate sensitivities using a forward finite difference scheme requires 8 iterations and 21 function evaluations to reach almost identical results.Fig. 23b depicts the displacement solution for the optimal shape with µ = 32.22 • and Table 3 summarizes the main results and computation times.It is remarked that the online computation time corresponds to the cluster with the maximum number of basis functions, thus the online cost might differ from one parameter to the other.

Conclusion
In this work we present a parametric model reduction framework for trimmed, multipatch isogeometric Kirchhoff-Love shells.The proposed strategy is suitable for fast simulations on parameterized geometries represented by multiple, non-conforming patches where the trimming interfaces change for different values of the geometric parameters.Following our previous work [30], efficient ROMs are constructed using a local reduced basis method and hyper-reduction with DEIM.The latter enables an efficient offline/online split with low online cost, which is advantageous for solving parametric optimization problems.
We have investigated the capabilities of the proposed framework through several numerical experiments.To this end, we considered trimmed, multi-patch geometries with parameterized interfaces, considering both conforming and non-conforming cases, and that were glued applying super-penalty and Nitsche coupling methods.We observed a high accuracy of the local ROMs, while the solution evaluation in the online phase was obtained with low computational cost.Moreover, we validated the proposed approach for parametric optimization problems and applied it to a complex geometry.
Overall, the application of the ROM framework to isogeometric Kirchhoff-Love shell analysis of complex geometries and optimization problems is a cost-effective alternative.The application to more complex optimization problems, including higher number of design parameters, and the extension to Reissner-Mindlin shell formulations are future research directions to explore.Regarding the ROM, a further subject of future work is related to error certification using greedy algorithms and tailored a posteriori error estimators.

Fig. 3 :
Fig. 3: Example 7.1: Singular values decay and relative error in H 2 norm vs. maximum number of reduced basis functions N over all the clusters, for different numbers of clusters.
(a) Evolution of the compliance (b) Final optimal shape

Fig. 6 :
Fig. 6: Example 7.1: Optimization results depicting: (a) the evolution of the relative compliance during the optimization comparing the FOM, the ROM with approximate sensitivities, and the ROM with exact sensitivities; and (b) the vertical displacement for the final optimal shape with µ = 0.0319 using the FOM (top) and the local ROM (bottom).

Fig. 8 :
Fig. 8: Example 7.2.1:Error decay of DEIM approximations in L ∞ -error norm for right-hand side vector and stiffness matrix using a POD tolerance of ε POD = 10 −7 .

Fig. 9 :
Fig. 9: Example 7.2.1:Singular values decay and relative error in H 2 norm vs. number of reduced basis functions N.

Fig. 12 :
Fig. 12: Example 7.2.2:Geometry setup and parameterization of the two planar trimmed patches for different values of µ.The trimming interface is a quadratic spline curve and denoted in red color.

8 Fig. 15 :
Fig. 15: Example 7.3.1:Geometrical parameterization for different values of µ (radius of skewed cylinders) for the joint of intersecting tubes.

Fig. 16 :
Fig. 16: Example 7.3.1:Exemplary side view for the joint of intersecting tubes with angle configuration.

Fig. 17 :
Fig. 17: Example 7.3.1:Trimming configuration for µ = 0.7465.On the left, trimming patches in the physical domain, on the center and right, trimmed parametric domains (trimmed elements are shaded).Black dots depict quadrature points in the cut elements, crosses (in physical domain) quadrature points at interfaces.

Fig. 20a .
Fig. 20a.Training and test samples are generated similarly to the previous test case for the construction of the ROM and the error analysis, respectively.The number of

Fig. 18 :
Fig. 18: Example 7.3.1:K-means variance over number of clusters N c and relative error in H 2 norm vs. maximum number of reduced basis functions N over all clusters.

Fig. 20 :
Fig. 20: Example 7.3.2:Geometrical parameterization and k-means variance over number of clusters N c .

Fig. 21 :
Fig. 21: Example 7.3.2:Absolute and relative errors in H 2 norm vs. maximum number of reduced basis functions N over all clusters.

Fig. 23 :
Fig. 23: Example 7.3.2:Optimization results: (a) the evolution of the relative compliance during the optimization for the ROM with exact sensitivities; and (b) the displacement for optimal shape with µ = 32.22 • using the ROM.

Table 1 :
Example 7.2.1:Number of basis functions and computational cost.

Table 2 :
Example 7.3.1:Number of basis functions and computational cost.

Table 3 :
Example 7.3.2:Number of basis functions and computational cost.