T Hirschler
email: thibaut.hirschler@utbm.fr

R Bouclier

P Antolin

A Buffa

Reduced Order Modeling based Inexact FETI-DP solver for lattice structures

Keywords: Multiscale Mechanics, Domain decomposition, Isogeometric analysis, Additive Manufacturing, Reduced basis, Architected materials

This paper addresses the overwhelming computational resources needed with standard numerical approaches to simulate architected materials. Those multiscale heterogeneous lattice structures gain intensive interest in conjunction with the improvement of additive manufacturing as they offer, among many others, excellent stiffness-to-weight ratios. We develop here a dedicated HPC solver that benefits from the specific nature of the underlying problem in order to drastically reduce the computational costs (memory and time) for the full fine-scale analysis of lattice structures. Our purpose is to take advantage of the natural domain decomposition into cells and, even more importantly, of the geometrical and mechanical similarities among cells. Our solver consists in a so-called inexact FETI-DP method where the local, cell-wise operators and solutions are approximated with reduced order modeling techniques. Instead of considering independently every cell, we end up with only few principal local problems to solve and make use of the corresponding principal cell-wise operators to approximate all the others. It results in a scalable algorithm that saves numerous local factorizations. Our solver is applied for the isogeometric analysis of lattices built by spline composition, which offers the opportunity to compute the reduced basis with macro-scale data, thereby making our method also multiscale and matrix-free. The solver is tested against various 2D and 3D analyses. It shows major gains with respect to black-box solvers; in particular, problems of several millions of degrees of freedom can be solved with a simple computer within few minutes.

Introduction

Additive Manufacturing (AM) and especially its metal variants constitute today a reality for the fabrication of industriallyrelevant high-performance parts and products [START_REF] Blakey-Milner | Metal additive manufacturing in aerospace: A review[END_REF]. In this context, AM has pushed the development of architected materials to a new stage. These can be viewed as specific multiscale heterogeneous lattice structures that are inspired, in particular, from biology such as bone architectures. In practice, they usually consist in repeating a pattern, the unit-cell, made of a combination of struts connected in different ways, all over a macro shape (see Figure 1). The cell topology, geometry and material can be tailored for achieving specific performance requirements [START_REF] Pan | Design and Optimization of Lattice Structures: A Review[END_REF]. Starting with unprecedented weight saving while maintaining the stiffness and strength of the structure, lattices can also accommodate novel behaviors (e.g., highly stretchable and auxetic) and multi-functionality, which make them attractive for lightweight components, energy absorbing, thermal management, design of medical implants, to name a few [START_REF] Benedetti | Architected cellular materials: A review on their mechanical properties towards fatigue-tolerant design and fabrication[END_REF][START_REF] Schaedler | Architected cellular materials[END_REF].

With the improvement of AM, representative lattices may now come with large numbers of cells of complex topology and shape, which make their numerical simulation very challenging. In particular, the computational cost and memory requirements become quickly tremendous or even intractable if standard methods are blindly applied. Thus, the efficient prediction of the structural response of lattice structures is becoming a hot topic nowadays.

Full lattice structure

Design space

Reference unit-cell

Figure 1: Example of a lattice structure (inspired from [START_REF]Element pro case study: F1 brake pedal[END_REF]): The pedal crank is composed of a lattice structure (right) which allows drastic weight saving compared to a classic design (left).

As a solution, one may be tended to draw inspiration from the multitude multiscale methods based on homogenization that have been developed over the years for multiscale heterogeneous structures. The goal of these approaches is to compute a homogenized macro-scale behavior from its underlying heterogeneous microstructure by different interchange of information between scales. Among others, let us mention the multilevel FEM (FE 2 [START_REF] Benaimeche | A k-means clustering machine learning-based multiscale method for anelastic heterogeneous structures with internal variables[END_REF][START_REF] Feyel | A multilevel finite element method (FE2) to describe the response of highly non-linear structures using generalized continua[END_REF][START_REF] Praster | An adaptive FE2 approach for fiber-matrix composites[END_REF][START_REF] Xu | Data-driven multiscale finite element method: From concurrence to separation[END_REF]) which computes homogenized fields from the micro-scale and assigns them at each integration point of the macro-scale model; the Multiscale FEM (MsFEM [START_REF] Hou | A Multiscale Finite Element Method for Elliptic Problems in Composite Materials and Porous Media[END_REF][START_REF] Yan | Multiscale analysis for 3D lattice structures based on parallel computing[END_REF]) where numerically computed basis functions encode the micro-scale heterogeneity; global/local coupling [START_REF] Wangermez | Interface coupling method for the global-local analysis of heterogeneous models: A second-order homogenization-based strategy[END_REF] that allows to replace the homogenized global model by the fine high-fidelity model in a certain local region; or direct numerical homogenization [START_REF] Grimm-Strele | FFT-based homogenization with mixed uniform boundary conditions[END_REF][START_REF] Khakalo | Modelling size-dependent bending, buckling and vibrations of 2D triangular lattices by strain gradient elasticity models: applications to sandwich beams and auxetics[END_REF][START_REF] Moulinec | A FFT-Based Numerical Method for Computing the Mechanical Properties of Composites from Images of their Microstructures[END_REF][START_REF] Torabi | Application of a non-conforming tetrahedral element in the context of the three-dimensional strain gradient elasticity[END_REF] where macro-scale material parameters are identified. However, these methods mostly rely on a true separation of scales to alleviate size effects [START_REF] Onck | Size effects in ductile cellular solids. Part I: modeling[END_REF][START_REF] Yoder | Size effects in lattice structures and a comparison to micropolar elasticity[END_REF]. This is usually not the case with lattice structures since 3D printers notwithstanding constraint the achievable length scale range of the cells. Furthermore, the macro-shape often follows a rather slender geometry (plate, shell), thereby exhibiting few cells in one direction. Alternatively, specific structural models such as sophisticated 3D beams [START_REF] Weeger | Isogeometric sizing and shape optimization of 3D beams and lattice structures at large deformations[END_REF][START_REF] Weeger | Fully isogeometric modeling and analysis of nonlinear 3D beams with spatially varying geometric and material parameters[END_REF] could be used to model the full architecture of the lattice while reducing the complexity of the problem. Yet, such methods may not be suitable in case of small cells (thus involving thick struts) and cannot be applied to more general lattices, (e.g., skeletal-TPSM or sheet-TPMS based lattices [START_REF] Benedetti | Architected cellular materials: A review on their mechanical properties towards fatigue-tolerant design and fabrication[END_REF]). As a result, it seems that instead of acting at the modeling stage, it would be desirable to work on the solution stage for lattices; that is, to derive HPC solvers capable of solving the high-fidelity, fine-scale problem efficiently.

In this spirit, immersed domain methods, such as the Finite Cell Method [START_REF] Parvizian | Finite cell method[END_REF][START_REF] Schillinger | The Finite Cell Method: A review in the context of higher-order structural analysis of cad and image-based geometric models[END_REF], have recently been applied for the full fine-scale analysis of lattice structures [START_REF] Korshunova | Image-based material characterization of complex microarchitectured additively manufactured structures[END_REF][START_REF] Korshunova | Bending behavior of octet-truss lattice structures: Modelling options, numerical characterization and experimental validation[END_REF][START_REF] Korshunova | Uncertainty quantification of microstructure variability and mechanical behavior of additively manufactured lattice structures[END_REF]. Such approaches switch the effort of mesh generation to the proper integration of cut elements and the sophisticated treatment of conditioning problems. Therefore, when used in combination with image-based modeling [START_REF] Düster | Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method[END_REF][START_REF] Rouwane | Adjusting fictitious domain parameters for fairly priced imagebased modeling: Application to the regularization of digital image correlation[END_REF][START_REF] Verhoosel | Image-based goal-oriented adaptive isogeometric analysis with application to the micro-mechanical modeling of trabecular bone[END_REF], they make it easier to take into account the (strut-level) geometrical defects due to the AM process (see [START_REF] Korshunova | Image-based material characterization of complex microarchitectured additively manufactured structures[END_REF][START_REF] Korshunova | Bending behavior of octet-truss lattice structures: Modelling options, numerical characterization and experimental validation[END_REF][START_REF] Korshunova | Uncertainty quantification of microstructure variability and mechanical behavior of additively manufactured lattice structures[END_REF] again). HPC extensions of such technologies exist [START_REF] Badia | Robust and scalable domain decomposition solvers for unfitted finite element methods[END_REF][START_REF] Jomo | Robust and parallel scalable iterative solutions for large-scale finite cell analyses[END_REF]. They are generic in terms of applications but they may not appear optimal in the specific context of lattice structures. In particular, they do not make use of the geometric similarity among different cells in the numerical solution, what allows a drastic reduction in terms of memory and computational cost, as demonstrated for the assembly of the numerical operators in [START_REF] Hirschler | Fast and multiscale formation of isogeometric matrices of microstructured geometric models[END_REF]. The objective of this work is to bring forward this idea of benefiting from the repetitive character of the lattice geometries to build a dedicated HPC solver.

In order to do so, our approach relies on three main ingredients. Firstly, it is based on the Dual-Primal variant of the Finite Element Tearing and Interconnecting method (FETI-DP), that, after its appearance in [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method-part I: A faster alternative to the two-level FETI method[END_REF] has been widely developed and adopted in high-performance computational structural mechanics (see [START_REF] Klawonn | A parallel implementation of dual-primal FETI methods for three-dimensional linear elasticity using a transformation of basis[END_REF][START_REF] Klawonn | Parallel adaptive FETI-DP using lightweight asynchronous dynamic load balancing[END_REF][START_REF] Toivanen | A multilevel FETI-DP method and its performance for problems with billions of degrees of freedom[END_REF][START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF] to name a few). This algorithm belongs to the class of non-overlapping Domain Decomposition (DD) solvers that appear highly relevant here given the natural splitting of the lattice structures into several subdomains (each cell can be associated to a subdomain). DD methods are preconditioned iterative algorithms where, at each iteration, subdomain-independent local systems are solved in parallel. The interface conditions are recovered through the iterative process. To obtain a scalable algorithm, a small coarse problem also needs to be solved at each iteration. We note that DD solvers have been used in the past to successfully compute multiscale heterogeneous structures (see, e.g., [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF][START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF][START_REF] Lloberas-Valls | On micro-to-macro connections in domain decomposition multiscale methods[END_REF]). In FETI-DP, the continuity constraints on the displacement at the subdomain corners are forced to hold throughout the iterative procedure, which naturally leads to solving a coarse problem at each iteration, while the other constraints are enforced by the use of Lagrange multipliers. The second ingredient of our approach is to consider the framework of inexact FETI-DP algorithms [START_REF] Klawonn | Inexact FETI-DP methods[END_REF][START_REF] Rheinbach | Parallel iterative substructuring in structural mechanics[END_REF]. In particular, we develop an inexact version that offers the opportunity to avoid solving numerous local (subdomain-wise) systems. For this, we iterate on the initial complete saddle point problem and design a block preconditioner consistent with the mathematical background established in [START_REF] Notay | A new analysis of block preconditioners for saddle point problems[END_REF][START_REF] Simoncini | Block triangular preconditioners for symmetric saddle-point problems[END_REF]. Thirdly, we investigate the repetitiveness of the cells by applying a Reduced Order Modeling (ROM) approach [START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF][START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF] (especially, a greedy technique) that extracts the "principal" cells in terms of stiffness. The resulting few local stiffness operators are then used in the preconditioner of the inexact FETI-DP to build reduced bases for the efficient approximate solutions of the numerous local systems. It results in a scalable algorithm that saves numerous local factorizations and thus offers the opportunity to deeply reduce the computational time with respect to the standard FETI-DP solver, in particular when the macro mapping is simple which is usually the case in real lattice structures. In this contribution, we place ourselves at the design stage of the product; that is, we only consider perfect lattices. Yet, we believe that the proposed framework owns all the necessary ingredients to be easily extended to treat lattices with geometrical imperfections due to AM. In addition, we assume here, as a starting point, that the lattice is linear elastic (at the scale of the cell-struts) and we place ourselves in the context of small displacements and rotations.

Furthermore, although not mandatory in our methodology, we adopt the Computer-Aided Design (CAD) paradigm based on spline composition [START_REF] Elber | Precise Construction of Micro-structures and Porous Geometry via Functional Composition[END_REF] to generate the lattice geometries. Such a geometrical modeling offers great flexibility to design lattices since the local (cell) and global (macro) geometries are naturally parameterized by splines. In addition, it is fully consistent with IsoGeometric Analysis (IGA) [START_REF] Cottrell | Isogeometric Analysis[END_REF][START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF], thus allowing to directly analyze such models (see, e.g., [START_REF] Antolin | Optimizing Micro-Tiles in Micro-Structures as a Design Paradigm[END_REF][START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF]) while avoiding possibly delicate meshing procedures. A first step toward the efficient analysis of lattices modeled by spline composition has been recently performed with the development of a fast multiscale assembly procedure for the IGA operators (see previous contribution [START_REF] Hirschler | Fast and multiscale formation of isogeometric matrices of microstructured geometric models[END_REF]). The strategy makes use of look-up tables with precomputed integrals associated to the cell pattern and macro-fields that encode the mechanical behavior related to the macro-shape. Benefiting from this data structure, we are actually able to extract the "principal" cells in our solver by acting only on the macro-fields, thereby avoiding forming all the local stiffness operators. In this context, our ROM-based inexact FETI-DP solver thus incorporates additional attractive features: It is multiscale and tends to be matrix-free. From this point of view, the present work also further advances the state-of-the-art of DD-IGA that represents an active research area (see [START_REF] Beirão Da Veiga | BDDC preconditioners for isogeometric analysis[END_REF][START_REF] Kleiss | IETI-isogeometric tearing and interconnecting[END_REF] for the origins and then [START_REF] Bosy | A domain decomposition method for Isogeometric multi-patch problems with inexact local solvers[END_REF][START_REF] Hirschler | A dual domain decomposition algorithm for the analysis of non-conforming isogeometric kirchhoff-love shells[END_REF][START_REF] Hofer | Inexact dual-primal isogeometric tearing and interconnecting methods[END_REF][START_REF] Widlund | BDDC deluxe algorithms for two-dimensional H(curl) isogeometric analysis[END_REF] to name a few).

The remainder of this paper is organized as follows. In Section 2, we briefly review the FETI-DP method starting from the complete saddle point formulation and give the key ingredients to build efficient block preconditioners which is the basis of our solver. Then, we describe in Section 3 the developed ROM-based inexact strategy to benefit from the repetitive character of lattice geometries in the FETI-DP solver. Making use of algebraic equations, the presentation is kept general at this stage to highlight the applicability of our solver to any discretization scheme. In Section 4, we specify the application of the methodology in case of a CAD modeling based on spline composition, which results in an isogeometric solver that is also multiscale and quasi matrix-free. Next, numerical results to assess the efficiency of our solver both in terms of memory and computational cost reduction are reported in Section 5 and, finally, Section 6 contains some conclusions and future directions of research.

FETI-DP and saddle point problems

Domain decomposition

A lattice structure is by construction made of multiple cells that are packed together and fill a particular area of a macro-geometry (recall Figure 1). In other words, the domain Ω describing a lattice structure can be naturally seen as a non-overlapping collection of subdomains Ω (𝑠) , each of them being associated to a lattice unit-cell. Thus, we define: and where 𝑁 is the total number of cells. Let us also introduce the subdomain interfaces:

Ω = 𝑁 ⋃ 𝑠=1 Ω(𝑠) , with Ω (𝑠 1 ) ∩ Ω (𝑠 2 ) = ∅ if 𝑠 1 ≠ 𝑠 2 , (1)
Ω (s 4 ) Ω (s 5 ) G (s 1 ) G (s 2 ) G (s 3 ) Ω ref Ω Γ N Γ D
Γ = ⋃ 𝑠 1 ≠𝑠 2 𝜕Ω (𝑠 1 ) ∩ 𝜕Ω (𝑠 2 ) . (2)
As mentioned in the introduction, there exist specific solvers which are particularly meant to take advantage of such a domain decomposition: the so-called non-overlapping DD algorithms. The basis of this work is thus to resort to the DD technology to build an efficient solver for lattice structures.

Remark 1. Note that in the general field of DD solvers, obtaining a domain decomposition of the form (1) actually requires a pre-processing step that can be rather complex and may involve sophisticated tools (see, e.g., [START_REF] Coupez | Parallel meshing and remeshing[END_REF][START_REF] Have | Algebraic domain decomposition methods for highly heterogeneous problems[END_REF][START_REF] Wang | A parallel approach for the generation of unstructured meshes with billions of elements on distributed-memory supercomputers[END_REF]. However, here, the geometrical structure of the domain leads in a straightforward manner to a non-overlapping decomposition, making the application of DD solvers a natural choice.

Dual-primal FETI method

More precisely, we investigate in this work a class of DD solvers known as FETI-DP [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method-part I: A faster alternative to the two-level FETI method[END_REF], where the starting point is the construction of a particular working space in which we will iteratively look for our (finite element) solution. Then, a specific algorithm is derived which takes advantage of the particular structure of this working space.

Model problem

We consider here the linear elasticity problem in Ω which is naturally decomposed in 𝑁 subdomains as expressed in Equation (1) (see Figure 2 for illustration). In each subdomain, the displacement field 𝑢 (𝑠) ∶ Ω(𝑠) → ℝ 𝑑 , 𝑑 being the dimension of the problem, satisfies the following equations:

div(σ (𝑠) ) + 𝑏 ext = 0 in Ω (𝑠) , (3a)
σ (𝑠) = C (𝑠) ε(𝑢 (𝑠) ) in Ω (𝑠) , (3b)
𝑢 (𝑠) = 𝑢 0 over Γ D ∩ 𝜕Ω (𝑠) , (3c)
σ (𝑠) 𝑛 (𝑠) = 𝑡 ext over Γ N ∩ 𝜕Ω (𝑠) , (3d)
where ε denotes the linearized Green-Lagrange strain tensor, σ the linearized Cauchy stress tensor, C the Hooke material tensor, and 𝑛 (𝑠) the outward unit normal to 𝜕Ω (𝑠) . The external loads are of two kinds: 𝑏 ext denotes body forces, and 𝑡 ext stands for tractions prescribed at the boundary Γ N of the body. Imposed displacements 𝑢 0 over the boundary Γ D are also considered (such that the total problem is well posed). In order to recover the linear elasticity problem posed over the full domain Ω, kinematic compatibility and equilibrium of the tractions along the subdomain interfaces Γ should be added. For every couple (𝑠 1 , 𝑠 2 ) of neighboring subdomains, this leads to write:

𝑢 (𝑠 1 ) -𝑢 (𝑠 2 ) = 0 over Γ (𝑠 1 ,𝑠 2 ) , (4a)
σ (𝑠 1 ) 𝑛 (𝑠 1 ) + σ (𝑠 2 ) 𝑛 (𝑠 2 ) = 0 over Γ (𝑠 1 ,𝑠 2 ) ,

with Γ (𝑠 1 ,𝑠 2 ) = 𝜕Ω (𝑠 1 ) ∩ 𝜕Ω (𝑠 2 ) ⊂ Γ, 𝑠 1 < 𝑠 2 , and ‖Γ (𝑠 1 ,𝑠 2 ) ‖ > 0. This coupled problem is solved approximately using a finite element approach. To this purpose, let us recall the variational forms that come with the problem defined in Equation (3):

a (𝑠) (𝑢, 𝑣) = ∫ Ω (𝑠) σ (𝑠) (𝑢) ∶ ε(𝑣) dx, (5a) b
(𝑠) (𝑣) = ∫ Ω (𝑠) 𝑏 ext ⋅ 𝑣 dx + ∫ Γ N ∩𝜕Ω (𝑠) 𝑡 ext ⋅ 𝑣 dl. (5b)
At this stage, let us make the following important assumption: We consider that the finite element spaces defined over the subdomains (denoted as 𝒲 ℎ(𝑠) and defined such that 𝒲 ℎ(𝑠) ⊂ [ 𝐻 1 (Ω (𝑠) ) ] 𝑑 ) are all built starting from the same functions. More specifically, we introduce a set of basis functions 𝑁 𝑖 ∶ Ω ref → ℝ, 𝑖 = 1, … , 𝑛, where Ω ref denotes the domain associated to the reference unit-cell, and we assume that for each subdomain there exits a smooth mapping 𝐺 (𝑠) ∶ Ω ref → Ω (𝑠) with a smooth inverse that maps each point 𝛏 in Ω ref to a point 𝐱 in Ω (𝑠) (see Figure 2(right)). Thus, we expand the local finite element spaces in terms of basis functions, as:

𝒲 ℎ(𝑠) = span { 𝑁 𝑖 ∘ 𝐺 (𝑠)-1 } 𝑖=1,…,𝑛 . (6)
Finally, let us decompose these local spaces as:

𝒲 ℎ(𝑠) = 𝒲 ℎ(𝑠) Γ ⊕ 𝒲 ℎ(𝑠) 𝐼 , (7)
where 𝒲 ℎ(𝑠) Γ

is the finite element trace space (i.e., made of the basis functions that do not vanish at the subdomain boundary) and 𝒲 ℎ(𝑠) 𝐼 is the finite element space associated to the interior part (i.e., ∀𝑣 ∈ 𝒲 ℎ(𝑠) 𝐼 , 𝑣 |𝜕Ω (𝑠) = 0).

Dual-Primal continuity constraints

The local spaces 𝒲 ℎ(𝑠) can be combined to define a working global space for the total solution:

𝒲 ℎ = 𝑁 ⨁ 𝑠=1 𝒲 ℎ(𝑠) = 𝒲 ℎ 𝐼 ⊕ 𝒲 ℎ Γ , with 𝒲 ℎ 𝐼 = 𝑁 ⨁ 𝑠=1 𝒲 ℎ(𝑠) 𝐼 , 𝒲 ℎ Γ = 𝑁 ⨁ 𝑠=1 𝒲 ℎ(𝑠) Γ . (8)
However, the members of 𝒲 ℎ are in general not continuous across the subdomain interfaces. Thus, an alternative choice would be to look for solutions in the following subset:

Ŵ ℎ = { 𝑢 ℎ ∈ 𝒲 ℎ ∶ 𝑢 ℎ is continuous across Γ } . (9)
The standard FETI method [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] iterates in space 𝒲 ℎ , whereas the alternative Balancing Domain Decomposition (BDD) method [START_REF] Mandel | Balancing domain decomposition[END_REF] iterates in space Ŵ ℎ . Within FETI-DP approaches [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method-part I: A faster alternative to the two-level FETI method[END_REF], the main idea consists in iterating in an intermediate space between 𝒲 ℎ and Ŵ ℎ , where so-called primal continuity constraints hold throughout the iterations. Those primal constraints are defined such that:

1. each subdomain problem becomes invertible (which is not the case when working in 𝒲 ℎ ), 2. good convergence bound can be obtained for the DD solver. We denote this space as W ℎ , for which we have Ŵ ℎ ⊂ W ℎ ⊂ 𝒲 ℎ . It can be defined, to some extent, as:

W ℎ = { 𝑢 ℎ ∈ 𝒲 ℎ ∶ 𝑢 ℎ verifies the primal continuity constraints across Γ } . ( 10
)
The introduced primal continuity constraints only affect a subset of interface variables (defined as primal variables).

The primal continuity constraints depend on the spatial dimension and on the heterogeneity of the problem that is considered (see, e.g., [START_REF] Klawonn | Dual-primal feti methods for linear elasticity[END_REF] for a theoretical study). In this work, we will simply consider the corners of the subdomains for the primal variables [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method-part I: A faster alternative to the two-level FETI method[END_REF][START_REF] Mandel | On the convergence of a dual-primal substructuring method[END_REF], which appears sufficient in our numerical experiments to obtain a scalable method. An illustration of the different solution spaces is given in Figure 3. We can eventually split the interface variables into two categories, the primal and dual variables, and, then, the FETI-DP space W ℎ can be viewed as:

W ℎ = 𝒲 ℎ 𝐼 ⊕ W ℎ Γ , with W ℎ Γ = 𝒲 ℎ 𝐷 ⊕ Ŵ ℎ 𝑃 . (11)
Hence, the interface space is composed of the continuous members of the primal space (i.e., space Ŵ ℎ 𝑃 ) and the standard dual space (i.e., space 𝒲 ℎ 𝐷 ) with possibly disconnected (or discontinuous) members.

Additionally, Lagrange multipliers are introduced in order to recover the continuity over the dual space along the iterative algorithm. Thus, we enforce continuity for all dual interface variables in a weak sense. As we placed ourselves in a fully matching setting (since the same reference mesh is used for each cell, see Figures 2 and3), we can directly measure the interface jump by considering the differences between couples of displacement Degrees Of Freedom (DOF), and assign a Lagrange multiplier to each of them to introduce a continuity constraint. Importantly, we also enforce the Dirichlet boundary conditions over the dual variables through Lagrange multipliers in order to keep local dual spaces with the same cardinality.

Matrix formulation and Saddle point problem

Let us now highlight how the imposition of the Primal-Dual coupling constraints is done in practice by employing a matrix notation. First, we partition every set of local displacement DOF into three categories: the internal, the dual, and the primal DOF denoted with subscripts ⋅ 𝑖 , ⋅ 𝑑 , and ⋅ 𝑝 , respectively. This reads:

𝐮 (𝑠) = ( 𝐮 (𝑠) 𝑖 , 𝐮 (𝑠) 𝑑 , 𝐮 (𝑠) 𝑝 ) , with 𝐮 (𝑠) 𝑖 ∈ ℝ 𝑛 𝑖 , 𝐮 (𝑠) 𝑑 ∈ ℝ 𝑛 𝑑 , 𝐮 (𝑠) 𝑝 ∈ ℝ 𝑛 𝑝 . ( 12
)
We additionally denote by 𝐮 𝑑 ) ∈ ℝ 𝑛 𝑟 the non-primal DOF (also called the remaining DOF), and we introduce a restriction operator 𝐓 𝑑𝑟 ∈ ℝ 𝑛 𝑑 ×𝑛 𝑟 that enables to extract the dual DOF from the remaining ones as 𝐮 (𝑠) 𝑑 = 𝐓 𝑑𝑟 𝐮 (𝑠) 𝑟 . All the local primal unknowns are assembled globally into 𝐮 𝑃 , and the remaining ones are simply concatenated which leads to define the total vector of unknowns 𝐮 ∈ ℝ 𝑛 𝑈 as:

𝐮 = ( 𝐮 (1) 𝑟 , 𝐮 (2) 𝑟 , … , 𝐮 (𝑁) 𝑟 , 𝐮 𝑃 ) . (13)
In the following we will also denote by 𝑛 𝑅 and 𝑛 𝑃 the total number of remaining and primal DOF (𝑛 𝑈 = 𝑛 𝑅 + 𝑛 𝑃 ), respectively.

In order to enforce the continuity between the dual interfaces, we introduce a (discrete) jump operator of the following form:

𝐵 ∶ W ℎ → ℝ 𝐿 ; 𝑢 ℎ ↦ 𝐵(𝑢 ℎ ) = 𝐁𝐮 -𝐝, ( 14
)
where 𝐿 is the total number of discrete continuity equations and Dirichlet constraints. Vector 𝐝 ∈ ℝ 𝐿 has zero entries where the coupling conditions concern a subdomain-to-subdomain constraint, and integrates the Dirichlet conditions (Equation (3c)) otherwise. The coupling matrix 𝐁 ∈ ℝ 𝐿×𝑛 𝑈 is formed by local (dual) Boolean operators 𝐁 (𝑠) 𝑑 ∈ ℝ 𝐿×𝑛 𝑑 which are rectangular sparse matrices with only 0, 1, -1 entries. More precisely, it reads:

𝐁 = [ 𝐁 𝑅 𝟎 𝐿𝑃 ] with 𝐁 𝑅 = [ 𝐁 (1) 𝑟 … 𝐁 (𝑁) 𝑟 ] , and 𝐁 (𝑠) 𝑟 = 𝐁 (𝑠) 𝑑 𝐓 𝑑𝑟 . (15)
The interested readers may refer to [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF][START_REF] Jolivet | Domain decomposition methods. Application to high-performance computing[END_REF] to find details on the construction of these operators. Finally, we define a set of pointwise Lagrange multipliers 𝛌 ∈ ℝ 𝐿 to enforce the continuity constraints 𝐁𝐮 -𝐝 = 𝟎. It leads to formulate the following saddle problem:

[ 𝐊 𝐁 ⊤ 𝐁 𝟎 ] ( 𝐮 𝛌) = ( 𝐟 𝐝) , (16)
where, due to the dual-primal domain decomposition, we have:

𝐊 = [ 𝐊 𝑅𝑅 𝐊 𝑅𝑃 𝐊 𝑃 𝑅 𝐊 𝑃 𝑃 ] , 𝐟 = ( 𝐟 𝑅 𝐟 𝑃 ) . ( 17
)
As already stated, the primal DOF are assembled globally within FETI-DP approaches. This is a common task in finite elements and can be done by employing a global-local assembly operator 𝐀 (𝑠) 𝑝 ∈ ℝ 𝑛 𝑃 ×𝑛 𝑝 per subdomain. Therefore, the quantities in Equation ( 17) are obtained by:

𝐊 𝑃 𝑃 = 𝑁 ∑ 𝑠=1 𝐀 (𝑠) 𝑝 𝐊 (𝑠) 𝑝𝑝 𝐀 (𝑠)⊤ 𝑝 , 𝐟 𝑃 = 𝑁 ∑ 𝑠=1 𝐀 (𝑠) 𝑝 𝐟 (𝑠) 𝑝 , (18)
and:

𝐊 𝑅𝑅 = ⎡ ⎢ ⎢ ⎣ 𝐊 (1) 𝑟𝑟 ⋱ 𝐊 (𝑁) 𝑟𝑟 ⎤ ⎥ ⎥ ⎦ , 𝐊 𝑅𝑃 = ⎡ ⎢ ⎢ ⎣ 𝐊 (1) 𝑟𝑃 ⋮ 𝐊 (𝑁) 𝑟𝑃 ⎤ ⎥ ⎥ ⎦ , with 𝐊 (𝑠) 𝑟𝑃 = 𝐊 (𝑠) 𝑟𝑝 𝐀 (𝑠)⊤ 𝑝 , 𝐊 𝑃 𝑅 = 𝐊 ⊤ 𝑅𝑃 . (19)
Hence, regarding the structure of the total matrix 𝐊, it is interesting to notice the presence of the block diagonal term 𝐊 𝑅𝑅 made of the local problems 𝐊 (𝑠)

𝑟𝑟 . In addition to being symmetric positive definite, this term is preponderant as 𝑛 𝑅 ≫ 𝑛 𝑃 (since the number of primal constraints is usually low). Thus, the main idea behind FETI-DP approaches (and DD methods in general), is to take advantage of this particular structure.

FETI-DP algorithm

The common starting point of FETI solvers consists in removing the displacement unknowns from system [START_REF] Elber | Precise Construction of Micro-structures and Porous Geometry via Functional Composition[END_REF] in order to formulate an interface problem where only the Lagrange multipliers remain. With the FETI-DP approach, this is done by firstly condensing the remaining DOF by using the relation:

𝐮 𝑅 = 𝐊 -1 𝑅𝑅 ( 𝐟 𝑅 -𝐊 𝑅𝑃 𝐮 𝑃 -𝐁 ⊤ 𝑅 𝛌 ) , (20)
and then the primal DOF by using the relation:

𝐮 𝑃 = 𝐒 -1 𝑃 𝑃 (( 𝐟 𝑃 -𝐊 𝑃 𝑅 𝐊 -1 𝑅𝑅 𝐟 𝑅) + [ 𝐊 𝑃 𝑅 𝐊 -1 𝑅𝑅 𝐁 ⊤ 𝑅 ] 𝛌 ) , (21)
where 𝐒 𝑃 𝑃 ∈ ℝ 𝑛 𝑃 ×𝑛 𝑃 denotes the primal Schur complement. It consists in the assembly of all the local primal Schur complements 𝐒 (𝑠) 𝑝𝑝 :

𝐒 𝑃 𝑃 = 𝑁 ∑ 𝑠=1 𝐀 (𝑠) 𝑝 𝐒 (𝑠) 𝑝𝑝 𝐀 (𝑠)⊤ 𝑝 , with 𝐒 (𝑠) 𝑝𝑝 = 𝐊 (𝑠) 𝑝𝑝 -𝐊 (𝑠) 𝑝𝑟 𝐊 (𝑠)-1 𝑟𝑟 𝐊 (𝑠) 𝑟𝑝 , (22)
and can be written in a more compact way as: 𝐒 𝑃 𝑃 = 𝐊 𝑃 𝑃 -𝐊 𝑃 𝑅 𝐊 -1 𝑅𝑅 𝐊 𝑅𝑃 . Using equations ( 20) and ( 21) enables to transform the continuity constraint between the dual interfaces (i.e., the second equation in the saddle point problem in Equation ( 16)) into the following interface problem:

𝐁𝐮 = 𝐝 ⇔ 𝐅𝛌 = d , (23)
with:

𝐅 = 𝐁 𝑅 𝐊 -1 𝑅𝑅 𝐁 ⊤ 𝑅 + 𝐁 𝑅 [ 𝐊 -1 𝑅𝑅 𝐊 𝑅𝑃 𝐒 -1 𝑃 𝑃 𝐊 𝑃 𝑅 𝐊 -1 𝑅𝑅 ] 𝐁 ⊤ 𝑅 , (24) d
= -𝐝 + 𝐁 𝑅 𝐊 -1 𝑅𝑅 𝐟 𝑅 -𝐁 𝑅 𝐊 -1 𝑅𝑅 𝐊 𝑅𝑃 𝐒 -1 𝑃 𝑃 ( 𝐟 𝑃 -𝐊 𝑃 𝑅 𝐊 -1 𝑅𝑅 𝐟 𝑅) . (25)
𝐅, which is called the global dual Schur complement, may also be expressed in a more compact form as:

𝐅 = 𝐁𝐊 -1 𝐁 ⊤ . (26)
The global dual Schur complement is made of two terms as shown in Equation [START_REF] Have | Algebraic domain decomposition methods for highly heterogeneous problems[END_REF]. The first term collects the local Neumann schur complements 𝐅 (𝑠) 𝑑𝑑 :

𝐁 𝑅 𝐊 -1 𝑅𝑅 𝐁 ⊤ 𝑅 = 𝑁 ∑ 𝑠=1 𝐁 (𝑠) 𝑑 𝐅 (𝑠) 𝑑𝑑 𝐁 (𝑠)⊤ 𝑑 , with 𝐅 (𝑠) 𝑑𝑑 = 𝐓 𝑑𝑟 𝐊 (𝑠)-1 𝑟𝑟 𝐓 ⊤ 𝑑𝑟 , (27)
and the second term is the coarse problem correction associated to the primal variables. Let us further detail key quantities in Equations ( 24) and ( 25) that involve the local dual and primal solutions denoted 𝐔 (𝑠)

𝑟𝑑 and 𝐔 (𝑠) 𝑟𝑝 :

𝐁 𝑅 𝐊 -1 𝑅𝑅 = [ 𝐁 (1) 𝑑 𝐔 (1)⊤ 𝑟𝑑 … 𝐁 (𝑁) 𝑑 𝐔 (𝑁)⊤ 𝑟𝑑 ] , with 𝐔 (𝑠) 𝑟𝑑 = 𝐊 (𝑠)-1 𝑟𝑟 𝐓 ⊤ 𝑑𝑟 , (28)
𝐊 -1 𝑅𝑅 𝐊 𝑅𝑃 = ⎡ ⎢ ⎢ ⎣ 𝐔 (1) 𝑟𝑝 𝐀 (1)⊤ 𝑝 ⋮ 𝐔 (𝑁) 𝑟𝑝 𝐀 (𝑁)⊤ 𝑝 ⎤ ⎥ ⎥ ⎦ , with 𝐔 (𝑠) 𝑟𝑝 = 𝐊 (𝑠)-1 𝑟𝑟 𝐊 (𝑠) 𝑟𝑝 . (29)
Secondly, the main idea is to solve the interface problem [START_REF] Harris | Array programming with NumPy[END_REF] with an iterative solver, as for instance the preconditioned conjugate gradient since 𝐅 is symmetric positive definite. This way, only matrix-vector products need to be performed. This can be done without building the dual Schur complement 𝐅. Indeed, in practice only the local primal solutions, the local primal Schur complements, and the global Schur complement are explicitly built. Furthermore, the application of 𝐅 to a vector involves the solution of local independent problems, which can be done in parallel.

The last thing to do thus consists in defining a preconditioner for system [START_REF] Harris | Array programming with NumPy[END_REF]. Usually, the so-called Dirichlet preconditioner is considered in linear elasticity [START_REF] Rixen | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF]. It is built as follows:

𝐌 -1 D = 𝑁 ∑ 𝑠=1 𝐃 (𝑠) 𝐁 (𝑠) 𝑑 𝐒 (𝑠) 𝑑𝑑 𝐁 (𝑠)⊤
𝑑 𝐃 (𝑠) , with 𝐒

(𝑠) 𝑑𝑑 = 𝐊 (𝑠) 𝑑𝑑 -𝐊 (𝑠) 𝑑𝑖 𝐊 (𝑠)-1 𝑖𝑖 𝐊 (𝑠) 𝑖𝑑 , (30)
where 𝐃 (𝑠) ∈ ℝ 𝐿×𝐿 are diagonal weight matrices that verify

∑ 𝑁 𝑠=1 𝐁 (𝑠) 𝑑 𝐁 (𝑠)⊤
𝑑 𝐃 (𝑠) = 𝐈 𝐿 and can be defined in order to better balance the interface correction in case of material heterogeneities and redundancies in the compatibility constraints, as explained in, e.g., [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF]. More elaborated strategies exist to handle very high heterogeneities, as for instance, the use of Generalized Eigenproblems in the Overlaps [START_REF] Spillane | Automatic spectral coarse spaces for robust finite element tearing and interconnecting and balanced domain decomposition algorithms[END_REF] or multipreconditioned Krylov solvers [START_REF] Bovet | Two-level adaptation for adaptive multipreconditioned FETI[END_REF][START_REF] Spillane | An adaptive MultiPreconditioned conjugate gradient algorithm[END_REF]. Once the interface problem ( 23) is solved, one can recover the displacement solution 𝐮 by using successively Equations ( 21) and [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF].

Overall, returning to the non-condensed form ( 16) and making use of ( 26), it is possible to rewrite this domain decomposition solution in a matrix from as follows:

( 𝐮 𝛌) = [ 𝐈 -𝐊 -1 𝐁 ⊤ 𝐈 ] ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ Recover full sol. [ 𝐊 -1 -𝐅 -𝚙𝚌𝚐 ] ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ Interface pb. [ 𝐈 -𝐁𝐊 -1 𝐈] ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ Condense ( 𝐟 𝐝) , (31)
where 𝐅 -𝚙𝚌𝚐 stands for the application of 𝐅 -1 via a preconditioned conjugate gradient (i.e., solving the interface problem).

Block preconditioners saddle point problems

Applying the standard algorithm explained in Section 2.2.4 requires to solve exactly (i.e., close to machine precision) the local problems (within 𝐊 -1 𝑅𝑅 ) and the coarse problem (𝐒 -1 𝑃 𝑃 ) at every iteration, see Equations ( 24) and [START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF]. In practice, it is commonly done by factorizing all these operators a priori, and then performing backward-forward substitutions during the iterative algorithm. It is worth nothing that 𝐊 -1 𝑅𝑅 involves the solution of local independent problems; therefore, it can be done in parallel. However, this still raises issues, for instance, in the case of large coarse problem (or large local problem) which leads to high memory usage (and bad parallel efficiency). In the specific case of lattice structures, one may seek to exploit the similarities between the local problems. If exact local solutions are mandatory, then it seems not straightforward to take advantage of these similarities (unless all cells are identical).

Therefore, we will rather consider an alternative approach that consists in working directly with the saddle point problem [START_REF] Elber | Precise Construction of Micro-structures and Porous Geometry via Functional Composition[END_REF], and to use a suitable iterative solver (GMRES for example) together with a block preconditioner. It is possible to build block preconditioners dedicated to the saddle point problem ( 16) that uses the same ingredients as those involved in standard FETI approaches (local and coarse solutions). However, since the treatment now appears in the preconditioner, it is possible to approximate the local problems (and also the coarse problem), solving them inexactly. This leads to the so-called inexact FETI-DP approaches [START_REF] Bosy | A domain decomposition method for Isogeometric multi-patch problems with inexact local solvers[END_REF][START_REF] Klawonn | Inexact FETI-DP methods[END_REF][START_REF] Rheinbach | Parallel iterative substructuring in structural mechanics[END_REF].

The starting point to build block preconditioners for saddle point problems of type ( 16) consists in performing the following block LDL ⊤ factorization:

[ 𝐊 𝐁 ⊤ 𝐁 𝟎 ] = [ 𝐈 𝐁𝐊 -1 𝐈] [ 𝐊 -𝐅] [ 𝐈 𝐊 -1 𝐁 ⊤ 𝐈 ] , (32)
which is possible as 𝐊 is a symmetric positive definite matrix and 𝐁 has full rank. Then, standard block preconditioners consist in using different sub-preconditioners, i.e., different approximations of each blocks, see for instance Notay [START_REF] Notay | A new analysis of block preconditioners for saddle point problems[END_REF].

Indeed, block preconditioners are of the form:

Å-1 = [ 𝐈 -Ů ⊤ 𝐈 ] [ K-1 -F -1 ] [ 𝐈 -V 𝐈] , (33)
where K-1 is a preconditioner for 𝐊, F -1 is a preconditioner for 𝐅, and Ů and V are approximations of 𝐁𝐊 -1 . Several choices for these sub-preconditioners can be considered. In the context of FETI-based approaches, we could select:

• Standard FETI-DP (see Equation ( 31)

): K-1 = 𝐊 -1 , F -1 = 𝐅 -𝚙𝚌𝚐 , Ů = V = 𝐁𝐊 -1 ,
• Block triangular Inexact FETI-DP from Klawonn and Rheinbach [START_REF] Klawonn | Inexact FETI-DP methods[END_REF]:

K-1 = 𝐊 -𝙰𝙼𝙶 , F -1 = 𝐌 -𝙰𝙼𝙶 D , Ů = 𝟎, V = 𝐁𝐊 -1
, where (⋅) -𝙰𝙼𝙶 means that the solution is obtained with an Algebraic MultiGrid method,

• Block diagonal Inexact IETI-DP based on Bosy et al. [START_REF] Bosy | A domain decomposition method for Isogeometric multi-patch problems with inexact local solvers[END_REF]:

K-1 = 𝐊 -𝙵𝙳 , F -1 = 𝐌 -𝙵𝙳 D , Ů = 𝟎, V = 𝟎
, where (⋅) -𝙵𝙳 means that the inverses are approximated via the Fast Diagonalization method.

In this work, we seek to introduce sub-preconditioners in Equation ( 33) that exploit the similarities between all local subproblems, i.e., the fact that:

𝐊 (1) ≃ 𝐊 (2) ≃ … ≃ 𝐊 (𝑁) .

(

) 34
The idea is therefore to use only a small set of local stiffness matrices to efficiently approximate all the local solutions, which is the object of next section. Let us now introduce the matrix manifold comprising all stiffness matrices associated to every sub-problem of our total DD problem: (1) , … , 𝐊 (𝑁) } .

Benefiting from the repetitiveness of lattices

ℳ K = { 𝐊
Due to the similarities between these matrices, we suspect the manifold ℳ K to be of low dimension; i.e., it can be well represented by the span of a low number 𝑁 𝔯𝔟 ≪ 𝑁 of reference matrices, denoted as 𝐊 [𝑘] , with 𝑘 = 1, … , 𝑁 𝔯𝔟 . Thus, we assume that there is a Reduced Basis (RB) space: 𝕂 = span { 𝐊 [1] , … , 𝐊 [𝑁 𝔯𝔟 ] } , [START_REF] Klawonn | Inexact FETI-DP methods[END_REF] such that:

∀𝐊 (𝑠) ∈ ℳ K , ∃ K(𝑠) ∈ 𝕂 𝑠.𝑡. ‖𝐊 (𝑠) -K(𝑠) ‖ 𝑣 ≤ tol, (37)
where ‖ ⋅ ‖ 𝑣 is an appropriate norm. In Section 4.2, we will explain in detail how to build efficiently the reduced space 𝕂. At this stage, let us just mention that we will make use of a greedy approach which generates an interpolatory basis; that is, the reference matrices 𝐊 [𝑘] correspond to local matrices 𝐊 (𝑠) of appropriately selected subdomains denoted 𝜎 𝑘 , 𝑘 = 1, … , 𝑁 𝔯𝔟 . In other words: where the 𝛼 (𝑠) 𝑘 ∈ ℝ are coefficients which depend on the subdomain 𝑠. Figure 4 illustrates the result of such a process, and, once again, Section 4.2 will detail the construction of these approximations, including the selection of the principal subdomains.

K(𝑠) = 𝑁 𝔯𝔟 ∑ 𝑘=1 𝛼 (𝑠) 𝑘 𝐊 (𝜎 𝑘 ) , (38)
Similarities among subdomains are also reflected in all the different local quantities involved in the FETI-DP method (see again Section 2. Based on the reduced basis of the local stiffness matrices and the selection of principal subdomains, we also consider reduced bases of all these local quantities by simply replacing 𝐊 (𝑠) by 𝐊 (𝜎 𝑘 ) in the above equations, which leads to:

𝕌 𝑟𝑝 = span { 𝐔 (𝜎 1 )
𝑟𝑝 , … , 𝐔 Then, the idea is to use these RB in the inexact FETI-DP preconditioner by approximating all the local operators.

Remark 2. All the reduced bases for the DD operators have the same size which is inherited from the reduced basis of the stiffness matrices [START_REF] Klawonn | Inexact FETI-DP methods[END_REF]. Indeed, we do not place ourselves in a standard ROM configuration where the reduced bases are built in an offline phase. We cannot pay the price of generating snapshots of all local DD operators as it would take more time than simply solving the FE system. Instead, we build these operators for the principal cells involved in 𝕂 and use them to generate the other reduced bases (see again Equation [START_REF] Kleiss | IETI-isogeometric tearing and interconnecting[END_REF]). Therefore, we will call upon a (global) iterative process to recover the exact solution (instead of in one go as for the standard FETI-DP approach, see Equation (31)).

Approximating the local operators

Let us invoke the main ideas behind ROM [START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF][START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF] in order to benefit from the RB spaces. We project each subdomain local operator into their corresponding reduced bases, or, in other words, we aim at finding the coefficients of the local operators in the reduced bases.

To this purpose, we solve two minimization problems per subdomain, which enable to approximate the primal and dual solutions denoted previously 𝐔 (𝑠) 𝑟𝑝 and 𝐔 (𝑠) 𝑟𝑑 , respectively:

Ů (𝑠) 𝑟𝑝 = argmin 𝐔 𝑟𝑝 ∈𝕌 𝑟𝑝 ‖ 𝐔 (𝑠) 𝑟𝑝 -𝐔 𝑟𝑝 ‖ 2 𝐸 (𝑠) and Ů (𝑠) 𝑟𝑑 = argmin 𝐔 𝑟𝑑 ∈𝕌 𝑟𝑑 ‖ 𝐔 (𝑠) 𝑟𝑑 -𝐔 𝑟𝑑 ‖ 2 𝐸 (𝑠) , (40)
where ‖ ⋅ ‖ 2 𝐸 (𝑠) denotes an energy norm which for a vector 𝐯 𝑟 ∈ ℝ 𝑟 and a column matrix 𝐕 𝑟𝑙 ∈ ℝ 𝑟×𝑙 and read as, respectively: In case of column matrices, the energy norm can also be seen as the sum of the energy norms of each column vectors. With such a choice, solving Equation ( 40) is equivalent to solve small 𝑁 𝔯𝔟 -by-𝑁 𝔯𝔟 linear systems: 𝐀 (𝑠) 𝛑 (𝑠) = 𝐛 (𝑠) and 𝐌 (𝑠) 𝛅 (𝑠) = 𝐝,

‖𝐯 𝑟 ‖ 2 𝐸 (𝑠) = 1 2 𝐯 ⊤ 𝑟 𝐊 (
where 𝐀 (𝑠) , 𝐌 (𝑠) ∈ ℝ 𝑁 𝔯𝔟 ×𝑁 𝔯𝔟 and 𝐛 (𝑠) , 𝐝 ∈ ℝ 𝑁 𝔯𝔟 are given by:

𝐀 (𝑠) [𝑖, 𝑗] = tr ( 𝐔 (𝜎 𝑖 )⊤ 𝑟𝑝 𝐊 (𝑠) 𝑟𝑟 𝐔 (𝜎 𝑗 )
𝑟𝑝 ) , 𝐛 (𝑠) [𝑖] = tr ( 𝐔

(𝜎 𝑖 )⊤ 𝑟𝑝 𝐊 (𝑠) 𝑟𝑝 ) , (42a)
𝐌 (𝑠) [𝑖, 𝑗] = tr ( 𝐔

(𝜎 𝑖 )⊤ 𝑟𝑑 𝐊 (𝑠) 𝑟𝑟 𝐔 (𝜎 𝑗 ) 𝑟𝑑 ) , 𝐝[𝑖] = tr ( 𝐔 (𝜎 𝑖 )⊤ 𝑟𝑑 𝐓 𝑟𝑑 ) . (42b)
We used the syntax [𝑖, 𝑗] and [𝑖] for the index notation in order to distinguish it from the size/DOF related notation needed for the description of the FETI algorithm.

Once the coefficients 𝛂 (𝑠) (defined in Equation ( 38)), 𝛑 (𝑠) , and 𝛅 (𝑠) are known, we can obtain the ROM of all the local DD operators as:

Ů (𝑠) 𝑟𝑝 = 𝕌 𝑟𝑝 𝛑 (𝑠) , S (𝑠) 𝑝𝑝 = sym ( K(𝑠) 𝑝𝑝 - K(𝑠)⊤ 𝑟𝑝 Ů (𝑠) 𝑟𝑝 ) , Ů (
𝑠) 𝑟𝑑 = 𝕌 𝑟𝑑 𝛅 (𝑠) , F (𝑠) 𝑑𝑑 = 𝔽 𝑑𝑑 𝛅 (𝑠) , S (𝑠) 𝑑𝑑 = 𝕊 𝑑𝑑 𝛂 (𝑠) . [START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF] In the above Equation, we made an abuse of notation for the sake of clarity: For instance, 𝕌 𝑟𝑝 𝛑 (𝑠) should be understood as the linear combination of the basis elements spanning the space 𝕌 𝑟𝑝 with the scalars collected in vector 𝛑 (𝑠) . And the same applies to 𝕌 𝑟𝑑 𝛅 (𝑠) , 𝔽 𝑑𝑑 𝛅 (𝑠) , and 𝕊 𝑑𝑑 𝛂 (𝑠) .

Algorithm 1

Steps when applying the iFETI-DP preconditioner to an input vector.

Input vector:

In 𝑑𝑑 are approximated directly through the ROM coefficients 𝛂 (𝑠) of the stiffness matrices. We could improve this approximation by employing an expression similar to the one of the local primal Schur complements. However, this would add significant extra computational costs (projection of 𝐊 (𝑠)-1 𝑖𝑖 𝐊 (𝑠) 𝑖𝑑 using a similar approach as in [START_REF] Korshunova | Image-based material characterization of complex microarchitectured additively manufactured structures[END_REF]) for little improvements as these operators are involved in the preconditioning of the interface problem only (see the expression of the Dirichlet preconditioner in Equation (30)). Our numerical experiments confirm that it is a sufficient choice.

Remark 4. Contrary to Remark 3, we could have chosen to approximate the local primal Schur complements as

S (𝑠) 𝑝𝑝 = 𝕊 𝑝𝑝 𝛂 (𝑠) as done for the local Dirichlet Schur complements. However, since we already paid the price to approximate the primal solutions properly by determining the coefficients 𝛑 (𝑠) , it appears natural to use these best approximations.

ROM-based inexact FETI-DP

We finally combine the FETI-DP approach together with the previously introduced ROM ingredients into the following block preconditioner:

Å-1 𝚒𝙵𝙴𝚃𝙸𝚍𝚙 = [ 𝐈 -Ů ⊤ 𝐈 ] [ K-𝚁𝙾𝙼 -F -𝚙𝚌𝚐 ] [ 𝐈 -Ů 𝐈] , (44)
where (⋅) -𝚁𝙾𝙼 indicates that the solution is obtained in a ROM-based, and (⋅) -𝚙𝚌𝚐 indicates that the solution is obtained via a preconditioned conjugate gradient (as it was the case for the standard FETI-DP approach, see again Equation ( 31)).

Algorithm 1 shows the main steps involved during the application of this inexact FETI-DP preconditioner to an input vectors, and the following sections detail the different sub-blocks.

Coarse problem

The global primal Schur complement 𝐒 𝑃 𝑃 (which appears in both K-𝚁𝙾𝙼 and F -𝚙𝚌𝚐 ) is approximated using the ROM of the local primal Schur complements S (𝑠) 𝑝𝑝 as defined in Equation ( 43). In the offline phase, once assembled, we factorize it so that only forward-backward substitutions are performed during the iterations of the solver (i.e., whenever one needs to solve the coarse problem). The factorization is done "exactly" via a Cholesky decomposition as for the standard FETI-DP approach, but we could eventually make use of approximative solvers for the coarse problem as we are in the context of inexact FETI methods. For instance, one would go for approximative solvers to maintain parallel scalabilty for large number of subdomains, as motivated in the initial works on inexact FETI-DP, see Klawonn and Rheinbach [START_REF] Klawonn | Inexact FETI-DP methods[END_REF]. This is for sure interesting for lattice structures as they have numerous subdomains, but this has not been investigated in the scope of this work.

Algorithm 2

Steps when applying the sub-preconditioner K-𝚁𝙾𝙼 to an input vector.

Input vector:

In , with 𝐱 𝑅 ∈ ℝ 𝑛 𝑅 , 𝐱 𝑃 ∈ ℝ 𝑛 𝑃 .

Inexact local solves

The sub-preconditioner K-𝚁𝙾𝙼 in Equation ( 44) reads as:

K-𝚁𝙾𝙼 = [ 𝐈 𝑅 -Ů 𝑅𝑃 𝐈 𝑃 ] [ K-𝚁𝙾𝙼 𝑅𝑅 S -1 𝑃 𝑃 ] [ 𝐈 𝑅 -Ů 𝑃 𝑅 𝐈 𝑃 ] . (45)
This block form originates from a block LDL ⊤ factorization of Equation ( 17), as done for the complete saddle point problem (see Equation ( 32)). One could rewrite Equations ( 20) and [START_REF] Greif | GMRES with multiple preconditioners[END_REF] in this block format. The novelty does not lie here but in the sub-blocks.

The coarse problem is handled as described previously (see Section 3.2.1). On the other hand, the column matrix Ů 𝑃 𝑅 assembles the ROM of the local primal solutions Ů (𝑠) 𝑟𝑝 , and thus is built quickly (once the ROM is performed). Finally, operator K-𝚁𝙾𝙼 𝑅𝑅 stands for the ROM-based solution of the local problems (recall the definition of 𝐊 𝑅𝑅 in Equation ( 19)). Indeed, a matrix-vector product of type 𝐱 𝑅 = K-𝚁𝙾𝙼 𝑅𝑅 𝐯 𝑅 involves the solution of 𝑁 independent local problems which, for each subdomain 𝑠, is done similarly as the approximation of the local DD operators described in Section 3.1:

𝐑 (𝑠) = [ 𝐊 (𝜎 1 )-1 𝑟𝑟 𝐯 (𝑠) 𝑟 , … , 𝐊 (𝜎 𝑁 𝔯𝔟 )-1 𝑟𝑟 𝐯 (𝑠) 𝑟 ] ∈ ℝ 𝑟×𝑁 𝔯𝔟 , (46a)
𝐱 (𝑠) 𝑟 = argmin 𝐱 𝑟 ∈span(𝐑 (𝑠) ) ‖ 𝐊 (𝑠)-1 𝑟𝑟 𝐯 (𝑠) 𝑟 -𝐱 𝑟 ‖ 𝐸 (𝑠) = 𝐑 (𝑠) [ 𝐑 (𝑠)⊤ 𝐊 (𝑠) 𝑟𝑟 𝐑 (𝑠) ] -1 𝐑 (𝑠)⊤ 𝐯 (𝑠) 𝑟 . (46b)
The first step given by Equation (46a) consists in generating a reduced basis for the local solution. This is done by using the 𝑁 𝔯𝔟 principal local problems. Those principal stiffness matrices are previously factorized during the offline phase, and thus, the generation of the reduced bases 𝐑 (𝑠) for the local solutions is done via forward-backward substitutions only. The second step given by Equation (46b) consists in approximating the solution within the reduced space, as already discussed in Section 3.1.2. Algorithm 2 summarizes the main steps involved in the product between the sub-preconditioner K-𝚁𝙾𝙼 and an input vector, and how it can be implemented in practice. [START_REF] Bridson | A multipreconditioned conjugate gradient algorithm[END_REF][START_REF] Greif | GMRES with multiple preconditioners[END_REF][START_REF] Spillane | An adaptive MultiPreconditioned conjugate gradient algorithm[END_REF] where the idea of generating multiple search directions per iteration is performed similarly to the process depicted in Equation [START_REF] Mandel | On the convergence of a dual-primal substructuring method[END_REF]. Indeed, one can consider using a multipreconditioned conjugate gradient solver [START_REF] Bridson | A multipreconditioned conjugate gradient algorithm[END_REF][START_REF] Spillane | An adaptive MultiPreconditioned conjugate gradient algorithm[END_REF] to solve the local problems up to a predefined tolerance with the 𝐊 (𝜎 𝑘 )-1 as the multiple preconditioners. This could lead to a K-𝙼𝙿𝙲𝙶 version of the subpreconditioner.

Remark 5. One application of K-𝚁𝙾𝙼 is quite expensive but presents several advantages. The main one is that only the principal local problems are involved: There is no need to factorize the 𝑁 local problems as needed with standard FETI-DP (or if one need to solve the local problems exactly at every iteration) but only the 𝑁 𝔯𝔟 principal local problems. This not only means less computational time in the preprocessing step due to the fewer factorizations, but also leads to less computer memory usage which is one of the main bottlenecks when dealing with the high-fidelity, fine-scale simulation of lattice structures.

Remark 6. Let us mention the existence of multipreconditioned iterative solvers

Fast interface solves

Then, we need to define the approximation of the global dual/Neumann Schur complement used in the sub-preconditioner (i.e., F -1 in Equation ( 33)). Following the inexact FETI-DP solver described in [START_REF] Klawonn | Inexact FETI-DP methods[END_REF], one may think of using a Dirichlet preconditioner:

F -1 = M-𝚁𝙾𝙼 D , (47)
where M-𝚁𝙾𝙼 D is the Dirichlet preconditioner (recall Equation ( 30)) but where the solutions of the local problems 𝐊 (𝑠)-1 𝑖𝑖 are approximated using a ROM-based solution similar to K-𝚁𝙾𝙼 in Equation [START_REF] Mandel | On the convergence of a dual-primal substructuring method[END_REF]. Although it is an interesting choice (as it reduces the number of operators to be factorized), we suggest another option in view of better approximating 𝐅 -1 and thus further reducing the number of iterations of the DD solver. This is where we use the ROM of the local Schur complements F (𝑠) 𝑑𝑑 and S (𝑠) 𝑑𝑑 built in Section 3.1.2. By recalling the expressions of the global dual Schur complement 𝐅 and the Dirichlet preconditioner M-1 D defined in Equations ( 24) and ( 30), respectively, we introduce their ROM-based approximations as:

F = 𝑁 ∑ 𝑠=1 𝐁 (𝑠) 𝑑 F (𝑠) 𝑑𝑑 𝐁 (𝑠)⊤ 𝑑 + 𝐁 𝑅 Ů 𝑅𝑃 S -1 𝑃 𝑃 Ů ⊤ 𝑅𝑃 𝐁 ⊤ 𝑅 , (48)
M-1 D = 𝑁 ∑ 𝑠=1 𝐃 (𝑠) 𝐁 (𝑠) 𝑑 S (𝑠) 𝑑𝑑 𝐁 (𝑠)⊤ 𝑑 𝐃 (𝑠) . (49)
Beyond being an approximation, the important difference between 𝐅 and F , and 𝐌 -1 D and M-1 D , is the computational cost required to perform the matrix-vector product with an input vector. Indeed, as the local Schur complements of the principal subdomains are explicitly built, the application of F and M-1 D to a vector only requires matrix-vector products instead of solving local problems, as requested with 𝐅 and 𝐌 -1 D . Based on this observation, we intend to use a preconditioned conjugate gradient to approximate and to solve the interfaces problem, i.e., we approximate 𝐅 -1 ≈ F -𝚙𝚌𝚐 , where F corresponds to [START_REF] Notay | A new analysis of block preconditioners for saddle point problems[END_REF] preconditioned with M-1 D , defined in Equation ( 49). [START_REF] Notay | A new analysis of block preconditioners for saddle point problems[END_REF] and (49) could be rewritten in terms of quantities defined over the principal subdomains. Instead of a sum of 𝑁 terms, it would be a sum of 𝑁 𝔯𝔟 terms. The key point to remember is that those operators are cost-effective, in the sense that matrix-vector products are fast to compute.

Remark 7. The formulations of the ROM-based approximations of the global dual Schur complement and the Dirichlet preconditioner given in Equations

Fast transitions between Lagrange multipliers and displacement DOF

Finally, we approximate the rectangular matrix as Ů ⊤ = K-𝚁𝙾𝙼 𝐁 ⊤ which is applied twice per call of the inexact FETI-DP preconditioner (see Equation ( 44)). It is involved when building the right-hand sides of the interface problems and when recovering the full solution as shown in Algorithm 1. The product between Ů (or Ů ⊤ ) and an input vector involves similar steps than the action of the sub-preconditioner K-𝚁𝙾𝙼 . The difference is that there is no need to generate reduced bases at every call. Indeed, we can use directly the ROM of the dual solutions Ů (𝑠) 𝑟𝑑 which leads to numerically cheap matrix-vector products. Algorithm 3 summarizes the main steps involved during those products between the sub-block Ů and input vectors, and how it can be implemented in practice (it would be similar for Ů ⊤ ).

Algorithm 3

Steps when performing matrix-vector products between Ů and an input vector.

Input vector:

In

Application to multiscale geometric models

Geometric modeling via composition

We consider in this work geometric models of lattice structures based on a multiscale geometric modeling approach introduced in [START_REF] Elber | Precise Construction of Micro-structures and Porous Geometry via Functional Composition[END_REF] and democratized since then (see for instance [START_REF] Antolin | Optimizing Micro-Tiles in Micro-Structures as a Design Paradigm[END_REF][START_REF] Bouclier | IGA: Non-conforming Coupling and Shape Optimization of Complex Multipatch Structures[END_REF][START_REF] Hirschler | Fast and multiscale formation of isogeometric matrices of microstructured geometric models[END_REF]). The idea is to model the two scales of the lattice structure separately first, and then to build the full detailed geometric model by composing the obtained micro and macro models. More specifically, within this approach, we consider that the reference microstructure is defined by a spline-based model 𝑇 , i.e., the domain Ω ref in Figure 2 is the image of the mapping 𝑇 :

Ω ref = 𝑇 [ Ω], (50)
where Ω is the parameter space of 𝑇 . For analysis purposes, we only consider analysis-suitable spline model made of multiple patches with matching interfaces. The spline basis functions of 𝑇 are used to defined the finite element spaces (the 𝑁 𝑖 in Equation ( 6)). Additionally, the mappings 𝐺 (𝑠) are given as Bezier multivariates, which might be obtained via a Bezier decomposition of a spline mesh. Finally, we combine these two mappings to build an explicit representation of the lattice cells. This is done via functional composition, meaning that the subdomains are the images of the composed mappings:

Ω (𝑠) = (𝐺 (𝑠) ∘ 𝑇 )[ Ω]. (51)
More information regarding this multiscale geometric representation can be found in [START_REF] Elber | Precise Construction of Micro-structures and Porous Geometry via Functional Composition[END_REF][START_REF] Hirschler | Fast and multiscale formation of isogeometric matrices of microstructured geometric models[END_REF] for instance.

Remark 8. Despite the fact that we leverage the use of spline functions for both, the geometry description and the solution discretization, following the isogeometric paradigm, the method presented is completely agnostic to type of finite element discretization applied to the PDE.

A multiscale Greedy procedure

Assembly via lookup tables

As quickly mentioned previously, in order to build the reduced basis of the local stiffness matrices as introduced in Equation ( 38), we employ a greedy approach. This enables to get a reduced basis that interpolates some local stiffness matrices, i.e., those forming the principal sub-problems. The crucial point here is that we do not build the reduced basis on the local stiffness matrices directly. This might be too expensive. Instead, the similarities between the cells can be observed by looking at the weak formulations given in Equation ( 5). After a change of coordinates by pulling back the macro-mappings 𝐺 (𝑠) , all local bilinear forms are defined over the reference domain Ω ref , as: where Ĉ (𝑠) ∶ Ω ref → ℝ 3×3×3×3 are fourth-order tensor fields which encapsulate all the dependencies on the macromappings. I.e., it encodes the material properties of the tensors C (𝑠) as well as all the quantities depending on the macro-mappings 𝐺 (𝑠) . So, one can take a look to these fields in order to build the reduced basis, instead of doing it over the local stiffness matrices. Furthermore, as all the other operators and quantities in Equation ( 52) are the same for every cell, one may resort to the approach developed in [START_REF] Hirschler | Fast and multiscale formation of isogeometric matrices of microstructured geometric models[END_REF]: It is possible to build lookup tables with precomputed integrals in order to assemble quickly local stiffness matrices. This approach consists in introducing suitable polynomial approximations of the tensor fields:

a (𝑠) (𝑢, 𝑣) = ∫ Ω ref ε (𝑢) ∶ Ĉ (𝑠) ∶ ε (𝑣) dξ, (52)
Ĉ 𝐻(𝑠) (ξ) = 𝐍 ⊤ (ξ)𝐀 (𝑠) , i.e., Ĉ 𝐻(𝑠) 𝑖𝑗𝑘𝑙 (ξ) = 𝑛 𝐴 ∑ 𝑞=1 𝑁 𝑞 (ξ)𝐴 (𝑠) 𝑖𝑗𝑘𝑙,𝑞 , (53)
where 𝐍 contains 𝑛 𝐴 polynomial basis functions and 𝐀 (𝑠) collects the associated polynomial coefficients. The latter ones are constants and encapsulate the dependencies on the macro-mappings 𝐺 (𝑠) and the material properties. Consequently, one can take them out of the integrals when substituting Equation (53) into Equation [START_REF] Parvizian | Finite cell method[END_REF], which leads to integrals that are independent of the macro-mappings and shared by all the cells. Those integrals are defined over the reference domain Ω ref , computed once and for all, and stored in lookup tables 𝐓 ref .

Building the stiffness matrices is then done by combining those precomputed integrals together with the polynomial coefficients, i.e., the assembly reads as:

𝐊 (𝑠) = 𝐓 ref ⊗ 𝐀 (𝑠) , (54)
where the operator ⊗ encodes multiplications, summations, and reordering of both quantities. We refer again to [START_REF] Hirschler | Fast and multiscale formation of isogeometric matrices of microstructured geometric models[END_REF] where this multiscale assembly strategy is detailed.

Remark 9.

It should be noticed that applying the fast assembly strategy via lookup tables [START_REF] Hirschler | Fast and multiscale formation of isogeometric matrices of microstructured geometric models[END_REF] is just a choice. The presented solver does not depend on this method. It however offers additional advantages as being able to perform matrix-vector products between the 𝐊 (𝑠) and input vectors in a matrix-free fashion (i.e., without explicitly building the matrices, but by using the lookup table and the polynomial coefficients directly, to save memory for instance). As will be presented in Section 4.2.2, it also simplifies the construction of the reduced basis of the stiffness matrices as one can operate on a small collections of coefficients 𝐀 (𝑠) instead of functions Ĉ (𝑠) (ξ) or large data 𝐊 (𝑠) . In this work, the procedure to identify the principal subdomains relies on a greedy algorithm (which can be viewed as a truncated Gram-Schmidt process), which is a common technique to build reduced bases, see for instance [START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF][START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF]. As discussed in Remark 9, the snapshot matrix could be formed using either the stiffness matrices 𝐊 (𝑠) , the functions Ĉ (𝑠) (ξ) in the bilinear forms (see Equation ( 52)), or the polynomial coefficients 𝐀 (𝑠) as defined by Equation [START_REF] Praster | An adaptive FE2 approach for fiber-matrix composites[END_REF]. Using directly the stiffness matrices is obviously the more expensive as it contains a lot of data. It is more efficient to consider one of the two other options. Here, we work directly with the polynomial coefficients (that come from the polynomial approximations of the fourth-order tensor fields Ĉ (𝑠) ). The greedy process for selecting the principal subdomains and performing the affine decomposition is given in Algorithm 4, and illustrated in Figure 5. Interestingly, due to its construction, the reduced space interpolates the principal subdomains, i.e.:

Å (𝜎 𝑖 ) = 𝐀 (𝜎 𝑖 ) , ∀𝑖 = 1, … , 𝑁 𝔯𝔟 . (55)
Consequently, we can perform a change of basis from the orthonormal basis ℤ formed during the greedy process to the basis formed by quantities defined over the principal subdomains 𝔸 = [𝐀 (𝜎 1 ) , … , 𝐀 (𝜎 𝑁 𝔯𝔟 ) ], i.e.:

Å (𝑠) = ‖𝐀 (𝑠) ‖ 2 ℤ𝛃 (𝑠) = 𝔸𝛂 (𝑠) . (56)
As a direct consequence, we end up with the approximate affine decompositions of the stiffness matrices as given in Equation ( 38) which is the starting point of our inexact FETI-DP solver. Instead of using the polynomial coefficients 𝐀 (𝑠) , one could perform the affine expansion of the functions Ĉ (𝑠) via the Empirical Interpolation Method. It consists in a similar greedy process as the one defined in Algorithm 4, but it involves an additional step where so-called magic points are selected. The interested readers may refer to [START_REF] Hesthaven | Certified Reduced Basis Methods for Parametrized Partial Differential Equations[END_REF][START_REF] Quarteroni | Reduced Basis Methods for Partial Differential Equations[END_REF] for more details regarding EIM and reduced basis space construction in general.

Numerical examples

All the following numerical examples were computed on a 2.50GHz i7-11850H (11th Gen Intel(R) Core(TM)) processor with 32 GB RAM (serial computing only). More precisely, three 2D examples and two 3D examples are investigated in this section. They all share same material parameters, i.e., an isotropic elastic material with Young's modulus 𝐸 = 5000 MPa and Poisson's ratio 𝜈 = 0.40. Regarding the solver, the global convergence relative tolerance was set to 𝑡𝑜𝑙 GMRES = 1𝑒 -5, while the local convergence relative tolerance for the interface problem was set to 𝑡𝑜𝑙 CG = 1𝑒 -11. The implementation has been done with Python using numpy/scipy packages [START_REF] Harris | Array programming with NumPy[END_REF][START_REF] Virtanen | SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF]. It also calls Cholesky factorization routines provided by scikit-sparse package which makes use of CHOLMOD from the SuiteSparse linear algebra package [START_REF] Chen | Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate[END_REF].

2D examples

Description

Let us start with the two-dimensional examples. The first inputs were the macro-geometries which are depicted in Figure 6. All these geometries were geometrically described via spline based models. They are of three kinds:

1. Fig. 6(a) -A rectangular domain subjected to an imposed displacement. The geometry was described with a linear B-Spline patch.

2. Fig. 6(b) -A curved beam under bending. The geometry was (exactly) described with a quadratic NURBS patch.

3. Fig. 6(c) -A brake pedal where the loading case mimics the operating conditions. The geometry was also described with a quadratic NURBS patch with knot-vectors 𝑈 = {0, 0, 0, 0.5, 0.5, 1, 1, 1} and 𝑉 = {0, 0, 0, 1, 1, 1}.

The coordinates and weights of the control points are given in the mentioned figure.

u x = 0 u y = 0 u x = 0.5 u y = 0 x y L = 8 h = 4
(a) Then, we associated a reference cell (microscale) to each of these macro-geometries. Again, we modelled these with spline based models. Two cells were considered as depicted in Figure 7:

u x = 0 u y = 0 t x = -1 x y R = 50 w = 50 (b) u = 0 u = 0 x y P =
1. Fig. 7(a) -A hollow square linked by a cross. It required twenty-four linear Bézier patches to build an analysissuitable model (i.e., matching interfaces between the patches, and no trimmed patches, see, e.g., [START_REF] Bouclier | IGA: Non-conforming Coupling and Shape Optimization of Complex Multipatch Structures[END_REF]). The model was further refined to generate the analysis model. This analysis discretization was made of 384 cubic elements, which led to 1866 DOF per cell.

2. Fig. 7(b) -An auxetic lattice pattern (obtained through a topology optimization procedure in [START_REF] Wang | Design of materials with prescribed nonlinear properties[END_REF]). We modelled the geometry with twenty quadratic matching NURBS patches. The refined version used during the analysis was made of 832 cubic elements for a total of 4368 DOF per cell.

We combined the crossed hollow square unit cell (Figure 7 Several results from these structural analyses are shown in Figure 8. We depict the displacement field and the von Mises stress field for the three test cases. As we perform full-scale simulations, we directly have access to these quantities over the different cell-struts.

Solver performance

We study the performance of the developed solver on these 2D numerical examples. A first question concerns the influence of the reduced bases. For instance, how many principal cells is it suitable to pick? What is the influence of the size of the reduced bases on the iterative convergence of the solver? Size of the reduced bases. Figure 9 gives an insight of the behavior of the solver w.r.t. the number of principal cells, (i.e., by varying on the tolerance 𝑡𝑜𝑙 RB during the greedy process, see Table 4). We performed this study on the curved beam problem (Figure 8(b)). More precisely, the corresponding results read as follows:

• 27 iterations when 6 principal cells out of the 512 cells were selected (𝑡𝑜𝑙 RB = 1𝑒 -1), Figure 9: Influence of the reduced basis on the convergence of the solver. The larger the reduced basis, the lower the number of iterations the solver needs.

• 23 iterations when 8 principal cells out of the 512 cells were selected (𝑡𝑜𝑙 RB = 1𝑒 -2),

• 8 iterations when 11 principal cells out of the 512 cells were selected (𝑡𝑜𝑙 RB = 1𝑒 -3),

• 6 iterations when 14 principal cells out of the 512 cells were selected (𝑡𝑜𝑙 RB = 1𝑒 -4),

• 4 iterations when 20 principal cells out of the 512 cells were selected (𝑡𝑜𝑙 RB = 1𝑒 -5),

• 4 iterations when 21 principal cells out of the 512 cells were selected (𝑡𝑜𝑙 RB = 1𝑒 -6),

• 3 iterations when 27 principal cells out of the 512 cells were selected (𝑡𝑜𝑙 RB = 1𝑒 -7).

This observed behavior was to be expected: the larger the RB (i.e., the lower 𝑡𝑜𝑙 RB ), the faster the convergence (in terms of number of global iterations). In terms of computational cost, it is straightforward however that increasing the number of principal cells leads to more memory usage. Indeed, it leads to compute and to store more quantities, as for instance, the Cholesky factorizations of the principal local problems. Regarding the computational time, selecting a smaller number of principal cells decreases the time needed to build the ROMs, and the time per global iteration. However, as it leads to a larger number of global iterations, the overall computational time might be above the one obtained with a larger number of principal cells. And this is somehow what we observed in the curved beam example. Indeed, the overall computational time was:

• around 25 s for the cases with 6 and 8 principal cells (respectively, 𝑡𝑜𝑙 RB = 1𝑒 -1 and 𝑡𝑜𝑙 RB = 1𝑒 -2),

• around 11 s for all the other cases, with a small decrease while increasing the RB size (i.e., lower tolerance while building the RBs in the greedy process).

We experienced the same behavior for other test cases: A larger number of principal cells leads to fewer global iterations, it however increases (a bit) the memory usage, but it leads to equivalent or less overall computational time. Thus, we advise to set a relative tolerance in the greedy selection algorithm (see again Table 4) between 1𝑒-5 and 1𝑒-7. This tolerance might seem low, but one should keep in mind that this tolerance is set over the local stiffness matrices only, and not over the local DD operators directly (e.g., the primal and dual Schur complements, see again Equation ( 39) and Remark 2). The particularity here is that we cannot build the ROMs a priori, and even more, we do not have access to the quantities for which we are building ROM-based approximations. For instance, we cannot run a similar greedy process on the Schur complements, directly, as it would mean to build them all, which would take more time than simply solving the problem with a standard FETI-DP approach. That is why, by selecting a sufficiently large number of principal cells, we ensure that the ROMs of the local DD operators are suitable.

Remark 10. The macro-geometry of the curved beam is modelled with a NURBS, see Figure 6(b). Hence, the quarter circles are exactly represented but not parameterized in the sense of the arc length. Thus, the cells do not follow exactly a circular repeating pattern. Furthermore, we do not apply a change of basis to rotate the cells prior the identification of the principal cells. Using an algebraic macro-mapping together with a change of basis could reduce the number of principal cells for this particular example. Scalability. In order to further investigate the performance of our ROM-based inexact FETI solver, we performed a scalability test for the three 2D examples. The results are presented in Table 1. The table also provides the results obtained with a standard FETI-DP solver, which means that the interface problem is "exactly" formed. In practice, it means that all the local stiffness matrices have been factorized prior to running the iterative algorithm. This is why, in the Memory column for the Standard FETI-DP data in Table 1, one can notice a substantial amount of memory usage in comparison with the Memory usage involved when employing the developed approach where we limit the construction of local operators to the picked principal cells. This is a great advantage of the proposed approach. Despite the use of full-scale models of lattice structures, the memory requirements were limited as we exploit the inherent repetitive pattern of those structures. On the contrary, black-box solvers are not aware of the structure of the problem at hand, and, hence, the available computational resources were rapidly surpassed by these (standard) approaches. This point is supported by the results: We were not able to compute the most refined version of the auxetic rectangular problem (64-by-32 cells for a total of 8.9M) with the standard FETI-DP approach and our computational resources, while we easily succeeded with the developed ROM-based FETI-DP approach.

The number of iterations required to solve the interface problems was the same with both approaches (Loc. It. columns). This was expected as the FETI-DP is known to be scalable for 2D problems when the so-called corners are chosen as the primal variables, see for instance [START_REF] Klawonn | Dual-primal feti methods for linear elasticity[END_REF]. The increase of the number of total cells has little influence on the number of global iterations for the ROM-based solver, as observed, for instance, for the brake pedal problem. However, this is minor as the solver still converged in very few global iterations (5 in the finest configuration). Furthermore, Table 1 provides in the columns Comp. Time the computational time in the following format:

• total time (time during preprocessing step, time for the iterative algorithm).

For the proposed ROM-based FETI-DP, the preprocessing step consists in building the reduced models (⋅) (𝑠) that approximate the DD local operators (see, e.g., Equation ( 43)), whereas in the standard FETI-DP, it consists mainly in factorizing all the local operators. Our test cases have an increasing level of complexity, in the sense of the geometry. For the rectangular problem, all the cells were actually the same. For the curved beam, they exhibited rather small geometric differences, whereas for the brake pedal, each cell had a different shape. Regarding the first example, the presented approach shows significant speedup in comparison with the standard approach. Even in the finest configuration with 8.9M DOF, the solver took only 5 s to get the solution.

Remark 11. The outstanding solver performance for the rectangular domain example is due to the fact that all the cells are geometrically identical. In this particular case, the ROM introduces no approximations, being the reduced bases composed of a single basis entry and all DD operators exactly recovered. Consequently, only one single global iteration is required. Indeed, the role of the global iterations in our solver is to compensate the ROM induced errors. If no approximations are introduced by the ROM, as it is the case for this particular example, then the inexact FETI-DP preconditioner is nothing else than the inverse of the system matrix, i.e., Å-1 = 𝐀 -1 .

For the curved beam, our solver took half the time needed with the standard FETI-DP, and for the brake pedal, it took a similar amount of time. For this last case, the number of principal cells was about 30 which explain the increase in computational time. However, as said previously, the memory usage was much less with the ROM-based approach. Applying blindly standard FETI-DP approach for large lattice structure would be unreasonable for this reason. This will become even clearer in the following 3D examples. We also think that some speedup with the ROM-based solver could be further obtained. Indeed, the code is partly built upon the sparse library from scipy which does not provide top performing matrix-matrix and matrix-vector products (contrary to the dense linear algebra which makes use of BLAS). Approximating the local solutions via the RB-based process described in Equation ( 46) took the major part of the time.

3D examples

Description

We now study two 3D examples which are depicted in Figure 10. The first test case consists in a straight beam under bending. Due to the symmetry of the problem, we considered only half of the problem (we could even consider only one quarter of the problem). The macro-geometry was modeled with a linear B-Spline patch with one single element. We associated to this macro-geometry a reference lattice cell. In this case, it was a Body-Centered Cubic cell. This BCC cell was modeled with 32 linear B-Spline patches. The geometric model was further refined for analysis purposes. The analysis discretization associated to one cell included 832 quadratic spline elements for a total of 9021 DOF. Several results from a structural analysis are given in Figure 11. Again, as we solved the full-scale problem, we had direct access to the displacement and the stress field over the different cell-struts.

The second test case (see again Figure 10(b)) consists in a curved beam which was subjected to bending too. The macro-geometry was modeled with a quadratic NURBS with one element. This enabled to exactly represent the circular shape of the beam. We associated to this macro-geometry a reference lattice cell whose shape looks like a cross with curved and circular arms. Initially, the geometric model of the reference cell was made of 7 cubic patches which were then further refined such that the analysis model included 56 cubic elements, and a total of 2175 DOF. We provide an illustration of the displacement and stress field for this test case in Figure 12.

Solver performances

We performed similar experiments as for the 2D examples. More specifically, we study the behavior of the presented ROM-based FETI-DP solver while increasing the number of cells for each problem. The results are summarized in Table 2, that presents the same structure as Table 1.

As already mentioned, an important characteristic of the proposed solver is the drastic reduction of the memory usage in comparison to "black box" solvers, as for instance a standard FETI-DP solver. This has been seen on the 2D examples, and it is even more obvious for the 3D examples. Here, only the coarser configurations could be solved with the standard method. Indeed, the memory limit was rapidly reached with the classic FETI-DP, whereas the developed ROM-based FETI-DP solver enabled us to simulate geometries with a higher number of lattice cells and low memory resources. Regarding the results presented in Table 2, enormous speedup was obtained for the straight beam problem. The reasons are exactly the same as those mentioned for the first of the 2D examples, i.e., the rectangle with the auxetic lattice pattern. Here, there were, by construction, no approximations of the local operators, and the reduced bases included one single element only as all the cells were the same. Thus, the solver took less than one minute to solve a problem with 6 × 6 × 24 cells and 7.8M DOF. Regarding the curved beam problem, it required 3 or 4 global iterations to converge. It is faster than the standard FETI-DP, and required much less memory. Thus, it was possible to run on a simple laptop a case with 6 × 6 × 56 cells and 3.7M DOF with the ROM-based FETI-DP solver. Indeed, only 20 principal cells were selected out of the 2016 cells. On the contrary, the standard approach reached the memory limit much before, and only a case with 4 × 4 × 32 cells and 1.1M DOF was trackable.

Conclusions

In this work, we developed a novel method for performing full fine-scale analysis of architected materials at very low computational cost. This dedicated solver takes advantage of the intrinsic geometrical nature of such structures; that is, the non-overlapping domain decomposition into cells that are all similar. The faced challenge was: Is it possible to leverage on those similarities to save memory usage and computational time during the simulation? Resorting to the class of DD methods, especially to the inexact variants, and complementing the approach with ROM, this turned out possible. More precisely, we built upon inexact FETI-DP algorithms that offer the possibility of not solving exactly all the local, cell-wise systems. Then, to design an efficient block preconditioner, we proposed to extract the principal cell-wise stiffness matrices with a greedy approach, and to use them to approximate all the local quantities occurring in a standard FETI-DP algorithm. This results in a solver which only needs, in terms of local solutions, the ones involving the principal stiffness matrices. Therefore, numerous local factorizations are saved compared to standard DD algorithms. Furthermore, although not mandatory in our methodology, we adopted a geometrical modeling based on spline composition and performed, following the isogeometric paradigm, analysis of lattices using look-up tables and macro-fields (see [START_REF] Hirschler | Fast and multiscale formation of isogeometric matrices of microstructured geometric models[END_REF]). This allowed us to extract the principal local operators by acting only on the macro-fields, which avoids forming all the local operators, and therefore makes our method also multi-scale and matrix-free.

Several numerical experiments in 2D and 3D and covering different micro-and macro-mappings were performed to assess the solver performance. Only serial computations were conducted at this stage. Comparisons with the standard FETI-DP were also carried out to highlight the interest of our solver. In any case, the memory usage was drastically reduced with our solver compared to the standard FETI-DP, which allowed us to compute problems of several millions of DOF within few minutes using an off-the-shelf laptop. Eventually, the computational time was also very low, in particular when the macro-mapping is affine (since all the cells become identical), which is usually the case in real architected material designs. The obtained results suggest that introducing multiscale modeling for lattice structures to keep computing resources reasonable could be circumvented with such HPC algorithms. In some sense, the multiscale, heterogeneous, and periodicity features of the structure are used directly in the solver, so we do not need to bring them at the modeling stage. In this respect, an efficient parallel implementation, and the extension of the method to the nonlinear regime, as well as to take into account manufacturing-induced defects, surely constitute interesting prospects to this work.

Figure 2 :

 2 Figure 2: Problem to be solved along with its domain decomposition. The finite element spaces associated to the different subdomains are constructed from the basis functions of the reference unit-cell.

Figure 3 :

 3 Figure 3: Illustration of the different solution spaces for the DD algorithms when applied to lattice structures. In our case of a FETI-DP based algorithm, the solution is sought in the intermediate space W ℎ , where the continuity of the local displacements is enforced at the corners of the subdomains.

Figure 4 :

 4 Figure 4: Example of results obtained with the reduced basis approach. 6 principal subdomains (marked as 𝜎 𝑘 , 𝑘 = 1, … , 6) and associated local stiffness matrices are selected (see top-left). Then, all the local stiffness matrices can be expressed as a linear combination of the principal ones.The exact values of the associated RB coordinates are given for subdomains 𝑠 1 and 𝑠 2 in the table at the bottom-left. For completeness, plots of the value of each RB coordinate for all subdomains are realized on the right.

  2.4), i.e., in:

(

6 Figure 5 :

 65 Figure 5: Iterative identification of the principal subdomains via the greedy procedure occurring at the macro-scale. The ROM Error corresponds to the norm of the residuals ‖𝚫 (𝑠) 𝑖 ‖ ∞ as written in Algorithm 4.

Algorithm 4

 4 Greedy algorithm for the identification of the principal subdomains. 𝛇 1 … 𝛇 𝑁 𝔯𝔟 ] Orthonormal basis S = {𝜎 1 , … , 𝜎 𝑁 𝔯𝔟 } Principal subdomain indices Å (𝑠) = ‖𝐀 (𝑠) ‖ 2 ℤ𝛃 (𝑠) , ∀𝑠 Certified affine decompositions 4.2.2. Algorithm to identify the principal subdomains

𝑤 4 Figure 6 :

 46 Figure 6: The macro-geometries and the loading scenarios that define the 2D examples: (a) a rectangular domain subjected to an imposed displacement, (b) a beam under bending, and (c) brake pedal subjected to operating conditions.

Figure 7 :

 7 Figure 7: The 2D cells used during the numerical studies: (a) a hollow square with a cross along its diagonal, and (b) an auxetic lattice pattern.

  (a)) with the curved beam problem (Figure 6(b)) and the brake pedal problem (Figure 6(c)). The auxetic lattice (Figure 7(b)) was combined with the rectangular macro-domain subjected to an imposed displacement (Figure 6(a)).

Figure 8 :

 8 Figure 8: Results of the structural analyses for the 2D test cases. Note that the auxetic behavior can be observed in (a): The initial macro rectangle inflates while it was subjected to traction.

Figure 10 :

 10 Figure 10: The macro-geometries, the loading scenarios, and the reference cells that define the 3D examples: (a) a straight beam filled by Body-Centered Cubic patterns, and (b) a curved beam filled by patterns initiated from a reference cell made of a cross with circular arms.

Figure 11 :

 11 Figure 11: Structural analysis results for the straight 3D beam problem.

Figure 12 :

 12 Figure 12: Structural analysis results for the curved 3D beam problem.

:

  𝐯 = ( 𝐯 𝑅 𝐯 𝑃 ), with 𝐯 𝑅 ∈ ℝ 𝑛 𝑅 , 𝐯 𝑃 ∈ ℝ 𝑛 𝑃 .

	v 𝑃	← 𝐯 𝑃 -Ů ⊤ 𝑅𝑃 𝐯 𝑅	Form the right-hand side of the coarse problem
	𝐳 𝑃 ←	S -1 𝑃 𝑃 v 𝑃	Solve the coarse problem (see Section 3.2.1)
	𝐳 ←	Ů ⊤ 𝑅𝐷 𝐯 𝑅 -𝐁 𝑅 Ů 𝑅𝑃 𝐳 𝑃	Approximate the interface vector
	Ouput vector:
	Out: 𝐳 = Ů 𝐯, 𝐳 ∈ ℝ 𝐿 .

Table 1 :

 1 Performance of the solver in the case of 2D problems.

	Test case			ROM-based inexact FETI-DP			Standard FETI-DP
	Macro-geom. #Cells	#DOFs	Memory Comp. Time Glo. It. Loc. It.	Memory	Comp. Time	Loc. It.
	Rectangle	16 × 08	0.5M	0.0 GiB 0.3s (0s, 0s)	1	28	1.6 GiB	8s ( 4s, 4s)	28
	Rectangle	32 × 16	2.2M	0.1 GiB	1s (0s, 1s)	1	29	6.4 GiB 38s (25s, 14s)	29
	Rectangle	64 × 32	8.9M	0.2 GiB	5s (1s, 4s)	1	29	> available	-	-
	Curved beam	16 × 08	0.2M	0.2 GiB	5s (3s, 2s)	3	35	0.8 GiB	6s ( 4s, 2s)	35
	Curved beam	32 × 16	1.0M	0.3 GiB 11s (3s, 8s)	3	35	3.2 GiB 22s (13s, 9s)	35
	Curved beam	64 × 32	3.8M	0.5 GiB 38s (3s, 35s)	3	35	13.4 GiB 73s (38s, 35s)	34
	Brake pedal	32 × 04	0.2M	0.3 GiB 10s (6s, 4s)	3	42	0.8 GiB	6s ( 3s, 3s)	42
	Brake pedal	64 × 08	1.0M	0.3 GiB 21s (6s, 15s)	4	43	3.2 GiB 25s (14s, 11s)	42
	Brake pedal	128 × 16	3.8M	0.5 GiB 86s (5s, 81s)	5	43	12.7 GiB 80s (38s, 42s)	43

Table 2 :

 2 Performance of the solver in the case of 3D problems.

	Test case		ROM-based inexact FETI-DP			Standard FETI-DP
	Macro-geom. #Cells	#DOFs	Memory Comp. Time Glo. It. Loc. It.	Memory	Comp. Time Loc. It.
	Straight beam 2 × 2 × 08	0.3M	0.5 GiB	0min 3s	1	103	4.2 GiB	0min 36s	103
				(1s, 2s)				(10s, 26s)
	Straight beam 3 × 3 × 12	1.0M	0.9 GiB	0min 7s	1	106	14.0 GiB	2min 4s	106
				(2s, 5s)				(31s, 93s)
	Straight beam 4 × 4 × 16	2.3M	1.7 GiB	0min 14s	1	108	> available	-	-
				(3s, 11s)
	Straight beam 5 × 5 × 20	4.5M	3.1 GiB	0min 26s	1	109	> available	-	-
				(4s, 22s)
	Straight beam 6 × 6 × 24	7.8M	5.1 GiB	0min 44s	1	110	> available	-	-
				(6s, 38s)
	Curved beam 3 × 3 × 24	0.5M	1.2 GiB	0min 36s	4	93	6.9 GiB	0min 49s	92
				(8s, 24s)				(23s, 26s)
	Curved beam 4 × 4 × 32	1.1M	1.6 GiB	1min 0s	3	103	15.7 GiB	1min 53s	102
				(12s, 48s)				(49s, 64s)
	Curved beam 5 × 5 × 40	2.2M	2.2 GiB	1min 57s	3	113	> available	-	-
				(19s, 98s)
	Curved beam 6 × 6 × 56	3.7M	3.1 GiB	3min 40s	4	120	> available	-	-
				(20s, 230s)

Acknowledgment. This research was supported by the European Union Horizon 2020 research and innovation program, under grant agreement No 862025 (ADAM2) and by the Swiss National Science Foundation through the project No 40B2-0_187094 (BRIDGE Discovery 2019). The second author, Robin Bouclier, also gratefully acknowledges support from the French National Research Agency under Grant ANR-22-CE46-0007 (AVATAR).