# Acknowledgments

The research leading to these results has received funding from, on the one hand, the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 320815, Advanced Grant Project COMP-DES-MAT, and, on the other hand, the Spanish Ministry of Science and Innovation under grant BIA2011-24258.

# Abstract

The present work is concerned with the application of projection-based, model reduction techniques to the efficient solution of the cell equilibrium equation appearing in (otherwise prohibitively costly) two-scale, computational homogenization problems. The main original elements of the proposed Reduced-Order Model (ROM) are fundamentally three. Firstly, the reduced set of empirical, globally-supported shape functions are constructed from pre-computed Finite Element (FE) snapshots by applying, rather than the standard Proper Orthogonal Decomposition (POD), a partitioned version of the POD that accounts for the elastic/inelastic character of the solution. Secondly, we show that, for purposes of fast evaluation of the nonaffine term (in this case, the stresses), the widely adopted approach of replacing such a term by a low-dimensional interpolant constructed from POD modes, obtained, in turn, from FE snapshots, leads invariably to ill-posed formulations. To safely avoid this ill-posedness, we propose a method that consists in expanding the approximation space for the interpolant so that it embraces also the gradient of the global shape functions. A direct consequence of such an expansion is that the spectral properties of the Jacobian matrix of the governing equation becomes affected by the number and particular placement of sampling points used in the interpolation. The third innovative ingredient of the present work is a points selection algorithm that does acknowledge this peculiarity and chooses the sampling points guided, not only by accuracy requirements, but also by stability considerations. The efficiency of the proposed approach is critically assessed in the solution of the cell problem corresponding to a highly complex porous metal material under plane strain conditions. Results obtained convincingly show that the computational complexity of the proposed ROM is virtually independent of the size and geometrical complexity of the considered representative volume, and this affords gains in performance with respect to finite element analyses of above three orders of magnitude without significantly sacrificing accuracy –-hence the appellation High-Performance ROM.

# 1 Introduction

## 1.1 Motivation

### 1.1.1 Moore's and Parkinson's laws in computational modeling

A basic precept in constructing a successful computational model –-one that strikes the right balance between accuracy and simplicity–- is, quoting M.Ashby [1], to unashamedly distort the inessentials in order to capture the features that really matter. Yet with the general availability of fast computers with large memories, and the exponential growth presaged by Moore's law of such capacities, it is becoming increasingly difficult to resist the temptation to violate this basic precept and, in the quest for higher accuracy, indiscriminately include features that do not contribute significantly to the overall response of the modeled system. As pointed out by Venkataraman et al.[2], this had led to the paradox that certain engineering problems that were computationally tackled for the first time 40 years ago appear to remain comparatively costly, even though the capacity of computers have increased one trillion-fold during this period. This apparent paradox is of course not exclusive of computational physical modeling, but rather just another manifestation of the general adage “software applications will invariably grow to fill up increased computer memory, processing capabilities and storage space”, known as the computerized Parkinson's law [3].

### 1.1.2 The two-scale homogenization problem

An engineering research area that seems to be falling into this paradigmatic trend is, as we intend to argue in the ensuing discussion, the two-scale modeling, via homogenization, of materials such as composites and polycrystalline metals, characterized by exhibiting a clear heterogeneous composition at some lower length scale –-the micro-, or meso-, scale–-, but that can be regarded as homogeneous at the length scale at which engineering predictions are needed –-the macro-scale. The actual challenge in the macro-scale continuum description of the mechanical behavior of such materials lies in the determination of a constitutive connection between macro-stresses and macro-strains that accurately reflects the material properties and geometrical arrangement of the distinct phases at the micro-/meso-scale. Under the hypotheses of either periodicity or statistical homogeneity, on the one hand; and scale separation, on the other hand, it is well-known [4] that this constitutive link can be systematically established by solving, for each point at the coarse scale, a boundary value problem (BVP) on a certain representative microscopic subdomain (the cell equilibrium problem). In a strain-driven formulation of this BVP, the macro-strain at a given point acts as “loading parameter”, in the form of appropriate essential boundary conditions, whereas the associated macro-stress is obtained through volume averaging –-i.e., homogenization–- of the corresponding micro-stress field.

### 1.1.3 Evolution of homogenization approaches

When the discipline of continuum micromechanics began to flourish in the 60's and 70's of the last century, research focus was fundamentally directed, absent powerful computational tools (if any), towards the development of approximated, closed-form solutions of this BVP for certain types of geometrically and constitutively simple micro-structures. To arrive at these solutions, pioneers such as R. Hill [5], Z. Hashin [6], and K. Tanaka [7] struggled to identify and retain, guided sometimes by an uncanny physical intuition, only those features of the microstructural response that has a significant impact on the accuracy of coarse-scale predictions –-filtering out, thus, the “inessentials”. The advent of increasingly fast computers, and the concomitant advancement of numerical methods to solve differential equations, progressively fostered a shift towards the development of homogenization techniques that rely less on analytical results and more on numerical solutions, widening thereby their scope; for comprehensive surveys of these semi-analytical homogenization methods, the reader is referred to Refs. [8,9,10]. The approach termed by some authors [11,12] computational homogenization (or direct computational homogenization ) can be regarded as the culmination of this shift towards purely numerical methods: in this approach, the microscopic boundary value problem at each coarse-scale point is attacked using no other approximation than the spatial discretization of the pertinent solution strategy (finite element method, for instance), thus, circumventing the need for simplifying assumptions regarding the topological arrangement of the micro-phases and/or their collective constitutive behavior.

### 1.1.4 Infeasibility of direct computational homogenization methods

Although no doubt the most versatile and accurate homogenization technique, with no other limitation in scope than the imposed by the aforementioned hypotheses of statistical homogeneity and scale separation, the direct computational homogenization approach violates squarely the modeling precept outlined at the outset –-it does not discriminate between essential and irrelevant features in solving the fine-scale BPVs–-, making the accuracy/parsimony balance to tilt unduly towards the accuracy side and far from the parsimony one. The consequence is its enormous computational cost (in comparison with analytical and semi-analytical homogenization techniques). For instance, when the Finite Element (FE) method is the strategy of choice for both scales –-the commonly known as multilevel finite element method, abbreviated FE${\textstyle ^{2}}$ method [13]–-, one has to solve, at every increment and every Newton-Raphson iteration of the global, coarse-scale FE analysis, and for each Gauss point of the coarse-scale mesh, a local, non-linear finite element problem which, in turn, may involve several thousand of degrees of freedom. To put it metaphorically, in the FE${\textstyle ^{2}}$ method, the complexity of the overall analysis is submitted to the “tyranny of scales” [14,15,16]–- it depends on the resolution of both coarse- and fine- scale grids. This explains why, presently, almost five decades after R.Hill [17] laid the foundations of non-linear continuum micro-mechanics, and despite the dizzying speed of today's computers, the routine application of this general and versatile homogenization method to model the inelastic behavior of large structural systems featuring complex micro-/meso-structures is still considered intolerably costly, especially when the system has to be analyzed for various configurations, as in design optimization or inverse analysis.

### 1.1.5 Shortcomings of semi-analytical approaches

The foregoing review and critical appraisal of how things have evolved in the field of continuum micromechanics clearly teaches us that, to defeat the tyranny of scales, and develop a successful homogenization method, one should strive to adhere to Ashby's modeling principle and introduce –-as actually done in analytical and semi-analytical approaches–- appropriate simplifications in dealing with the fine-scale BVPs. Negating the need for such simplifications in the hope that the upcoming generations of peta-, exa- (and so forth) flops computers will eventually come to the rescue and cope with the resulting complexity, is, in the authors' opinion, a mistaken view, at odds with the true spirit of physical modeling, and the culprit for the paradoxical trend alluded to earlier.

Yet finding, for a given, arbitrarily complex microstructure, a suitable set of simplifying assumptions, let alone incorporating in a consistent manner such assumptions in the formulation of the local BVPs, is in general a formidably challenging endeavor. The admittedly brilliant idea underlying many advanced semi-analytical homogenization methods, such as the Transformation Field Analysis (TFA) [18] and variants thereof [19,20,21,22], of pre-computing certain characteristic operators (strain localization and influence tensors) by solving a carefully chosen battery of fine-scale BPVs has only partially relieved modelers from this burden: these methods are still predicated, to a lesser or greater extent, on ad-hoc assumptions connected with the constitutive description of the involved phases. Consideration of new materials with unstudied compositions will thereby require additional research efforts by specialists in the field and eventual modifications of the corresponding mathematical and numerical formulations –-in contrast to direct computational homogenization approaches, such as the FE${\textstyle ^{2}}$ method, in which the formulation is “material-independent”, and hence more versatile.

## 1.2 Goal

The current state of affairs in the field of two-scale homogenization seems to call, thus, for a unified homogenization approach that combines somewhat the advantages of direct computational homogenization and semi-analytical techniques. It would be desirable to have a homogenization method with a computational cost virtually independent of the geometric complexity of the considered representative volume, as in semi-analytical techniques –-i.e., one that defies the tyranny of scales. At the same time, it would be also interesting to arrive at a method whose mathematical formulation dispenses with ad-hoc, simplifying assumptions related with the composition of the heterogeneous material; i.e, one enjoying the versatility, unrestricted applicability and “user-friendliness” –-insofar as it would totally relieve the modeler from the often exceedingly difficult task of visualizing such assumptions –- of direct computational homogenization methods. The goal of the present work is to show that these desirable, apparently conflicting attributes can be conceivably achieved, for arbitrarily complex heterogeneous materials well into the inelastic range, by using the so-called [23] Reduced-Basis (RB) approximation in the solution of the cell BVPs.

## 1.3 Reduced-basis approach

### 1.3.1 Essence of the approach

Generally speaking, the reduced-basis approximation is simply a class of Galerkin approximation procedure that employs, as opposed to the FE method, but similarly to classical Rayleigh-Ritz solution techniques [24], globally supported basis functions. The main difference with respect to classical Rayleigh-Ritz schemes is that these basis functions or modes are not constructed from globally supported polynomials or transcendental functions (sines, cosines ...), but rather are determined from a larger set of previously computed, using the finite element (FE) method or other classical solution techniques, solutions of the BVP at appropriately selected values of the input of interest. These functions are commonly termed empirical basis functions [25], the qualifier empirical meaning “derived from computational experiments”.

As noted earlier, the input of interest or “loading” parameter in the fine-scale cell problem is the macro-scale strain tensor. Accordingly, the starting point for constructing the basis functions in the case under study would consist in solving a battery of cell BVPs for various, judiciously chosen macro-strain values. In the linear elastic regime, for instance, the displacement solution depends linearly on the prescribed macro-strain tensor, and hence it would suffice to perform 6 linear FE analysis (stretch in three directions and three shears); the corresponding modes would simply arise from orthonormalizing the displacement solutions of these problems. By constraining the cell to deform only into these 6 pre-defined, shape functions, one automatically obtains a genuine reduced-order model (ROM) of the cell. Note that the dimension of this ROM is totally independent of the spatial discretization (and, therefore, of the geometric complexity of the cell) used in the preliminary or offline FE analyses.

### 1.3.2 Dimensionality reduction

Things get a little bit more complicated in the inelastic range. The solution of the problem not only ceases to bear a linear relation to the prescribed macro-strain tensor: it also depends in general on the entire history of this coarse-scale kinematic variable. As a consequence, instead of single linear FE analyses, it becomes necessary to perform nonlinear FE studies on the cell subjected to various, representative macro-strain histories. The outcome of these FE calculations is a data set comprising an ensemble of hundred or even thousand (depending on the number of time steps into which the strain histories are discretized) displacement field solutions (also called snapshots). Were all these snapshots barely correlated with each other, the dimension of the manifold spanned by them would prove overly high, rendering the entire approach impractical –-it would no longer qualify as a truly reduced basis method. Fortunately, as we show in the present work, in general, most of these snapshots do display strong linear correlations between each other –-i.e., they have redundant information–-, and, in addition, contain deformation modes that are irrelevant to the quality of coarse-scale predictions –-in the sense that their associated stress fields have vanishing or negligible average volumes. Accordingly, all that is required to obtain a much lower dimensional representation of the solution data set, and therewith the desired reduced basis, is an automatic means to identify and remove this redundant and irrelevant information, while preserving, as much as possible, its essential features. This problem of removing unnecessary complexity from huge data sets so as to uncover dominant patterns is the central concern of disciplines such as digital image and video compression [26], and patter recognition [27], to name but a few, and thereby many efficient dimensionality reduction (or data compression, in more common parlance) algorithms already exist to deal with it. In the present work, we employ the arguably simplest and most popular of such algorithms: the Proper Orthogonal Decomposition (POD).

It becomes clear from the above discussion that, in the reduced basis approach, the inescapable task of discriminating what is essential and what is not is automatically carried out by these dimensionality reduction methods. In other words, the “burden” of simplification of the cell BVP in the RB approach is entirely borne by the computer, and not by the modeler, as it occurs in analytical and, to a lesser extent, in semi-analytical homogenization methods. It is precisely this feature that confers the advantages of versatility and “user-friendliness” alluded to earlier.

### 1.3.3 Numerical integration

Once the global shape functions have been computed, the next step in the construction of the reduced-order model of the cell is to introduce an efficient method for numerically evaluating the integrals appearing in the weak form of the cell BVP. Of course one can simply use the same Gauss quadrature formulae and the same sampling points (a total number of ${\textstyle n_{g}={\mathcal {O}}(n)}$, ${\textstyle n}$ being the number of mesh nodes) as the underlying finite element model. But this would be akin to integrating, say, a third-order polynomial function using thousand of sampling points–-a profligate waste of computational resources. Since displacement solutions for the cell BVP are constrained to lie in a reduced-order space of dimension ${\textstyle n_{u}<, it is reasonable to expect that the corresponding stresses, internal forces and Jacobians will also reside in reduced-order spaces of dimensions of order ${\textstyle {\mathcal {O}}(n_{u})}$, and consequently, only ${\textstyle p={\mathcal {O}}(n_{u})< sampling points would suffice in principle to accurately evaluate the corresponding integrals. The challenging questions that have to be confronted are where to locate these ${\textstyle p}$ sampling points and, loosely speaking, how to determine their associated weighting functions so that maximum accuracy in the integration is attained.

Approaches found in the model reduction literature that, directly or indirectly, deal with these fundamental questions can be broadly classified either as interpolatory approaches [28,29,30,31,32] or Gauss-type quadrature approaches [33,34], the former having a wider scope, for they also serve to reduce models derived from finite difference approximations [32,31]. The starting point for both is to collect, during the alluded to earlier finite element calculations, snapshots of the integrand or part of the integrand. In interpolatory approaches, the resulting ensemble of snapshots is submitted to a dimensionality reduction process, in the manner described previously for the solution basis functions, in order to compute a reduced set of dominant, orthogonal modes. Such orthogonal modes are employed to construct an interpolant of the integrand, using as interpolation points (the desired sampling points) those minimizing the interpolation error over the finite element snapshots. In the spirit of classical, interpolatory quadrature schemes such as Newton-Cotes [35], the resulting quadrature formula, and therefore, the weighting functions, emerges from approximating the integrand by this reduced-order interpolant. In Gauss-type quadrature procedures [33,34], by contrast, the selection of sampling points and the calculation of the accompanying weighting factors are simultaneously carried out, guided by a criterion of minimum integration error over the snapshots–-in a vein, in turn, similar to that used in classical Gauss-Legendre quadrature rules [35]. A more comprehensive review on both type of integration methods can be found in Ref. [36].

In the BVP under consideration, the output of interest is the volume average of the stresses over the cell domain and, therefore, accuracy is required not only in the integration of the equilibrium equation, but also on the approximation of the stresses themselves. This is the reason why, notwithstanding the presumably higher efficiency of Gauss-type quadrature (less integration error for the same number of sample points), attention is focused in the present work on interpolatory integration strategies, the variable subject to spatial interpolation being precisely the stresses.

## 1.4 Originality of this work

The idea of exploiting the synergistic combination of multiscale modeling and reduced basis approximation is admittedly not new. In the specific context of two-scale homogenization, it has been recently explored by Boyaval [37], Yvonnet et al. [38], and Monteiro et al. [39]. Traces of this idea can also be found in articles dealing with more general hierarchical multiscale techniques –-that do not presuppose either scale separation or periodicity/statistical homogeneity, or both–-, namely, in the multiscale finite element method [40,41,42] and in the heterogeneous multiscale method [43,44]. However, it should be noted that none of the above cited papers confronts the previously described, crucial question of how to efficiently integrate the resulting reduced-order equations, simply because, in most of them [37,40,41,42,43,44], integration is not an issue –- the fine-scale BVPs addressed in these works bear an affine relation with the corresponding coarse-scale, input parameter, as in linear elasticity, and, consequently, all integrals can be pre-computed, i.e., evaluated offline, with no impact in the online computational cost. Thus, the development of cell reduced-order models endowed with efficient, mesh-size independent integration schemes –-able to handle any material composition–- is a research area that, to the best of the authors' knowledge, still remains uncharted.

### 1.4.1 Main original contribution

The theory underlying ROMs that incorporate efficient interpolatory integration schemes, henceforth termed High-Performance ROMs (HP-ROMs), is still at its embryonic stage of development –-the first general proposal for parametrized BVPs dates back to 2004 [28]–- and many fundamental issues remain to be addressed. Foremost among these is the crucial question of well-posedness of the resulting system of algebraic equations: does the replacement of the integrand, or nonaffine term in the integrand, by a reduced-order interpolant always lead to a well-posed, discrete problem ? Examination of the reduced basis literature indicates that apparently no researcher has so far been confronted with ill-posed reduced-order equations, a fact that might certainly promote the view that uniqueness of solution can be taken for granted whenever the full-order model is well-posed. Unfortunately, this is not always so: we demonstrate in this work that the choice of the reduced-order space in which the interpolant of the integrand resides has a profound impact on the well-posedness of the discrete problem. In particular, we show that, in the case of the cell boundary-value problem, the widely adopted [29] approach of determining the basis functions for this space from (converged) FE snapshots leads invariably to ill-posed, discrete formulations. The main original contribution of the present work to the field of reduced-order modeling is the development of an interpolatory integration method that safely overcomes this type of ill-posedness. The gist of the method is to expand the interpolation space so that it embraces, aside from the span of the POD stress basis functions, the space generated –-and herein lies the novelty–- by the gradient of the (reduced-order) shape functions.

### 1.4.2 Other original contributions

An inevitable consequence of adopting the aforementioned expanded space approach is that, in contrast to the situation encountered when using standard interpolatory schemes in other parametrized BVPs [29], in the proposed method, the number and particular placement of sampling points within the integration domain influence notably the spectral properties (positive definiteness) of the Jacobian matrix of the governing equation, and therefore, the convergence characteristics of the accompanying Newton-Raphson solution algorithm. Another innovative ingredient of the present work is a points selection algorithm that does acknowledge this peculiarity and chooses the desired sampling points guided, not only by accuracy requirements (minimization of the interpolation error over the FE stress snapshot), but also by stability considerations.

Lastly, a further original contribution of the present work is the strategy for computing the global shape functions and stress basis functions. Instead of directly applying the POD over the snapshots to obtain the dominant modes, we first decompose these snapshots into (mutually orthogonal) elastic and inelastic components, and then apply separately the standard POD to the resulting elastic and inelastic snapshots. In so doing, the resulting reduced-order model is guaranteed to deliver linear elastic solutions with the same accuracy as the underlying (full-order) finite element model.

## 1.5 Organization of the document

The remaining of this document is organized as follows. Chapter 2is devoted to the formulation and finite element implementation of the cell equilibrium problem. Chapter 3, on the other hand, is concerned with the issue of the offline computation of the POD reduced basis. In Section 3.2, the Galerkin projection of the cell equilibrium equation onto the space spanned by the POD basis is presented. Chapter 4 outlines the procedure for efficiently integrating the reduced-order equilibrium equation. In Section 4.3, we discuss the crucial issue of where the low-dimensional approximation of the stress field should lie in order to obtain an accurate and at the same time well-posed reduced-order problem; in Section 4.3.3, the original proposal of expanding the stress basis with the gradient of the (reduced-order) shape functions is put forward. Section 4.4 delineates the derivation of the modal coefficients in the approximation of the stress field in terms of the stress values computed at the set of pre-specified sampling points. The determination of the optimal location of these sampling points is addressed in Section 4.5. For the reader's convenience and easy reference, in Section 4.6, both the offline and online steps leading to the proposed reduced-order model are conveniently summarized. Chapter 5 is dedicated to numerically assess the efficiency of the proposed model reduction strategy. Finally, in Chapter 6, some concluding remarks are presented.

In order to preserve the continuity of the presentation, details concerning the algorithmic implementation of the computation of the reduced basis and the selection of sampling points are relegated to the appendices.

# 2 First-order homogenization

## 2.1 Basic assumptions

The fundamental assumptions upon which the homogenization approach followed in this work rests are presented below. For a more in-depth description of the underlying axiomatic framework, the reader is referred to Refs. [45,46,47,48].

### 2.1.1 Existence of a representative subvolume

The homogenization approach employed in this work –-commonly known as first-order homogenization–- is only valid for materials that either display statistical homogeneity or spatial periodicity [45]. In both type of materials, it is possible to identify a subvolume ${\textstyle \Omega \subset \mathbb {R} ^{d}\,(d=2,3)}$, of characteristic length ${\textstyle l}$, that is representative, in a sense that will be properly defined later, of the heterogeneous material as a whole. Furthermore, this subvolume ${\textstyle \Omega }$ has to be small enough that it can approximately be regarded as a point at the coarse-scale level [4] (i.e, ${\textstyle l<, ${\textstyle l_{M}}$ being the characteristic length of the macro-continuum ${\textstyle \Omega _{M}}$, see Figure 1). This is the so-called scale separation hypothesis. In micro-structures that exhibit statistical homogeneity, this domain receives the name of Representative Volume Element (RVE), whereas in micro-structures that display periodicity, it is commonly known as repeating unit cell (RUC), or simply unit cell [45]. In the sequel, both the acronym RVE and the more generic term “cell” will be used interchangeably to refer to ${\textstyle \Omega }$.

 Figure 1: First-order homogenization.

### 2.1.2 Decomposition into macroscopic and microscopic contributions

The displacement ${\textstyle {\boldsymbol {u}}_{\mu }}$ at any point ${\textstyle \mathbf {x} \in \Omega }$ is assumed to be decomposed into macroscopic and microscopic parts; under the hypothesis of infinitesimal deformations, this decomposition can be written as:

 ${\displaystyle {\boldsymbol {u}}_{\mu }(\mathbf {x} )={\boldsymbol {\epsilon _{M}}}\,\mathbf {x} +{\boldsymbol {u}}(\mathbf {x} ),}$
(2.1)

where ${\textstyle {\boldsymbol {\epsilon _{M}}}}$ stands1 for the macroscopic strain tensor (the input parameter in the BVP we wish to efficiently solve) and ${\textstyle {\boldsymbol {u}}(\mathbf {x} )}$ denotes the so-called displacement fluctuation field (in turn, the basic unknown of this BVP). The macroscopic term ${\textstyle {\boldsymbol {\epsilon _{M}}}\,\mathbf {x} }$ represents the displacements that would have been observed had the material been homogeneous, whereas the fluctuating contribution ${\textstyle {\boldsymbol {u}}}$ accounts for the deviations from this homogeneous state due to the presence of heterogeneities [49].

The decomposition of the microscopic strain tensor ${\textstyle {\boldsymbol {\epsilon }}}$ follows from simply differentiating Eq.(2.1):

 ${\displaystyle {\boldsymbol {\epsilon }}(\mathbf {x} )={\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}(\mathbf {x} )}.}$
(2.2)

Notice that, in writing Eq.(2.2), one is tacitly assuming that the macroscopic strain tensor ${\displaystyle {\boldsymbol {\epsilon _{M}}}}$ is uniform over the spatial length associated to the cell size . This proviso renders the employed homogenization approach –-commonly termed first-order homogenization [12,48]–- not suitable for appropriately handling localization problems. For this reason, in the present work, we shall presuppose that localization of strains does not take place within the deforming cell.

Implicit in the scale separation assumption is the fact that fine-scale deformations only influence coarse-scale behavior through its volume average over the RVE; this implies that

 ${\displaystyle {\boldsymbol {\epsilon _{M}}}={\dfrac {1}{V}}\displaystyle \int _{\Omega }{\boldsymbol {\epsilon }}\,d\Omega ,}$
(2.3)

where ${\textstyle V}$ stands for the volume of the RVE:

 ${\displaystyle V=\displaystyle \int _{\Omega }\,d\Omega {.}}$
(2.4)

By virtue of Eq.(2.2), this condition is equivalent to impose the volume average of the fluctuation contribution ${\textstyle \nabla ^{s}{{\boldsymbol {u}}(\mathbf {x} )}}$ to vanish:

 ${\displaystyle \displaystyle \int _{\Omega }\nabla ^{s}{{\boldsymbol {u}}(\mathbf {x} )}\,d\Omega ={\boldsymbol {0}}.}$
(2.5)

(1) Some remarks concerning notation are in order here. Firstly, macroscopic variables will be identified by appending a subscript “M”, while variables associated to the fine scale will be designated by bare symbols. For instance, we shall write ${\displaystyle {\boldsymbol {\epsilon _{M}}}}$ and ${\displaystyle {\boldsymbol {\epsilon }}(\mathbf {x} )}$ to denote the macroscopic strain tensor and the fine-scale strain field, respectively. Secondly, readers accustomed to classical continuum mechanics notation are reminded that, in this work, the symbol ${\displaystyle {\boldsymbol {u}}}$ does not represent the displacement field, but rather the fluctuating part –-due to the presence of heterogeneities in the concerned RVE–- of such a field.

### 2.1.3 Hill-Mandel principle of macro-homogeneity

The scale bridging is completed by the Hill-Mandel principle of macro-homogeneity, that states that the stress power at any point ${\textstyle \mathbf {x} _{M}}$ of the macro-continuum must equal the volume average (over the RVE) of the corresponding fine-scale, stress power field (making, thus, both coarse- and fine-scale, continuum description energetically equivalent [8]). The variational statement of this principle is as follows [47]: let ${\textstyle {\boldsymbol {\sigma }}}$ be a statically admissible, fine-scale stress field, and ${\textstyle {\boldsymbol {\sigma _{M}}}}$ the associated macroscopic stress tensor, then, the identity

 ${\displaystyle {\boldsymbol {\sigma _{M}}}:{\boldsymbol {{\dot {\epsilon }}_{M}}}={\dfrac {1}{V}}\displaystyle \int _{\Omega }{\boldsymbol {\sigma }}(\mathbf {x} ):{\dot {\boldsymbol {\epsilon }}}(\mathbf {x} )\,d\Omega }$
(2.6)

must hold for any kinematically admissible strain rate ${\textstyle {\dot {\boldsymbol {\epsilon }}}}$. Inserting the rate form of Eq.(2.2) into the above equation, and by virtue of Eq.(2.3), it follows that, for identity (2.6) to hold, the macroscopic stress tensor ${\textstyle {\boldsymbol {\sigma _{M}}}}$ must be the volume average of the microscopic stress field ${\textstyle {\boldsymbol {\sigma }}={\boldsymbol {\sigma }}(\mathbf {x} ),\;\mathbf {x} \in \Omega }$:

 ${\displaystyle {\boldsymbol {\sigma _{M}}}={\dfrac {1}{V}}\displaystyle \int _{\Omega }{\boldsymbol {\sigma }}\,d\Omega {.}}$
(2.7)

Another necessary condition for Eq.(2.6) to hold is that:

 ${\displaystyle \displaystyle \int _{\Omega }{\boldsymbol {\sigma }}:\nabla ^{s}{\boldsymbol {\dot {u}}}\,d\Omega =0}$
(2.8)

for any kinematically admissible displacement fluctuation field ${\textstyle {\boldsymbol {\dot {u}}}}$. It is shown in Ref. [47] that this condition amounts to requiring that the external surface traction and body force field in the RVE be purely reactive –-i.e., a reaction to the kinematic constraints imposed upon the RVE. This is why the upcoming cell equilibrium equation does not contain neither external boundary traction nor body force terms.

## 2.2 RVE equilibrium problem

### 2.2.1 Boundary conditions

To address the issue of boundary conditions (BCs), it proves first convenient to, using Gauss' theorem, recast equation (2.5) in terms of boundary displacement fluctuations:

 ${\displaystyle \displaystyle \int _{\Omega }\nabla ^{s}{{\boldsymbol {u}}(\mathbf {x} )}\,d\Omega =\displaystyle \int _{\partial \Omega }{\boldsymbol {u}}(\mathbf {x} )\otimes ^{s}{\boldsymbol {n}}(\mathbf {x} )\,d\Omega =0,}$
(2.9)

where ${\textstyle \partial \Omega }$ represents the boundary of ${\textstyle \Omega }$, ${\textstyle {\boldsymbol {n}}}$ is the outer unit normal vector to ${\textstyle \partial \Omega }$ and the symbol ${\textstyle \otimes ^{s}}$ denotes the symmetrized tensor product. Any type of boundary conditions prescribed on the cell must obey this condition.

The natural choice for a repeating unit cell (RUC) –-defined as the smallest volume that allows one to generate, by periodic repetition, the entire microstructure of a periodic media [49]–- is to employ periodic boundary conditions for the displacement fluctuations. By definition, the boundary of an RUC comprises ${\textstyle K}$ pairs of identical –-appropriately shifted in space–- surfaces ${\textstyle \partial \Omega _{k}^{+}}$ and ${\textstyle \partial \Omega _{k}^{-}}$ (${\textstyle k=1,2\ldots K}$), with the property that for every ${\textstyle \mathbf {x} _{k}^{+}\in \partial \Omega _{k}^{+}}$, ${\textstyle {\boldsymbol {n}}(\mathbf {x} _{k}^{+})=-{\boldsymbol {n}}(\mathbf {x} _{k}^{-})}$, ${\textstyle \mathbf {x} _{k}^{-}}$ being, loosely speaking, the “counterpart” of ${\textstyle \mathbf {x} _{k}^{+}}$ in ${\textstyle \partial \Omega _{k}^{-}}$. Periodicity of displacement fluctuations implies that for every ${\textstyle \mathbf {x} _{k}^{+}\in \partial \Omega _{k}^{+}}$, ${\textstyle {\boldsymbol {u}}(\mathbf {x} _{k}^{+})={\boldsymbol {u}}(\mathbf {x} _{k}^{-})}$. (See Refs. [8,49] for more details on this type of BCs).

In statistically homogeneous micro-structures, by contrast, the situation is not so clear-cut. The concept of RVE admits various, alternative interpretations and, as a consequence, there is a certain latitude in the choice of boundary conditions (vanishing fluctuations, uniform tractions, quasi-periodic conditions …); the reader interested in this topic is urged to consult Refs. [45,50,51]. Arguably, the most practical choice –-from an implementational point of view–- in a strain-driven finite element context (and the one adopted here) is to use vanishing boundary conditions for the displacement fluctuations (${\textstyle {\boldsymbol {u}}(\mathbf {x} )=0,\forall \mathbf {x} \in \partial \Omega }$), and correspondingly determine the RVE as the smallest subvolume of the statistically homogeneous microstructure whose mechanical response is, under this type of boundary conditions, indistinguishable from that of the material-at-large.

### 2.2.2 Trial and test function spaces

It is not hard to see that both periodic and vanishing BCs are just particular instances of the more general case of homogeneous boundary conditions [52], i.e, conditions of the form:

 ${\displaystyle {\boldsymbol {A}}_{0}{\boldsymbol {u}}={\boldsymbol {0}},\;\;on\;\;\;\partial \Omega ,}$
(2.10)

${\textstyle {\boldsymbol {A}}_{0}}$ being the corresponding linear operator. An immediate corollary of this observation is that the set of all kinematically admissible displacement fluctuation fields, henceforth denoted by ${\textstyle {\mathcal {V}}_{u}}$, will form for both type of boundary conditions a vector space; this space is defined formally as:

 ${\displaystyle {\mathcal {V}}_{u}=\left\lbrace {\boldsymbol {u}}\in H^{1}(\Omega )^{d}\;|\;{\boldsymbol {A}}_{0}{\boldsymbol {u}}={\boldsymbol {0}},\;\;on\;\;\;\partial \Omega \right\rbrace .}$
(2.11)

Here, ${\textstyle H^{1}(\Omega )^{d}}$ stands for the Sobolev space of functions possessing square integrable derivatives over ${\textstyle \Omega }$. Since the test functions ${\textstyle {\boldsymbol {\eta }}}$ are kinematically admissible variations, we can write

 ${\displaystyle {\boldsymbol {\eta }}={u}-{v},\;\;\;{u},{v}\in {\mathcal {V}}_{u},}$
(2.12)

from which it follows that ${\textstyle {\boldsymbol {\eta }}\in {\mathcal {V}}_{u}}$, i.e., the spaces of trial and test functions coincide.

Observation 1: As will be further explained later, ${\textstyle {\mathcal {V}}_{u}}$ having structure of linear space confers a unique advantage –-not enjoyed by ROMs of BVPs with general inhomogeneous boundary conditions [25]–- for model reduction purposes: reduced-order responses will invariably and automatically conform to the imposed boundary conditions, regardless of the level of approximation.

### 2.2.3 Variational statement

Consider a time discretization of the interval of interest ${\textstyle [t_{0},t_{f}]=\bigcup _{n=1}^{n_{stp}}[t_{n},t_{n+1}]}$. The current value of the microscopic stress tensor ${\textstyle {\boldsymbol {\sigma }}_{n+1}}$ at each ${\textstyle \mathbf {x} \in \Omega }$ is presumed to be entirely determined by, on the one hand, the current value of the microscopic strain tensor ${\textstyle {\boldsymbol {\epsilon }}_{n+1}(\mathbf {x} )={\boldsymbol {\epsilon _{M}}}_{n+1}+\nabla ^{s}{{\boldsymbol {u}}_{n+1}(\mathbf {x} )}}$, and, on the other hand, a set of microscopic internal variables ${\textstyle {\boldsymbol {\xi }}_{n+1}}$ –-that encapsulate the history of microscopic deformations. The relationship between these variables is established by (phenomenological) rate constitutive equations; these equations may vary from point to point within the cell (multiphase materials). Likewise, the considered RVE may contain also voids distributed all over the domain. The (incremental) RVE equilibrium problem at time ${\textstyle t_{n+1}}$ can be stated as follows (see Ref. [47]): given the initial data ${\textstyle \{{\boldsymbol {u}}_{n}(\mathbf {x} ),{\boldsymbol {\epsilon _{M}}}_{n},{\boldsymbol {\xi }}_{n}\!(\mathbf {x} )\}}$ and the prescribed macroscopic strain tensor ${\textstyle {\boldsymbol {\epsilon _{M}}}_{n+1}}$, find ${\textstyle {\boldsymbol {u}}_{n+1}\in {\mathcal {V}}_{u}}$ such that

 ${\displaystyle \displaystyle \int _{\Omega }\nabla ^{s}{\boldsymbol {\eta }}:{\boldsymbol {\sigma }}_{n+1}({\boldsymbol {\epsilon _{M}}}_{n+1}+\nabla ^{s}{{\boldsymbol {u}}_{n+1}},{\boldsymbol {\xi }}_{n+1})\,d\Omega =0,}$
(2.13)

for all ${\textstyle {\boldsymbol {\eta }}\in {\mathcal {V}}_{u}}$. The actual output of interest in this fine-scale BVP is not the displacement fluctuation field per se, but rather the macroscopic stress tensor ${\textstyle \left.{\boldsymbol {\sigma _{M}}}\right|_{n+1}}$:

 ${\displaystyle \left.{\boldsymbol {\sigma _{M}}}\right|_{n+1}={\dfrac {1}{V}}\displaystyle \int _{\Omega }{\boldsymbol {\sigma }}_{n+1}\,d\Omega {.}}$
(2.14)

In order to keep the notation uncluttered, the superindex “n+1” will be hereafter dropped out and all quantities will be assumed to be evaluated at time ${\textstyle t_{n+1}}$; only when confusion is apt to show up, the pertinent distinction will be introduced.

## 2.3 Finite element formulation

Let ${\textstyle \Omega =\bigcup _{n=1}^{n_{e}}\Omega ^{e}}$ be a finite element discretization of the cell. It will be assumed that this discretization is fine enough to consider the exact and FE approximated solutions indistinguishable at the accuracy level of interest. Let ${\textstyle \left\lbrace N_{1}(\mathbf {x} ),N_{2}(\mathbf {x} )\ldots N_{N}(\mathbf {x} )\right\rbrace }$ (${\textstyle n}$ denotes the number of nodes of the discretization) be a set of shape functions associated to this discretization such that

 ${\displaystyle {\mathcal {V}}_{u}^{h}={\textrm {span}}\!\left\lbrace N_{I}\right\rbrace _{I=1}^{n}\subset {\mathcal {V}}_{u}.}$
(2.15)

In the case of vanishing displacement fluctuations, this can be achieved by simply ensuring that:

 ${\displaystyle N_{I}(\mathbf {x} )=0,\;\;\;\forall \mathbf {x} _{I}\in \partial \Omega .\;\;}$
(2.16)

As for periodic boundary conditions, let us suppose, for simplicity, that the spatial grid is such that every node ${\textstyle \mathbf {x} ^{+}\in \partial \Omega _{k}^{+}}$ (${\textstyle k=1,2\ldots K}$) has a “counterpart” node ${\textstyle \mathbf {x} ^{-}\in \partial \Omega _{k}^{-}}$ (and vice versa), and that, in addition, no nodes are placed at the intersection of adjoining surfaces. Such being the case, it can be readily shown that the conditions that the shape functions of a given node ${\textstyle \mathbf {x} _{I}\in \partial \Omega _{k}^{+}}$ and its counterpart ${\textstyle \mathbf {x} _{L}\in \partial \Omega _{k}^{-}}$ have to satisfy for proviso (2.15) to hold are:

 ${\displaystyle N_{I}(\mathbf {x} _{J})=\delta _{IJ}+\delta _{LJ}\;\;\;\;(J=1,2\ldots n)}$
(2.17)

 ${\displaystyle N_{L}(\mathbf {x} _{J})=0\;\;\;(J=1,2\ldots n).}$
(2.18)

The reader is referred to Refs. [49,53] for guidelines on how to enforce periodicity conditions in a more general scenario.

Now we approximate ${\textstyle {\boldsymbol {u}}\in {\mathcal {V}}_{u}}$ and ${\textstyle {\boldsymbol {\eta }}\in {\mathcal {V}}_{u}}$ as:

 ${\displaystyle {\boldsymbol {u}}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}})\approx {\boldsymbol {u}}^{(h)}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}})=\displaystyle \sum _{I=1}^{n}\!N_{I}(\mathbf {x} ){\boldsymbol {U}}_{I}({\boldsymbol {\epsilon _{M}}}),}$
(2.19)

 ${\displaystyle {\boldsymbol {\eta }}(\mathbf {x} )\approx {\boldsymbol {\eta ^{(h)}}}(\mathbf {x} )=\displaystyle \sum _{I=1}^{n}\!N_{I}(\mathbf {x} ){\boldsymbol {\eta }}_{I},}$
(2.20)

where ${\textstyle {\boldsymbol {U}}_{I}\in \mathbb {R} ^{d}}$ and ${\textstyle {\boldsymbol {\eta }}_{I}\in \mathbb {R} ^{d}}$ (${\textstyle I=1,2\ldots n}$) denote the nodal values of the displacement fluctuations and test functions, respectively. Inserting these approximations in Eq.(2.13), and exploiting the arbitrariness of coefficients ${\textstyle {\boldsymbol {\eta }}_{I}}$ (${\textstyle I=1,2\ldots n}$), one arrives at the following set of discrete equilibrium equations (repeated indices implies summation):

 ${\displaystyle \displaystyle \int _{\Omega }{\dfrac {\partial N_{I}}{\partial x_{j}}}{\boldsymbol {\sigma }}_{ij}({\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{(h)}},{\boldsymbol {\xi }})\,d\Omega ={0}\;\;(i=1\ldots d;\;I=1\ldots n).}$
(2.21)

Introducing Voigt's notation1, the above equation can be expressed in matrix format as:

 ${\displaystyle \displaystyle \int _{\Omega }{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}({\boldsymbol {\epsilon _{M}}}+{\boldsymbol {B}}{\boldsymbol {U}},{\boldsymbol {\xi }})\,d\Omega ={\boldsymbol {0}}.}$
(2.22)

where ${\textstyle {\boldsymbol {\sigma }}:\Omega \rightarrow \mathbb {R} ^{s}}$ represents, with a slight abuse of notation2, the column matrix form of the stress tensor ( ${\textstyle s=4}$ and ${\textstyle s=6}$ for plane and 3D problems, respectively), and ${\textstyle {\boldsymbol {B}}:\Omega \rightarrow \mathbb {R} ^{s\times n\cdot d}}$ is the classical, global ${\textstyle {\boldsymbol {B}}}$-matrix" connecting strains at a given point with the vector ${\textstyle {\boldsymbol {U}}\in \mathbb {R} ^{n\cdot d}}$ containing all nodal displacements:

 ${\displaystyle {\boldsymbol {U}}{\mathrel {\mathop {:}}}={\begin{bmatrix}{\boldsymbol {U}}_{1}\\{\boldsymbol {U}}_{2}\\\vdots \\{\boldsymbol {U}}_{n}\end{bmatrix}}.}$
(2.23)

As usual, numerical evaluation of the integral in Eq.(2.22) is carried out by Gaussian quadrature:

 ${\displaystyle \displaystyle \int _{\Omega }{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}\,d\Omega \approx \displaystyle \sum _{g=1}^{n_{g}}w_{g}{{\boldsymbol {B}}^{T}}\!\!(\mathbf {x} _{g}){\boldsymbol {\sigma }}(\mathbf {x} _{g},;)={\boldsymbol {0}}.}$
(2.24)

Here, ${\textstyle n_{g}={\mathcal {O}}(n)}$ stands for the total number of Gauss points of the mesh; ${\textstyle w_{g}}$ denotes the weight associated to the ${\textstyle g-th}$ Gauss point ${\textstyle \mathbf {x} _{g}}$ (this weight includes both the quadrature weight itself and the corresponding Jacobian determinant.); and ${\textstyle {\boldsymbol {B}}(\mathbf {x} _{g})}$ and ${\textstyle {\boldsymbol {\sigma }}(\mathbf {x} _{g},;)}$ stand for the B-matrix and the stress vector at Gauss point ${\textstyle \mathbf {x} _{g}}$, respectively.

(1) Here, it is convenient to use the so-called modified Voigt's notation rather than the standard one. In the modified Voigt's notation, both stress ${\displaystyle {\boldsymbol {\sigma }}}$ and strain ${\displaystyle {\boldsymbol {\epsilon }}}$ tensors are represented as column vectors (${\displaystyle \left\lbrace {\boldsymbol {\sigma }}\right\rbrace }$ and ${\displaystyle \left\lbrace {\boldsymbol {\epsilon }}\right\rbrace }$, respectively ) in which the shear components are multiplied by ${\displaystyle {\sqrt {2}}}$. The advantage of this notation over the conventional, engineering Voigt's notation is the equivalence between norms; viz., ${\displaystyle \Vert {\boldsymbol {\sigma }}\Vert ={\sqrt {{\boldsymbol {\sigma }}:{\boldsymbol {\sigma }}}}=\Vert \left\lbrace {\boldsymbol {\sigma }}\right\rbrace \Vert ={\sqrt {\left\lbrace {\boldsymbol {\sigma }}\right\rbrace ^{T}\left\lbrace {\boldsymbol {\sigma }}\right\rbrace }}}$. The reader is urged to consult [54] for further details on this notation.

(2) The same symbol is used to denote the stress tensor and its counterpart in Voigt's notation.

# 3 Reduced-order model of the RVE

## 3.1 Computation of reduced basis

A basic, intuitive picture of the strategy for computing the reduced basis onto which to project the cell equilibrium equation (2.13) was already given in the introductory Chapter. In the following, we put the idea behind this strategy on a more rigorous footing. We begin by noting that, from a functional analysis standpoint, the term model reduction is conceptually akin to the more common term model discretization, since both connote transitions from higher-dimensional to lower-dimensional solution spaces. Indeed, whereas model discretization is used to refer to the (classical) passage from the infinite dimensional space ${\textstyle {\mathcal {V}}_{u}}$ to the finite element subspace ${\textstyle {\mathcal {V}}_{u}^{h}\subset {\mathcal {V}}_{u}}$, model reduction denotes a transition from this finite dimensional space ${\textstyle {\mathcal {V}}_{u}^{h}}$ to a significantly smaller manifold ${\textstyle {\mathcal {V}}_{u}^{*}\subset {\mathcal {V}}_{u}^{h}}$ –-the reduced-order space. This latter transition is not carried out directly, but in two sequential steps, namely, sampling of the parameter space and dimensionality reduction. The precise meaning of these terms is explained below.

### 3.1.1 Sampling of the input parameter space

In constructing the finite element space of kinematically admissible functions ${\textstyle {\mathcal {V}}_{u}^{h}}$, the only restrictions placed on the motion of the mesh nodes are those imposed at the boundaries through conditions (2.16) (for RVEs) and (2.17), (2.18) (for RUCs). The finite element solution space, thus, does not presuppose any constraint on the motion of the interior nodes of the mesh.

However, in actuality, interior nodes cannot fluctuate freely, independently from each other, but they rather move according to deformational patterns dictated by the constitutive laws that govern the mechanical behavior of the distinct phases in the cell –-as noted by Lubliner [55], constitutive laws can be regarded as internal restrictions on the kinds of deformation a body can suffer. This means that the solution of the finite element equilibrium equation (2.13) for given values of the macro-strain tensor ${\textstyle {\boldsymbol {\epsilon _{M}}}}$ actually lives in a smaller subspace ${\textstyle {\mathcal {V}}_{u}^{\epsilon }\subset {\mathcal {V}}_{u}^{h}}$ (in the parlance of model reduction [23,56], ${\textstyle {\mathcal {V}}_{u}^{\epsilon }}$ is the manifold induced by the parametric dependence of the BVP on the input variables).

Yet, in general, this manifold cannot be precisely determined –-such a task would require finite element analyses of the cell under all conceivable strain paths. Rather, one has to be content to construct an approximation of it as the span of the displacement fluctuation solutions obtained for a judiciously chosen set of ${\textstyle n_{hst}}$ input strain histories ${\textstyle \{{^{\,t}}{\boldsymbol {\epsilon _{M}}}{_{}}^{1},{^{\,t}}{\boldsymbol {\epsilon _{M}}}{_{}}^{2},\ldots {^{\,t}}{\boldsymbol {\epsilon _{M}}}{_{}}^{n_{hst}}\}}$. Suppose, for simplicity, that each of these strain histories is discretized into equal number of steps ${\textstyle n_{stp}}$, and let

 ${\displaystyle {\boldsymbol {u}}^{k}(\mathbf {x} )={\boldsymbol {u}}(\mathbf {x} ;{^{\,t}}{\boldsymbol {\epsilon _{M}}}{_{j}}^{i}),\;\;k=(i-1)n_{hst}+j}$
(3.1)

denote the displacement fluctuation solution at the ${\textstyle j-th}$ time step of the ${\textstyle i-th}$ strain history (${\textstyle i=1,2\ldots n_{hst}}$, ${\textstyle j=1,2\ldots n_{stp}}$). The approximating space for ${\textstyle {\mathcal {V}}_{u}^{\epsilon }}$, henceforth called the snapshots space, is then defined as:

 ${\displaystyle {\mathcal {V}}_{u}^{snp}={\textrm {span}}\left\lbrace {\boldsymbol {u}}^{1}(\mathbf {x} ),{\boldsymbol {u}}^{2}(\mathbf {x} ),\ldots {\boldsymbol {u}}^{n_{snp}}(\mathbf {x} )\right\rbrace \subseteq {{\mathcal {V}}_{u}^{\epsilon }},}$
(3.2)

${\textstyle n_{snp}=n_{stp}n_{hst}}$ being the total number of snapshots. Likewise, the matrix containing, in columns, the nodal values of these displacement fluctuations solutions:

 ${\displaystyle {\boldsymbol {X}}_{u}={\begin{bmatrix}{\boldsymbol {U}}^{1}&{\boldsymbol {U}}^{2}&\cdots &{\boldsymbol {U}}^{n_{snp}}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {U}}_{1}^{1}&{\boldsymbol {U}}_{1}^{2}&\cdots &{\boldsymbol {U}}_{1}^{n_{snp}}\\{\boldsymbol {U}}_{2}^{1}&{\boldsymbol {U}}_{2}^{2}&\cdots &{\boldsymbol {U}}_{2}^{n_{snp}}\\\vdots &\vdots &\vdots &\vdots \\{\boldsymbol {U}}_{n}^{1}&{\boldsymbol {U}}_{n}^{2}&\cdots &{\boldsymbol {U}}_{n}^{n_{snp}}\\\end{bmatrix}}\in \mathbb {R} ^{n\cdot d\times n_{snp}}}$
(3.3)

will correspondingly be termed the (displacement fluctuations) snapshot matrix.

Remark 1: The first error incurred in solving the cell equilibrium equations using the reduced basis approach arises in approximating ${\textstyle {\mathcal {V}}_{u}^{\epsilon }}$ by this space of snapshots ${\textstyle {\mathcal {V}}_{u}^{snp}}$. In order to keep this error to a minimum, one should strive to select the set of strain histories in such a way that the span of the corresponding displacement fluctuation solutions cover as much as possible the space ${\textstyle {\mathcal {V}}_{u}^{\epsilon }}$ (or at least, the region or regions of particular interest), while, at the same time, trying to keep the total number of snapshots in check –-the computational cost of the subsequent dimensionality reduction process grows considerably with the size of the snapshot matrix. In this respect, it may be interesting to note that this task of sampling the input parameter space (also known as “training”1) is somehow akin to the experimental process whereby material parameters of standard phenomenological models are calibrated in a laboratory. In this analogy, the RVE plays the role of the corresponding experimental specimen, whereas the macro-strain training trajectories represent the loading paths of the pertinent calibration tests. As opposed to the situation encountered in standard laboratory experiments, however, in the training process, one has “privileged” information regarding the phenomenological behavior of the constituents. Hindsight and elementary physical considerations can therefore aid in restricting the number of strain histories (and hence of snapshots) necessary to characterize the response. For instance, if the behavior of the materials that compose the cell is governed by rate-independent constitutive models, we know beforehand that it is not necessary to study the response under varying rates of deformation. Strategies for efficiently sampling the input parameter space in general model reduction contexts can be found in Refs. [58,59,60,61].

(1) The term “training”, which, incidentally, is borrowed from the neural network literature [57], is used throughout the text to refer to the offline generation of snapshots.

### 3.1.2 Dimensionality reduction

The next and definitive step in the transition from the high-dimensional finite element space ${\textstyle {\mathcal {V}}_{u}^{h}}$ to the desired reduced-order space ${\textstyle {\mathcal {V}}_{u}^{*}}$ –-in which the cell BVP is to be finally posed–- is the dimensionality reduction process, in which, as pointed out in the introductory Chapter, the dominant deformational patterns of the cell response are identified and unveiled by washing out the “inessentials”.

#### 3.1.2.1 Proper Orthogonal Decomposition

To accomplish this central task, we employ here the Proper Orthogonal Decomposition (POD). The formal statement of the POD problem goes as follows: given the ensemble of snapshots ${\textstyle \{{\boldsymbol {u}}^{1},{\boldsymbol {u}}^{2},\ldots {\boldsymbol {u}}^{n_{snp}}\}}$, find a set of ${\textstyle n_{u} orthogonal basis functions ${\textstyle \{{\boldsymbol {\mathit {\Phi }}}_{1},{\boldsymbol {\mathit {\Phi }}}_{2},\ldots {\boldsymbol {\mathit {\Phi }}}_{n_{u}}\}}$ (${\textstyle {\boldsymbol {\mathit {\Phi }}}_{i}\in {\mathcal {V}}_{u}^{snp}}$) such that the error defined as

 ${\displaystyle {e}_{u}(n_{u}){\mathrel {\mathop {:}}}={\sqrt {\displaystyle \sum _{k=1}^{n_{snp}}\Vert {\boldsymbol {u}}^{k}-\mathbf {P} ^{*}{\boldsymbol {u}}^{k}\Vert _{L_{2}(\Omega )}^{2}}}}$
(3.4)

is minimized. Here, ${\textstyle \mathbf {P} ^{*}{\boldsymbol {u}}^{k}}$ represents the projection of ${\textstyle {\boldsymbol {u}}^{k}}$ onto the subspace spanned by the basis functions ${\textstyle \{{\boldsymbol {\mathit {\Phi }}}_{i}\}_{i=1}^{n_{u}}}$, and ${\textstyle \Vert \cdot \Vert _{L_{2}(\Omega )}}$ symbolizes the ${\textstyle L_{2}}$ norm. It is shown in Ref. [62] that the solution of this optimization problem can be obtained by first solving the following eigenvalue problem:

 ${\displaystyle {\boldsymbol {L}}{v}=\lambda {v},\;\;\;\Vert {v}\Vert =1,}$
(3.5)

where ${\textstyle {\boldsymbol {L}}\in \mathbb {R} ^{n_{snp}\times n_{snp}}}$ is a symmetric matrix defined as:

 ${\displaystyle {\boldsymbol {L}}_{ij}{\mathrel {\mathop {:}}}=\left\langle {{\boldsymbol {u}}^{i}},{{\boldsymbol {u}}^{j}}\right\rangle _{L_{2}(\Omega )}=\displaystyle \int _{\Omega }{\boldsymbol {u}}^{i}{\boldsymbol {u}}^{j}\,d\Omega ,}$
(3.6)

i.e., ${\textstyle {\boldsymbol {L}}_{ij}}$ is the ${\textstyle L_{2}}$ inner product between snapshots ${\textstyle {\boldsymbol {u}}^{i}}$ and ${\textstyle {\boldsymbol {u}}^{j}}$. In statistical terms, ${\textstyle {\boldsymbol {L}}}$ can be interpreted as a covariance matrix: the off-diagonal entries captures the degree of linear correlation or redundancy between pair of snapshots (the covariance), whereas the diagonal terms are the variances [63]. The goal in diagonalizing ${\textstyle {\boldsymbol {L}}_{ij}}$ by solving the eigenvalue problem (3.5) is to re-express the snapshots data in axes that filter out the redundancies and reveals the actual dominant displacement fluctuation patterns.

Since ${\textstyle {\boldsymbol {u}}^{i}(\mathbf {x} )=\sum _{I=1}^{n}N_{I}(\mathbf {x} ){\boldsymbol {U}}_{I}^{i}}$, expression (3.6) for the covariance matrix can be rewritten in terms of the snapshot matrix as follows:

 ${\displaystyle {\boldsymbol {L}}_{ij}=\displaystyle \sum _{I=1}^{n}\sum _{J=1}^{n}{\boldsymbol {U}}_{I}^{i}{\boldsymbol {U}}_{J}^{j}\int _{\Omega }N_{I}(\mathbf {x} )N_{J}(\mathbf {x} )\,d\Omega ,}$
(3.7)

or in matrix form:

 ${\displaystyle {\boldsymbol {L}}={\boldsymbol {X}}_{u}^{T}\mathbf {M} {\boldsymbol {X}}_{u},}$
(3.8)

where

 ${\displaystyle \mathbf {M} _{IJ}{\mathrel {\mathop {:}}}=\displaystyle \int _{\Omega }N_{I}(\mathbf {x} )N_{J}(\mathbf {x} )\,d\Omega \;\;\;\;\;\;\;I,J=1,2\ldots n.}$
(3.9)

Note that, except for the density factor, this matrix ${\textstyle \mathbf {M} }$ is similar to the “mass matrix” appearing in finite element implementations of dynamical problems.

Once the eigenvalues problem (3.5) has been solved, the desired reduced basis ${\textstyle \{{\boldsymbol {\mathit {\Phi }}}_{1},{\boldsymbol {\mathit {\Phi }}}_{2},\ldots {\boldsymbol {\mathit {\Phi }}}_{n_{u}}\}}$ is calculated from the ${\textstyle n_{u}}$ largest eigenvalues ${\textstyle \lambda _{1}\geq \lambda _{2}\geq \ldots \lambda _{n_{u}-1}\geq \lambda _{n_{u}}>0}$ and associated eigenvectors through the following expression:

 ${\displaystyle {\boldsymbol {\mathit {\Phi }}}_{i}(\mathbf {x} )=\displaystyle \sum _{k=1}^{n_{snp}}{\boldsymbol {u}}^{k}(\mathbf {x} ){\dfrac {{v}_{i}}{\sqrt {\lambda _{i}}}},\;\;\;i=1,2\ldots n_{u},}$
(3.10)

(modes are judged to be essential or dominant, and hence worthy of being included in the reduced basis set, if their associated eigenvalues have relatively large magnitudes). Substitution of ${\textstyle {\boldsymbol {u}}^{k}(\mathbf {x} )=\sum _{I=1}^{n}N_{I}(\mathbf {x} ){\boldsymbol {U}}_{I}^{k}}$ into the above equation yields:

 ${\displaystyle {\boldsymbol {\mathit {\Phi }}}_{i}(\mathbf {x} )=\displaystyle \sum _{I=1}^{n}N_{I}(\mathbf {x} )\sum _{k=1}^{n_{snp}}{\boldsymbol {U}}_{I}^{k}{\dfrac {{v}_{i}}{\sqrt {\lambda _{i}}}}=\sum _{I=1}^{n}N_{I}(\mathbf {x} ){\boldsymbol {\Phi }}_{Ii},\;\;\;i=1,2\ldots n_{u},}$
(3.11)

where1 ${\textstyle {\boldsymbol {\Phi }}_{Ii}\in \mathbb {R} ^{d}}$ is given by

 ${\displaystyle {\boldsymbol {\Phi }}_{Ii}=\sum _{k=1}^{n_{snp}}{\boldsymbol {U}}_{I}^{k}{\dfrac {{v}_{i}}{\sqrt {\lambda _{i}}}}\;\;\;I=1,2\ldots n,\;i=1,2\ldots n_{u},}$
(3.12)

and stands for the value of the basis function ${\textstyle {\boldsymbol {\mathit {\Phi }}}_{i}(\mathbf {x} )}$ at the ${\textstyle I-th}$ node of the finite element grid. The matrix ${\textstyle {\boldsymbol {\Phi }}\in \mathbb {R} ^{n\cdot d\times n_{u}}}$ defined by the above equation will be hereafter called the reduced basis matrix. Each column ${\textstyle {\boldsymbol {\Phi }}_{i}\in \mathbb {R} ^{n\cdot d}}$ (${\textstyle i=1,2\ldots n_{u}}$) of this matrix can be compactly expressed in terms of the snapshot matrix as follows:

 ${\displaystyle {\boldsymbol {\Phi }}_{i}={\boldsymbol {X}}_{u}{\dfrac {{v}_{i}}{\sqrt {\lambda _{i}}}},\;\;\;\;i=1,2\ldots n_{u}.}$
(3.13)

#### 3.1.2.2 Singular value decomposition

Instead of solving the eigenvalue problem (3.5), and then obtaining the reduced basis matrix ${\textstyle {\boldsymbol {\Phi }}}$ from expression (3.12), one can alternatively compute this basis matrix using the Singular Value Decomposition. Indeed, let ${\textstyle \mathbf {M} ={\bar {\mathbf {M} }}^{T}{\bar {\mathbf {M} }}}$ be the Cholesky decomposition of ${\textstyle \mathbf {M} }$, and let ${\textstyle {\boldsymbol {\bar {X}}}_{u}}$ denote the matrix defined as:

 ${\displaystyle {\boldsymbol {\bar {X}}}_{u}{\mathrel {\mathop {:}}}={\bar {\mathbf {M} }}{\boldsymbol {X}}_{u}.}$
(3.14)

It can be shown (see Appendix A) that the ${\textstyle i-th}$ column of the reduced basis matrix ${\textstyle {\boldsymbol {\Phi }}}$ is related to the ${\textstyle i-th}$ left singular vector of ${\textstyle {\boldsymbol {\bar {X}}}_{u}}$, denoted by ${\textstyle {\boldsymbol {\bar {U}}}_{i}}$, through expression

 ${\displaystyle {\boldsymbol {\Phi }}_{i}={\bar {\mathbf {M} }}^{-1}{\boldsymbol {\bar {U}}}_{i},\;\;\;\;i=1,2\ldots n_{u}.}$
(3.15)

#### 3.1.2.3 Elastic/Inelastic reduced basis functions

The POD can be viewed as a multidimensional data fitting procedure intended to obtain a sequence of orthogonal basis functions whose span best approximate the space of snapshots. As such, the POD is a purely data-driven process, “agnostic” to the physical origin of the data [64,63]. For instance, for POD basis construction purposes, it is completely immaterial whether a given snapshot corresponds to a purely linear elastic solution or to a solution well into the inelastic regime. The task of discriminating which features of the cell response are essential and which are not is exclusively guided by statistical considerations: if the elastic response happens to be poorly represented within the snapshot ensemble, the POD may regard as unimportant the contribution of these snapshots, and, as a consequence, the basis functions with largest associated eigenvalues –-i.e., the essential modes–- would hardly contain any information of this range. To accurately replicate the apparently trivial linear elastic behavior, thus, one may be forced to take a relatively large number of basis functions, and this may translate into a significant increase in the overall online computational cost. This fact certainly places the POD-based reduced basis approach at a competitive disadvantage compared with semi-analytical homogenization approaches such as the Nonlinear Transformation Field Analysis [20], which do capture exactly (and effortlessly) the linear elastic response of the cell.

To eliminate this shortcoming, we propose here a slightly different strategy for constructing the reduced basis. The essence of the proposal is to partition the space of snapshots ${\textstyle {\mathcal {V}}_{u}^{snp}}$ into elastic (${\textstyle {\mathcal {V}}_{u,el}^{snp}}$) and inelastic (${\textstyle {\mathcal {V}}_{u,inel}^{snp}}$) subspaces:

 ${\displaystyle {\mathcal {V}}_{u}^{snp}={\mathcal {V}}_{u,el}^{snp}\oplus {\mathcal {V}}_{u,inel}^{snp},}$
(3.16)

(${\textstyle \oplus }$ symbolizes direct sum of subspaces [65]) and then obtain the reduced basis as the union of the bases for both subspaces. Below, we describe this strategy more in detail.

The first step is to determine an orthogonal basis for ${\textstyle {\mathcal {V}}_{u,el}^{snp}}$. One can do this by simply performing ${\textstyle m_{e}}$ independent, linear elastic finite element analysis of the cell (${\textstyle m_{e}=6}$ for 3D problems, and ${\textstyle m_{e}=3}$ for plane strain), and then orthonormalizing the resulting displacement fluctuation fields. These ${\textstyle m_{e}}$ elastic modes will be considered as the first ${\textstyle m_{e}}$ basis functions of the reduced basis:

 ${\displaystyle {\textrm {span}}\!\left\lbrace {\boldsymbol {\mathit {\Phi }}}_{1},{\boldsymbol {\mathit {\Phi }}}_{2},\ldots ,{\boldsymbol {\mathit {\Phi }}}_{m_{e}}\right\rbrace ={\mathcal {V}}_{u,el}^{snp}.}$
(3.17)

Once we have at our disposal this set of elastic basis functions, we compute the (orthogonal) projection of each snapshot ${\textstyle {\boldsymbol {u}}^{k}}$ onto the orthogonal complement of ${\textstyle {\mathcal {V}}_{u,el}^{snp}}$ (which is precisely the aforementioned inelastic space ${\textstyle {\mathcal {V}}_{u,inel}^{snp}}$):

 ${\displaystyle {\boldsymbol {u}}_{inel}^{k}{\mathrel {\mathop {:}}}={\boldsymbol {u}}^{k}-\displaystyle \sum _{i=1}^{m_{e}}\left\langle {{\boldsymbol {\mathit {\Phi }}}_{i}},{{\boldsymbol {u}}^{k}}\right\rangle _{L_{2}(\Omega )}{\boldsymbol {\mathit {\Phi }}}_{i},\;\;\;\;k=1,2\ldots n_{snp}.}$
(3.18)

It is now on this ensemble of inelastic snapshots ${\textstyle \{{\boldsymbol {u}}_{inel}^{k}\}_{k=1}^{n_{snp}}}$ that the previously described POD is applied to obtain the remaining ${\textstyle n_{u}-m_{e}}$ basis functions. Thus, we finally have:

 ${\displaystyle {\mathcal {V}}_{u}^{*}={\mathcal {V}}_{u,el}^{snp}\oplus {\mathcal {V}}_{u,inel}^{snp}={\textrm {span}}\{\overbrace {{\boldsymbol {\mathit {\Phi }}}_{1},{\boldsymbol {\mathit {\Phi }}}_{2},\ldots ,{\boldsymbol {\mathit {\Phi }}}_{6}} ^{\textrm {Elasticmodes}},\overbrace {{\boldsymbol {\mathit {\Phi }}}_{7},\ldots ,{\boldsymbol {\mathit {\Phi }}}_{n_{u}}} ^{\textrm {'Essential'Inelasticmodes}}\}.}$
(3.19)

for 3D problems, and

 ${\displaystyle {\mathcal {V}}_{u}^{*}={\textrm {span}}\{\overbrace {{\boldsymbol {\mathit {\Phi }}}_{1},{\boldsymbol {\mathit {\Phi }}}_{2},{\boldsymbol {\mathit {\Phi }}}_{3}} ^{\textrm {Elasticmodes}},\overbrace {{\boldsymbol {\mathit {\Phi }}}_{4},\ldots ,{\boldsymbol {\mathit {\Phi }}}_{n_{u}}} ^{\textrm {'Essential'inelasticmodes}}\}.}$
(3.20)

for plane strain.

Observation 2: In placing the ${\textstyle m_{e}}$ elastic modes within the first ${\textstyle m_{e}}$ positions, the reduced-order model is guaranteed to deliver linear elastic solutions with the same accuracy as the underlying (full-order) finite element model (obviously, provided that ${\textstyle n_{u}\geq m_{e}}$).

Further details concerning the numerical implementation of this apparently novel –-to the best of the authors' knowledge–- basis construction strategy can be found in Appendix B.

(1) It is necessary at this point to further clarify the notation employed for referring to the basis functions. A bold italic “phi” symbol with one subscript is employed to denote the basis function itself, i.e., ${\displaystyle {\boldsymbol {\mathit {\Phi }}}_{i}:\Omega \rightarrow \mathbb {R} ^{d}}$ (${\displaystyle i=1,2\ldots n_{u}}$). Bold normal “phi” symbols, on the other hand, are employed to represent values of such basis functions at the nodes of the underlying finite element mesh. For instance, the value of the ${\displaystyle i-th}$ basis functions at the ${\displaystyle I-th}$ node is symbolized as ${\displaystyle {\boldsymbol {\Phi }}_{Ii}={\boldsymbol {\mathit {\Phi }}}_{i}(\mathbf {x} _{I})\in \mathbb {R} ^{d}}$. When accompanied by only one (lowercase) subscript, the bold normal “phi” symbol denotes a column vector containing the nodal values of the pertinent basis functions at all gauss points: ${\displaystyle {\boldsymbol {\Phi }}_{i}={{\begin{bmatrix}{\boldsymbol {\Phi }}_{1i}^{T}&{\boldsymbol {\Phi }}_{2i}^{T}&\ldots &{\boldsymbol {\Phi }}_{ni}\end{bmatrix}}^{T}}\in \mathbb {R} ^{n\cdot d}}$ (${\displaystyle i=1,2\ldots n_{u}}$). Lastly, when no subscript is attached, ${\displaystyle {\boldsymbol {\Phi }}}$ represents the reduced basis matrix: ${\displaystyle {\boldsymbol {\Phi }}={\begin{bmatrix}{\boldsymbol {\Phi }}_{1}&{\boldsymbol {\Phi }}_{2}&\cdots &{\boldsymbol {\Phi }}_{n_{u}}\end{bmatrix}}\in \mathbb {R} ^{n\cdot d\times n_{u}}}$.

## 3.2 Galerkin projection onto the reduced subspace

We now seek to pose the boundary-value problem represented by Eq.(2.13) in the reduced-order space ${\textstyle {\mathcal {V}}_{u}^{*}\subseteq {\mathcal {V}}_{u}^{h}}$ spanned by the basis functions ${\textstyle \{{\boldsymbol {\mathit {\Phi }}}_{1},{\boldsymbol {\mathit {\Phi }}}_{2},\ldots ,{\boldsymbol {\mathit {\Phi }}}_{n_{u}}\}}$. To this end, we approximate both test ${\textstyle {\boldsymbol {\eta }}\in {\mathcal {V}}_{u}}$ and trial ${\textstyle {\boldsymbol {u}}\in {\mathcal {V}}_{u}}$ functions by the following linear expansions:

 ${\displaystyle {\boldsymbol {u}}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}})\approx {\boldsymbol {u}}^{*}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}})=\displaystyle \sum _{i=1}^{n_{u}}{\boldsymbol {\mathit {\Phi }}}_{i}(\mathbf {x} )U_{i}^{*}({\boldsymbol {\epsilon _{M}}}),}$
(3.21)

 ${\displaystyle {\boldsymbol {\eta }}(\mathbf {x} )\approx {\boldsymbol {\eta }}^{*}(\mathbf {x} )=\displaystyle \sum _{i=1}^{n_{u}}{\boldsymbol {\mathit {\Phi }}}_{i}(\mathbf {x} ){\eta }_{i}^{*},}$
(3.22)

${\textstyle {\boldsymbol {u}}^{*}(\mathbf {x} )}$ and ${\textstyle {\boldsymbol {\eta }}^{*}(\mathbf {x} )}$ being the low-dimensional approximations of trial and test functions, respectively (hereafter, asterisked symbols will be used to denote low-dimensional approximations of the associated variables). Inserting Eqs. (3.21) and (3.22) into Eq.(2.13), and exploiting the arbitrariness of coefficients ${\textstyle {\eta }_{i}^{*}}$ (${\textstyle i=1,2\ldots n_{u}}$), we arrive at the following set of ${\textstyle n_{u}}$ equilibrium equations:

 ${\displaystyle \displaystyle \int _{\Omega }\nabla ^{s}{{\boldsymbol {\mathit {\Phi }}}_{i}(\mathbf {x} )}:{\boldsymbol {\sigma }}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }})\,d\Omega =0,\;\;\;i=1,2\ldots n_{u}.}$
(3.23)

Expressing now the reduced basis functions in the above equation in terms of finite element shape functions (through expression ${\textstyle {\boldsymbol {\mathit {\Phi }}}_{i}(\mathbf {x} )=\sum _{I=1}^{n}N_{I}(\mathbf {x} ){\boldsymbol {\Phi }}_{Ii}}$), we get (in Voigt's notation):

 ${\displaystyle \displaystyle \int _{\Omega }{{\boldsymbol {B}}_{i}^{*}}^{T}\!\!(\mathbf {x} )\,{\boldsymbol {\sigma }}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}}+{{\boldsymbol {B}}^{*}}{\boldsymbol {U}}^{*},{\boldsymbol {\xi }})\,d\Omega ={\boldsymbol {0}},\;\;\;i=1,2\ldots n_{u},}$
(3.24)

or more compactly:

 ${\displaystyle \displaystyle \int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}\!\!(\mathbf {x} )\,{\boldsymbol {\sigma }}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}}+{{\boldsymbol {B}}^{*}}{\boldsymbol {U}}^{*},{\boldsymbol {\xi }})\,d\Omega ={\boldsymbol {0}}.}$
(3.25)

Here, ${\textstyle {\boldsymbol {U}}^{*}\in \mathbb {R} ^{n_{u}}}$ denotes the vector containing the reduced displacement fluctuations ${\textstyle U_{i}^{*}\in \mathbb {R} ^{}}$ ${\textstyle (i=1,2\ldots n_{u})}$ –-the basic unknowns of the reduced-order problem:

 ${\displaystyle {\boldsymbol {U}}^{*}{\mathrel {\mathop {:}}}={\begin{bmatrix}U_{1}^{*}\\U_{2}^{*}\\\vdots \\U_{n_{u}}^{*}\end{bmatrix}}}$
(3.26)

and ${\textstyle {{\boldsymbol {B}}^{*}}:\Omega \rightarrow \mathbb {R} ^{s\times n_{u}}}$ stands for the reduced “B-matrix”, defined as:

 ${\displaystyle {{\boldsymbol {B}}^{*}}(\mathbf {x} ){\mathrel {\mathop {:}}}={\boldsymbol {B}}(\mathbf {x} ){\boldsymbol {\Phi }},}$
(3.27)

and that connects the gradient of the displacement fluctuation field with the vector of reduced displacement fluctuations, i.e.:

 ${\displaystyle \nabla ^{s}{{\boldsymbol {u}}^{*}}=\displaystyle \sum _{i=1}^{n_{u}}{{\boldsymbol {B}}_{i}^{*}\!}U_{i}^{*}=\overbrace {\begin{bmatrix}{{\boldsymbol {B}}_{1}^{*}\!}&{{\boldsymbol {B}}_{2}^{*}\!}&\ldots &{{\boldsymbol {B}}_{n_{u}}^{*}\!}\end{bmatrix}} ^{{\boldsymbol {B}}^{*}}\overbrace {\begin{bmatrix}U_{1}^{*}\\U_{2}^{*}\\\vdots \\U_{n_{u}}^{*}\end{bmatrix}} ^{{\boldsymbol {U}}^{*}}}$ ${\displaystyle ={{\boldsymbol {B}}^{*}}{\boldsymbol {U}}^{*}={\boldsymbol {B}}{\boldsymbol {\Phi }}{\boldsymbol {U}}^{*}.}$
(3.28)

For implementational purposes, it is more expedient to express Eq.(3.27) in terms of elemental ${\textstyle B-}$matrices. To this end, we write:

 ${\displaystyle {\boldsymbol {B}}(\mathbf {x} )=\left\lbrace {\begin{array}{l}{\boldsymbol {B}}^{e}(\mathbf {x} ),\;\;\;\;{\textrm {if}}\,\mathbf {x} \in \Omega ^{e}\\{\boldsymbol {0}},\;\;\;\;\;\;\;{\textrm {otherwise}}\\\end{array}}\right.}$
(3.29)

where ${\textstyle {\boldsymbol {B}}^{e}\in \mathbb {R} ^{s\times d\cdot {\bar {n}}_{e}}}$ denotes the local ${\textstyle B}$-matrix of element ${\textstyle \Omega ^{e}}$ (${\textstyle {\bar {n}}_{e}}$, in turn, is the number of nodes in ${\textstyle \Omega ^{e}}$). Thus,

 ${\displaystyle {{\boldsymbol {B}}^{*}}(\mathbf {x} )={\boldsymbol {B}}(\mathbf {x} ){\boldsymbol {\Phi }}={\boldsymbol {B}}^{e}(\mathbf {x} ){\boldsymbol {\Phi }}^{e}.}$
(3.30)

In the above equation, ${\textstyle {\boldsymbol {\Phi }}^{e}\in \mathbb {R} ^{d{\bar {n}}_{e}\times n_{u}}}$ represents the block matrix of ${\textstyle {\boldsymbol {\Phi }}}$ corresponding to the ${\textstyle {\bar {n}}_{e}}$ nodes of finite element ${\textstyle \Omega ^{e}}$ (${\textstyle e=1,2\ldots n_{e}}$).

# 4 Numerical integration

## 4.1 Classical Gauss quadrature: the standard ROM

A straightforward –-but, as already mentioned, ostensibly inefficient–- route for numerically evaluating the integral appearing in the reduced-order equilibrium equation (3.24) is to simply use the same Gauss quadrature formulae and the same set of Gauss points as the underlying finite element model (see Eq.(2.24)):

 ${\displaystyle \displaystyle \int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}{\boldsymbol {\sigma }}\,d\Omega \approx \displaystyle \sum _{g=1}^{n_{g}}w_{g}{{{\boldsymbol {B}}^{*}}^{T}}\!\!(\mathbf {x} _{g}){\boldsymbol {\sigma }}(\mathbf {x} _{g};\cdot )={\boldsymbol {0}}.}$
(4.1)

Low-rank approximations that employ the underlying finite element Gauss points for numerically evaluating integrals in the weak statement of the problem are commonly known as standard reduced-order models [66].

## 4.2 Efficient numerical integration: the High-Performance ROM (HP-ROM)

### 4.2.1 Overview

As outlined in the introductory Chapter, approaches found in the model reduction literature that, directly or indirectly, deal with the still underexplored question of how to efficiently –-at an online computational cost independent of the dimension ${\textstyle n}$ of the underlying finite element model–- integrate reduced-order equations can be broadly classified either as interpolatory approaches [28,29,30,31,32] or Gauss-type quadrature methods [33,34]. The integration strategy proposed in the present work falls into the former category of interpolatory approaches. Recall that the basic idea in such approaches is to replace the nonaffine term in the BVP by a low-dimensional approximation. In our case, a glance at the reduced-order equilibrium equation (3.24) reveals that such “offending”, nonaffine term is the stress field –-the reduced ${\textstyle B}$-matrix ${\textstyle {{\boldsymbol {B}}^{*}}={{\boldsymbol {B}}^{*}}(\mathbf {x} )}$ is independent of the input parameter ${\textstyle {\boldsymbol {\epsilon _{M}}}}$ and hence need not be subject to approximation. The proposed integration scheme, thus, is predicated on the assumption that, not only the displacement fluctuations, but also the stress field over the cell admits an accurate, low-dimensional approximation. Numerical experiments confirm (see Chapter 5) that, luckily, this premise generally holds whenever the displacement fluctuations field itself lives in a low-dimensional space.

Let ${\textstyle {\boldsymbol {\mathit {\Psi }}}_{i}(\mathbf {x} )\in L_{2}(\Omega )^{s}}$ (${\textstyle i=1,2\ldots n_{\sigma }}$) denote a set of ${\textstyle n_{\sigma }={\mathcal {O}}(n_{u})}$ orthogonal basis functions for such low-dimensional, stress approximation space. Then, the reduced-order approximation of the stress field may be written as:

 ${\displaystyle {\boldsymbol {\sigma }}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})\approx {\boldsymbol {\sigma }}^{*}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})=\displaystyle \sum _{i=1}^{n_{\sigma }}{\boldsymbol {\mathit {\Psi }}}_{i}(\mathbf {x} )c_{i}({\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})}$
(4.2)

(notice that, in keeping with the notational convention introduced in Section 3.2, the low-dimensional approximation of the stress field is represented by attaching an asterisk to the stress symbol). In the spirit of classical polynomial quadrature (such as Newton-Cotes formulae [35]), coefficients ${\textstyle c_{i}\in \mathbb {R} ^{}}$ (${\textstyle i=1,2\ldots p}$) in Eq.(4.2) are calculated by fitting this linear expansion to the stress values computed at a set of ${\textstyle p={\mathcal {O}}(n_{u})\geq n_{\sigma }}$ pre-specified sampling points:

 ${\displaystyle {\mathcal {X}}=\{\mathbf {x} _{1}^{s},\mathbf {x} _{2}^{s},\ldots \mathbf {x} _{p}^{s}\},\;\;\;\;\mathbf {x} _{i}^{s}\in \Omega .}$
(4.3)

Approximation (4.2) becomes therefore expressible as:

 ${\displaystyle {\boldsymbol {\sigma }}^{*}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})=\displaystyle \sum _{i=1}^{p}{\boldsymbol {\mathcal {R}}}_{i}(\mathbf {x} ){\boldsymbol {\sigma }}(\mathbf {x} _{i}^{s};{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{}),}$
(4.4)

where ${\textstyle {\boldsymbol {\mathcal {R}}}_{i}(\mathbf {x} )}$ (${\textstyle i=1,2\ldots p}$) stands for the interpolation or, more generally, reconstruction operator1 at sampling point ${\textstyle \mathbf {x} _{i}^{s}}$, whereas

 ${\displaystyle {\boldsymbol {\sigma }}(\mathbf {x} _{i}^{s};{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})={\mathcal {F}}(\mathbf {x} _{i}^{s};{\boldsymbol {\epsilon _{M}}}+{{\boldsymbol {B}}^{*}}(\mathbf {x} _{i}^{s}){\boldsymbol {U}}^{*},{\boldsymbol {\xi }}_{}(\mathbf {x} _{i}^{s}))}$
(4.5)

represents the stress vector evaluated at sampling point ${\textstyle \mathbf {x} _{i}^{s}}$ through the pertinent constitutive relation ${\textstyle {\mathcal {F}}(\mathbf {x} _{i}^{s};\cdot )}$.

Substitution of the above approximation into equation Eq.(3.25) leads to:

 ${\displaystyle \displaystyle \int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}\!\!(\mathbf {x} )\,{\boldsymbol {\sigma }}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})\,d\Omega }$ ${\displaystyle \approx \displaystyle \sum _{i=1}^{p}\overbrace {\left(\int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}\!\!(\mathbf {x} ){\boldsymbol {\mathcal {R}}}_{i}(\mathbf {x} )\,d\Omega \right)} ^{{\boldsymbol {Q}}_{i}^{T}}{\boldsymbol {\sigma }}(\mathbf {x} _{i}^{s};{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})={\boldsymbol {0}}.}$
(4.6)

The bracketed integral (denoted by ${\textstyle {\boldsymbol {Q}}_{i}^{T}\in \mathbb {R} ^{n_{u}\times p}}$) in the above equation is independent of the input parameter –-the macroscopic strain ${\textstyle {\boldsymbol {\epsilon _{M}}}}$–-, and, hence, it can be entirely pre-computed offline, using, for instance, the full set of finite element Gauss points:

 ${\displaystyle {\boldsymbol {Q}}_{i}^{T}{\mathrel {\mathop {:}}}=\int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}\!\!(\mathbf {x} ){\boldsymbol {\mathcal {R}}}_{i}(\mathbf {x} )\,d\Omega \approx \displaystyle \sum _{g=1}^{n_{g}}w_{g}{{\boldsymbol {B}}^{*}}^{T}\!\!(\mathbf {x} _{g}){\boldsymbol {\mathcal {R}}}_{i}(\mathbf {x} _{g}),\;\;\;\;i=1,2\ldots p.}$
(4.7)

Introducing the above definition into Eq.(4.6) finally yields the quadrature formula for the reduced equilibrium equation:

 ${\displaystyle \displaystyle \int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}\!\!(\mathbf {x} )\,{\boldsymbol {\sigma }}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})\,d\Omega \approx \sum _{i=1}^{p}{\boldsymbol {Q}}_{i}^{T}\,{\boldsymbol {\sigma }}(\mathbf {x} _{i}^{s};{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})={\boldsymbol {0}}.}$
(4.8)

It is noteworthy that this quadrature formula requires evaluation of stresses at only ${\textstyle p={\mathcal {O}}(n_{u})}$ sampling points within the cell domain (in contrast to the standard ROM (4.1), that needs ${\textstyle n_{g}={\mathcal {O}}(n)>>n_{u}}$ Gauss points). Furthermore, inserting Eq.(4.4) into Eq.(2.14), we get:

 ${\displaystyle {\boldsymbol {\sigma _{M}}}={\dfrac {1}{V}}\displaystyle \int _{\Omega }\!\!{\boldsymbol {\sigma }}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})\,d\Omega }$ ${\displaystyle \approx \displaystyle \sum _{i=1}^{p}\overbrace {\left({\dfrac {1}{V}}\int _{\Omega }{\boldsymbol {\mathcal {R}}}_{i}(\mathbf {x} )\,d\Omega \right)} ^{{\boldsymbol {T}}_{i}}{\boldsymbol {\sigma }}(\mathbf {x} _{i}^{s};{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{})}$ ${\displaystyle =\sum _{i=1}^{p}{\boldsymbol {T}}_{i}\,{\boldsymbol {\sigma }}(\mathbf {x} _{i}^{s};{\boldsymbol {\epsilon _{M}}}+\nabla ^{s}{{\boldsymbol {u}}^{*}},{\boldsymbol {\xi }}_{}),}$
(4.9)

where

 ${\displaystyle {\boldsymbol {T}}_{i}{\mathrel {\mathop {:}}}={\dfrac {1}{V}}\int _{\Omega }{\boldsymbol {\mathcal {R}}}_{i}(\mathbf {x} )\,d\Omega \approx {\dfrac {\displaystyle \sum _{g=1}^{n_{g}}w_{g}{\boldsymbol {\mathcal {R}}}_{i}(\mathbf {x} _{g})}{\displaystyle \sum _{g=1}^{n_{g}}w_{g}}},\;\;\;\;i=1,2\ldots p.}$
(4.10)

Note that ${\textstyle {\boldsymbol {T}}_{i}}$ (${\textstyle i=1,2\ldots p}$) can also be pre-computed offline.

In summary, in projecting the cell equilibrium equation onto the displacement fluctuations, reduced space, and in additionally adopting the above described integration scheme, one automatically arrives at a reduced-order model in which the operation count –-the complexity–- in both solving the cell equilibrium equation and calculating the macroscopic stress tensor depends exclusively on the dimension ${\displaystyle n_{u}}$ of the reduced basis. We shall refer to this model as the High-Performance, Reduced-Order Model (HP-ROM), to highlight the tremendous gains in performance that affords this model over the previously described standard ROM, let alone over the full-order, finite model, discussed in Section 2.3–-in the numerical example shown in Chapter 5, we report speedup factors of above three order of magnitudes.

(1) The term interpolation usually connotes that the fit is exact at the sampling points, and this only occurs when ${\displaystyle p=n_{\sigma }}$.

## 4.3 Stress approximation space

Two crucial aspects of the integration scheme sketched in the foregoing remains to be addressed, namely, the determination of the vector space (hereafter denoted by ${\textstyle {\mathcal {V}}_{\sigma }^{apr}}$) in which the low-dimensional approximation of the stress field should lie in order to obtain an accurate and at the same time well-posed HP-ROM; and the calculation of the optimal location of the sampling or integration points at which the stress tensor is to be evaluated. Attention here and in the next Section is confined to the aspect related to the stress approximation space, while the discussion of the issue related to the selection of sampling points is deferred to Section 4.5.

### 4.3.1 The reduced-order subspace of statically admissible stresses (Vσ*)

Similarly to the problem addressed in Chapter 3 concerning the reduced basis for the displacement fluctuations, the problem of constructing a ${\textstyle {\mathcal {O}}(n_{u})}$-dimensional representation of the stress field reduces, in principle, to finding a set of orthogonal basis functions ${\textstyle \{{\boldsymbol {\mathit {\Psi }}}_{1}(\mathbf {x} ),{\boldsymbol {\mathit {\Psi }}}_{2}(\mathbf {x} )\ldots {\boldsymbol {\mathit {\Psi }}}_{n_{\sigma }}(\mathbf {x} )\}}$ (${\textstyle n_{\sigma }={\mathcal {O}}(n_{u})}$) such that its span accurately approximates the set of all possible stress solutions –-that is, the set of all statically admissible stresses. Accordingly, the procedure to compute the reduced basis for the stress field would be, mutatis mutandis, formally identical to that explained earlier for the displacement fluctuations. Firstly, finite element, stress distributions over the cell are computed for representative, input macro-strain histories (the most practical and somehow consistent choice regarding these strain trajectories is to use the same as in the computation of the displacement fluctuations snapshots). Then, the elastic/inelastic dimensionality reduction process set forth in Section 3.1.2 is applied to the resulting ensemble of stress solutions ${\textstyle \{{\boldsymbol {\sigma }}^{1}(\mathbf {x} ),{\boldsymbol {\sigma }}^{2}(\mathbf {x} )\ldots {\boldsymbol {\sigma }}^{n_{snp}}(\mathbf {x} )\}}$, in order to identify both the elastic and the essential inelastic stress modes. The space spanned by these modes will be denoted hereafter by ${\textstyle {\mathcal {V}}_{\sigma }^{*}}$ and termed the reduced-order subspace of statically admissible stresses:

 ${\displaystyle {\mathcal {V}}_{\sigma }^{*}={\textrm {span}}\{\overbrace {{\boldsymbol {\mathit {\Psi }}}_{1}(\mathbf {x} ),{\boldsymbol {\mathit {\Psi }}}_{2}(\mathbf {x} ),\ldots ,{\boldsymbol {\mathit {\Psi }}}_{m_{e}}(\mathbf {x} )} ^{\hbox{Elastic stress modes}},\overbrace {{\boldsymbol {\mathit {\Psi }}}_{m_{e}+1}(\mathbf {x} ),{\boldsymbol {\mathit {\Psi }}}_{m_{e}+2}(\mathbf {x} ),\ldots ,{\boldsymbol {\mathit {\Psi }}}_{n_{\sigma }}(\mathbf {x} )} ^{\hbox{'Essential', inelastic stress modes}}\}.}$
(4.11)

### 4.3.2 Ill-posedness of the HP-ROM

At first sight, it appears reasonable to simply construct the low-dimensional approximation ${\textstyle {\boldsymbol {\sigma }}^{*}}$ required in the proposed integration method as a linear combination of the above described stress reduced basis–- hence making ${\textstyle {\mathcal {V}}_{\sigma }^{apr}={\mathcal {V}}_{\sigma }^{*}}$–-; i.e.,

 ${\displaystyle {\boldsymbol {\sigma }}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})\approx {\boldsymbol {\sigma }}^{*}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})=\displaystyle \sum _{i=1}^{n_{\sigma }}{\boldsymbol {\mathit {\Psi }}}_{i}(\mathbf {x} )c_{i}({\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*}),}$
(4.12)

where ${\textstyle c_{i}\in \mathbb {R} ^{}}$ (${\textstyle i=1,2\ldots n_{\sigma }}$). This strategy of approximating the offending, nonaffine term in the BVP by a linear combination of pre-computed basis functions –-obtained, in turn, from samples of the nonaffine term evaluated at the solution–- has been successfully applied by several authors, with no apparent –-or at least not reported–- computational pitfalls, to a wide gamut of problems: nonlinear monotonic elliptic and nonlinear parabolic BPVs [67,29], nonlinear miscible viscous fingering in porous media [68,31], uncertainty quantification in inverse problems [69], and nonlinear heat conduction problems [32,66], to cite but a few.

However, a closer examination of the the cell equilibrium problem reveals that, in this case, this “standard” strategy proves completely fruitless, for it leads to patently ill-posed reduced-order equations. To show this, let us first substitute approximation (4.12) into Eq.(3.24):

 ${\displaystyle \displaystyle \int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}(\mathbf {x} )\,{\boldsymbol {\sigma }}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})\,d\Omega \approx \displaystyle \int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}(\mathbf {x} )\,{\boldsymbol {\sigma }}^{*}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})=}$ ${\displaystyle \displaystyle \sum _{i=1}^{n_{\sigma }}\left(\int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}(\mathbf {x} ){\boldsymbol {\mathit {\Psi }}}_{i}(\mathbf {x} )\,d\Omega \right)c_{i}({\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})={\boldsymbol {0}}.}$
(4.13)

By virtue of Eq.(3.27), the bracketed integral in the preceding equation can be rephrased as:

 ${\displaystyle \displaystyle \int _{\Omega }{{\boldsymbol {B}}^{*}}^{T}(\mathbf {x} ){\boldsymbol {\mathit {\Psi }}}_{i}(\mathbf {x} )\,d\Omega ={\boldsymbol {\Phi }}^{T}\!\!\left(\int _{\Omega }{\boldsymbol {B}}^{T}(\mathbf {x} ){\boldsymbol {\mathit {\Psi }}}_{i}(\mathbf {x} )\,d\Omega \right),\;\;\;i=1,2\ldots n_{\sigma }.}$
(4.14)

Each basis function ${\textstyle {\boldsymbol {\mathit {\Psi }}}_{i}(\mathbf {x} )}$ (${\textstyle i=1,2\ldots n_{\sigma }}$) is, by construction (see Chapter 3), a linear combination of the stress snapshots collected during the offline, finite element analysis; thus, we can write:

 ${\displaystyle {\boldsymbol {\mathit {\Psi }}}_{i}=\displaystyle \sum _{j=1}^{n_{snp}}\beta _{ij}{\boldsymbol {\sigma }}^{j}\;\;\;i=1,2\ldots n_{\sigma },}$
(4.15)

${\textstyle \beta _{ij}\in \mathbb {R} ^{}}$ being the corresponding coefficients in the linear combination. Inserting the above equation into Eq.(4.14) and considering that ${\textstyle {\boldsymbol {\sigma }}^{j}}$ (${\textstyle j=1,2\ldots n_{snp}}$) are finite element stress solutions –-and therefore fulfill the finite element equilibrium equation (2.22)–-, we finally arrive at:

 ${\displaystyle {\boldsymbol {\Phi }}^{T}\!\!\left(\int _{\Omega }{\boldsymbol {B}}^{T}(\mathbf {x} ){\boldsymbol {\mathit {\Psi }}}_{i}(\mathbf {x} )\,d\Omega \right)={\boldsymbol {\Phi }}^{T}\!\!\sum _{j=1}^{n_{snp}}\beta _{ij}\overbrace {\left(\int _{\Omega }{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}^{j}\,d\Omega \right)} ^{=0}={\boldsymbol {0}},\;\;\;i=1,2\ldots n_{\sigma },}$
(4.16)

that is, the integral (4.14) appearing in the equilibrium equation (4.13), and hence, the left-hand side of the equation itself, vanishes identically regardless of the value of the modal coefficients ${\textstyle c_{i}\in \mathbb {R} ^{}}$ (${\textstyle i=1,2\ldots n_{\sigma }}$), and therefore, regardless of the value of the reduced displacement fluctuations ${\textstyle {\boldsymbol {U}}^{*}}$–-hence the ill-posedness.

### 4.3.3 Proposed remedy: the expanded space approach

It is clear from the foregoing discussion that the root cause of the ill-posedness lies in the fact that the set of all admissible stress fields (${\textstyle {\mathcal {V}}_{\sigma }}$) forms a vector space, and, consequently, the POD stress modes ${\textstyle {\boldsymbol {\mathit {\Psi }}}_{i}\in {\mathcal {V}}_{\sigma }}$ (${\textstyle i=1,2\ldots n_{\sigma }}$) –-and any linear combination of them–- turn out to be self-equilibrated fields. Thus, for the reduced-order problem to be well-posed, the approximation space ${\textstyle {\mathcal {V}}_{\sigma }^{apr}}$ cannot be only formed by statically admissible stresses, but it must also include statically inadmissible fields –-i.e. stress functions that do not satisfy the reduced-order equilibrium equation (3.24).

One plausible route for determining a low-dimensional approximation space that embraces both statically admissible and statically inadmissible stresses might be to collect, during the offline finite element calculations, not only converged stresses, but also the unconverged ones –-i.e., those generated during the corresponding iterative algorithm–-, and then perform the POD-based dimensionality reduction over the whole ensemble of snapshots1. In the present work, however, we pursue an approach that precludes the necessity of undertaking this computationally laborious and in some aspects objectionable –-there is no guarantee that the span of selected, unconverged stress snapshots covers the entire space of statically inadmissible stresses–-process. The idea behind the employed approach was originally conceived, but not fully developed, by the authors in a recent report [36]. Here, the theory underlying such an idea is further elaborated and cast into the formalisms of functional analysis.

(1) Incidentally, this way of proceeding has the flavor of the nonlinear model reduction strategy advocated by Carlberg and co-workers [78,80], in which the reduction is carried over the linearized form of the pertinent governing equation

#### 4.3.3.1 Continuum formulation

To originate our considerations from a general standpoint, it proves convenient first to rephrase the left-hand side of the reduced-order equilibrium equation Eq.(3.24) as the action of a certain linear operator ${\textstyle {\boldsymbol {G}}:L_{2}(\Omega )^{s}\rightarrow \mathbb {R} ^{n_{u}}}$ on the stress field over the cell:

 ${\displaystyle \displaystyle \int _{\Omega }{{\boldsymbol {B}}_{i}^{*}}^{T}\!\,{\boldsymbol {\sigma }}\,d\Omega =\left\langle {{\boldsymbol {B}}_{i}^{*}\!},{\boldsymbol {\sigma }}\right\rangle _{L_{2}(\Omega )}=({\boldsymbol {G}})_{i}\;\;\;i=1,2\ldots n_{u}.}$
(4.17)

Invoking now the orthogonal decomposition of ${\textstyle L_{2}(\Omega )^{s}}$ induced by this operator, one obtains:

 ${\displaystyle L_{2}(\Omega )^{s}={\mathcal {N}}({\boldsymbol {G}})\oplus {\textrm {span}}\!\left\lbrace {{\boldsymbol {B}}_{i}^{*}\!}\right\rbrace _{i=1}^{n_{u}},}$
(4.18)

where ${\textstyle {\mathcal {N}}({\boldsymbol {G}})}$ stands for the nullspace of ${\textstyle {\boldsymbol {G}}}$. Since the cell equilibrium equation has a vanishing right-hand side term, it follows that ${\textstyle {\mathcal {N}}({\boldsymbol {G}})}$ is actually the space of statically admissible stress fields. Its orthogonal complement, ${\textstyle {\textrm {span}}\!\left\lbrace {{\boldsymbol {B}}_{i}^{*}\!}\right\rbrace _{i=1}^{n_{u}}}$, can be therefore construed as the aforementioned space of statically inadmissible stresses. The key fact here is that such a space is inherently ${\displaystyle n_{u}}$-dimensional and, thus, there is no need to perform any dimensionality reduction whatsoever over unconverged snapshots to arrive at the desired basis: the strain-displacement functions ${\textstyle \{{{\boldsymbol {B}}_{1}^{*}\!},{{\boldsymbol {B}}_{2}^{*}\!}\ldots {{\boldsymbol {B}}_{n_{u}}^{*}\!}\}}$ themselves are linearly independent (albeit not orthogonal) and can thereby serve this very purpose.

According to the preceding decomposition, any ${\textstyle {\boldsymbol {\sigma }}\in L_{2}(\Omega )^{s}}$ can be resolved as (see Figure 2):

 ${\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {\sigma }}^{ad}+{\boldsymbol {\sigma }}^{in},\;\;\;{\textrm {with}}\left\langle {{\boldsymbol {\sigma }}^{ad}},{{\boldsymbol {\sigma }}^{in}}\right\rangle _{L_{2}(\Omega )}=0,}$
(4.19)

where ${\textstyle {\boldsymbol {\sigma }}^{ad}\in {\mathcal {N}}({\boldsymbol {G}})}$ and ${\textstyle {\boldsymbol {\sigma }}^{in}\in {\textrm {span}}\!\left\lbrace {{\boldsymbol {B}}_{i}^{*}\!}\right\rbrace _{i=1}^{n_{u}}}$ stand for the statically admissible and statically inadmissible components of ${\textstyle {\boldsymbol {\sigma }}}$, respectively. Following the standard approach, the statically admissible component ${\textstyle {\boldsymbol {\sigma }}^{ad}}$ –-i.e., the stress solution we wish to calculate for a given input ${\textstyle {\boldsymbol {\epsilon _{M}}}}$–- is forced to lie in the span of the POD modes ${\textstyle {\boldsymbol {\mathit {\Psi }}}_{i}}$ (${\textstyle i=1,2\ldots n_{\sigma }}$) obtained from converged snapshots:

 ${\displaystyle {\boldsymbol {\sigma }}^{ad}\approx {\boldsymbol {\sigma }}^{*}=\sum _{i=1}^{n_{\sigma }}{\boldsymbol {\mathit {\Psi }}}_{i}c_{i}^{ad},}$
(4.20)

${\textstyle c_{i}^{ad}\in \mathbb {R} ^{}}$ (${\textstyle i=1,2\ldots n_{\sigma }}$) being the corresponding modal coefficients. The non-equilibrated component ${\textstyle {\boldsymbol {\sigma }}^{in}}$, on the other hand, resides naturally in the span of the reduced strain-displacement functions, so we can directly write–-i.e., without introducing further approximations–-:

 ${\displaystyle {\boldsymbol {\sigma }}^{in}=\sum _{i=1}^{n_{u}}{{\boldsymbol {B}}_{i}^{*}\!}c_{i}^{in},}$
(4.21)

with ${\textstyle c_{i}^{in}\in \mathbb {R} ^{}}$ (${\textstyle i=1,2\ldots n_{u}}$). The low-dimensional approximation required in the proposed integration method, denoted in what follows by ${\textstyle {\boldsymbol {\sigma }}^{ex*}}$ (the appended superscript “ex” means stress approximated in the expanded space), is finally obtained as the sum of Eq.(4.20) and Eq.(4.21) :

 ${\displaystyle {\boldsymbol {\sigma }}^{ex*}=\displaystyle \sum _{i=1}^{n_{\sigma }}{\boldsymbol {\mathit {\Psi }}}_{i}c_{i}^{ad}+\sum _{j=1}^{n_{u}}{{\boldsymbol {B}}_{j}^{*}\!}c_{j}^{in}.}$
(4.22)
 Figure 2: Expanded space approach. The stress approximation space is expanded so that it embraces, not only the span of the stress POD modes, but also the span of the reduced strain-displacement functions ${\displaystyle \{{{\boldsymbol {B}}_{1}^{*}\!},{{\boldsymbol {B}}_{2}^{*}\!}\ldots {{\boldsymbol {B}}_{n_{u}}^{*}\!}\}}$. The reduced-order cell equilibrium problem boils down to find the reduced displacement fluctuations vector ${\displaystyle {\boldsymbol {U}}^{*}}$ that makes the non-equilibrated component ${\displaystyle {\boldsymbol {\sigma }}^{in}}$ to vanish (${\displaystyle {\boldsymbol {\sigma }}^{in}({\boldsymbol {U}}^{*},{\boldsymbol {\epsilon _{M}}})={\boldsymbol {0}}}$ ).

Substituting the above approximation into the equilibrium equation, one gets:

 ${\displaystyle \displaystyle \int _{\Omega }{{\boldsymbol {B}}_{i}^{*}}^{T}\!(\mathbf {x} )\,{\boldsymbol {\sigma }}^{ex*}(\mathbf {x} ;{\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})}$ ${\displaystyle =\displaystyle \sum _{j=1}^{n_{\sigma }}\overbrace {\left(\int _{\Omega }{{\boldsymbol {B}}_{i}^{*}}^{T}\!(\mathbf {x} ){\boldsymbol {\mathit {\Psi }}}_{j}(\mathbf {x} )\,d\Omega \right)} ^{=0}c_{j}^{ad}({\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})+\sum _{j=1}^{n_{u}}\left(\int _{\Omega }{{\boldsymbol {B}}_{i}^{*}}^{T}\!(\mathbf {x} ){{\boldsymbol {B}}_{j}^{*}\!}(\mathbf {x} )\,d\Omega \right)c_{j}^{in}({\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})}$ ${\displaystyle =\sum _{j=1}^{n_{u}}\left(\int _{\Omega }{{\boldsymbol {B}}_{i}^{*}}^{T}\!(\mathbf {x} ){{\boldsymbol {B}}_{j}^{*}\!}(\mathbf {x} )\,d\Omega \right)c_{j}^{in}({\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})={\boldsymbol {0}},\;\;\;\;i=1,2\ldots n_{u}.}$
(4.23)

Since ${\textstyle \{{{\boldsymbol {B}}_{1}^{*}\!},{{\boldsymbol {B}}_{2}^{*}\!}\ldots {{\boldsymbol {B}}_{n_{u}}^{*}\!}\}}$ are linearly independent functions, it becomes immediately clear that the above equations holds only if:

 ${\displaystyle c_{j}^{in}({\boldsymbol {\epsilon _{M}}},{\boldsymbol {U}}^{*})=0,\;\;j=1,2\ldots n_{u},}$
(4.24)

i.e., if the ${\textstyle n_{u}}$ coefficients multiplying ${\textstyle {{\boldsymbol {B}}_{i}^{*}\!}\in L_{2}(\Omega )^{s}}$ (${\textstyle i=1,2\ldots n_{u}}$) are identically zero. In adopting the proposed integration approach, thus, the reduced-order cell equilibrium problem (3.24) is transformed into the problem of finding, for a given input macroscopic strain tensor ${\textstyle {\boldsymbol {\epsilon _{M}}}}$, the reduced displacement fluctuations vector ${\textstyle {\boldsymbol {U}}^{*}\in \mathbb {R} ^{n_{u}}}$ that makes the non-equilibrated component ${\textstyle {\boldsymbol {\sigma }}^{in}}$ (defined in Eq.(4.21)) to vanish.

In a nutshell, the ill-posedness exhibited by the discrete problem when adopting the standard approach of using only POD modes is eliminated by expanding the stress approximation space so that it embraces also the span of the reduced strain-displacement functions (or strain modes2) ${\textstyle {{\boldsymbol {B}}_{i}^{*}\!}\in L_{2}(\Omega )^{s}}$ (${\textstyle i=1,2\ldots n_{u}}$):

 ${\displaystyle {\mathcal {V}}_{\sigma }^{apr}={\mathcal {V}}_{\sigma }^{*}\,\oplus \,{\textrm {span}}\!\left\lbrace {{\boldsymbol {B}}_{i}^{*}\!}\right\rbrace _{i=1}^{n_{u}}={\textrm {span}}\{\overbrace {{\boldsymbol {\mathit {\Psi }}}_{1},{\boldsymbol {\mathit {\Psi }}}_{2}\ldots {\boldsymbol {\mathit {\Psi }}}_{n_{\sigma }}} ^{n_{\sigma }\;{\hbox{stress modes}}},\overbrace {{{\boldsymbol {B}}_{1}^{*}\!},{{\boldsymbol {B}}_{2}^{*}\!}\ldots {{\boldsymbol {B}}_{n_{u}}^{*}\!}} ^{n_{u}\;{\hbox{strain modes}}}\}.}$
(4.25)

#### 4.3.3.2 Discrete formulation

In typical finite element implementations, both stresses and gradients of shape functions are only calculated and stored at the Gauss points of the underlying spatial discretization. For practical reasons, thus, it proves imperative to reformulate the above explained expanded space strategy and treat both magnitudes as spatially discrete variables, defined only at such Gauss points.

The discrete counterparts of the continuously defined fields ${\textstyle {\boldsymbol {\sigma }}\in L_{2}(\Omega )^{s}}$ and ${\textstyle {{\boldsymbol {B}}_{i}^{*}\!}\in L_{2}(\Omega )^{s}}$ (${\textstyle i=1,2\ldots n_{u}}$) will be denoted by ${\textstyle {\boldsymbol {\mathcal {S}}}\in \mathbb {R} ^{n_{g}\cdot s}}$ and ${\textstyle {\boldsymbol {{\mathcal {B}}^{*}}}={\begin{bmatrix}{\boldsymbol {{\mathcal {B}}_{1}^{*}}}&{\boldsymbol {{\mathcal {B}}_{2}^{*}}}&\cdots &{\boldsymbol {{\mathcal {B}}_{n_{u}}^{*}}}\end{bmatrix}}\in \mathbb {R} ^{n_{g}\cdot s\times n_{u}}}$, and termed the global stress vector, and the global matrix of strain modes, respectively. The global stress vector ${\textstyle {\boldsymbol {\mathcal {S}}}}$ is constructed by stacking the stress vectors ${\textstyle {\boldsymbol {\sigma }}(\mathbf {x} _{g};\cdot )\in \mathbb {R} ^{s}}$ (${\textstyle g=1,2\ldots n_{g}}$) at the Gauss points of the finite element grid into a single column vector:

 ${\displaystyle {\boldsymbol {\mathcal {S}}}{\mathrel {\mathop {:}}}={\begin{bmatrix}{\boldsymbol {\sigma }}(\mathbf {x} _{1};\cdot )\\{\boldsymbol {\sigma }}(\mathbf {x} _{2};\cdot )\\\vdots \\{\boldsymbol {\sigma }}(\mathbf {x} _{n_{g}};\cdot )\end{bmatrix}}.}$
(4.26)

Similarly, the global matrix of strain modes ${\textstyle {\boldsymbol {{\mathcal {B}}^{*}}}}$ is constructed as:

 ${\displaystyle {\boldsymbol {{\mathcal {B}}^{*}}}{\mathrel {\mathop {:}}}={\begin{bmatrix}{{\boldsymbol {B}}^{*}}(\mathbf {x} _{1})\\{{\boldsymbol {B}}^{*}}(\mathbf {x} _{2})\\\vdots \\{{\boldsymbol {B}}^{*}}(\mathbf {x} _{n_{g}})\end{bmatrix}}={\begin{bmatrix}{{\boldsymbol {B}}_{1}^{*}\!}(\mathbf {x} _{1})&{{\boldsymbol {B}}_{2}^{*}\!}(\mathbf {x} _{1})&\cdots &{{\boldsymbol {B}}_{n_{u}}^{*}\!}(\mathbf {x} _{1})\\{{\boldsymbol {B}}_{1}^{*}\!}(\mathbf {x} _{2})&{{\boldsymbol {B}}_{2}^{*}\!}(\mathbf {x} _{2})&\cdots &{{\boldsymbol {B}}_{n_{u}}^{*}\!}(\mathbf {x} _{2})\\\vdots &\vdots &\vdots &\vdots \\{{\boldsymbol {B}}_{1}^{*}\!}(\mathbf {x} _{n_{g}})&{{\boldsymbol {B}}_{2}^{*}\!}(\mathbf {x} _{n_{g}})&\cdots &{{\boldsymbol {B}}_{n_{u}}^{*}\!}(\mathbf {x} _{n_{g}})\end{bmatrix}}.}$
(4.27)

With definitions (4.26) and (4.27) at hand, the standard ROM equilibrium equation (4.1) can be readily rephrased in the following compact, matrix form:

 ${\displaystyle \left\langle {{\boldsymbol {B}}_{i}^{*}\!},{\boldsymbol {\sigma }}\right\rangle _{L_{2}(\Omega )}=\displaystyle \int _{\Omega }{{\boldsymbol {B}}_{i}^{*}}^{T}\!(\mathbf {x} ){\boldsymbol {\sigma }}(\mathbf {x} ;\cdot )\,d\Omega \approx \displaystyle \sum _{g=1}^{n_{g}}w_{g}{{{\boldsymbol {B}}_{i}^{*}}^{T}\!}\!(\mathbf {x} _{g}){\boldsymbol {\sigma }}(\mathbf {x} _{g};\cdot )=0}$ ${\displaystyle \Rightarrow {\boldsymbol {{\mathcal {B}}_{i}^{*}}}^{T}{\boldsymbol {W}}{\boldsymbol {\mathcal {S}}}=0,\;\;\;\;\;\;i=1,2\ldots n_{u},}$
(4.28)

where ${\textstyle {\boldsymbol {W}}}$ is a diagonal matrix containing the weights at each Gauss point:

 ${\displaystyle {\boldsymbol {W}}{\mathrel {\mathop {:}}}={\begin{bmatrix}w_{1}{\boldsymbol {I}}&{\boldsymbol {0}}&{\boldsymbol {0}}&\cdots &{\boldsymbol {0}}\\{\boldsymbol {0}}&w_{2}{\boldsymbol {I}}&{\boldsymbol {0}}&\cdots &{\boldsymbol {0}}\\\vdots &\vdots &\vdots &\vdots &\vdots \\{\boldsymbol {0}}&{\boldsymbol {0}}&{\boldsymbol {0}}&{\boldsymbol {0}}&w_{n_{g}}{\boldsymbol {I}}\\\end{bmatrix}}^{}.}$
(4.29)

(here, ${\textstyle {\boldsymbol {I}}}$ denotes the ${\textstyle s\,{\textrm {x}}\,s}$ identity matrix). Assuming that ${\textstyle w_{g}>0}$ (${\textstyle g=1,2\ldots n_{g}}$) –-Gauss quadrature rules with negative weights are excluded from our considerations–-, one can reexpress Eq.(4.28) as:

 ${\displaystyle \left\langle {{\boldsymbol {B}}_{i}^{*}\!},{\boldsymbol {\sigma }}\right\rangle _{L_{2}(\Omega )}\approx \left\langle {\boldsymbol {{\mathcal {B}}_{i}^{*}}},{\boldsymbol {\mathcal {S}}}\right\rangle _{\boldsymbol {W}}=0,\;\;\;\;\;\;i=1,2\ldots n_{u},}$
(4.30)

where ${\textstyle \left\langle {\cdot },{\cdot }\right\rangle _{\boldsymbol {W}}}$ symbolizes the following ${\textstyle \mathbb {R} ^{n_{g}\cdot s}}$ inner product:

 ${\displaystyle \left\langle {a},{b}\right\rangle _{\boldsymbol {W}}={a}^{T}{\boldsymbol {W}}{b},\;\;\;\;{a},{b}\in \mathbb {R} ^{n_{g}\cdot s}.}$
(4.31)

Equation (4.30) reveals that any statically admissible global stress vector ${\textstyle {\boldsymbol {\mathcal {S}}}}$ is orthogonal to the global strain modes vectors ${\textstyle {\boldsymbol {{\mathcal {B}}_{i}^{*}}}}$ (${\textstyle i=1,2\ldots n_{u}}$) in the sense of the inner product induced by ${\textstyle {\boldsymbol {W}}}$. In other words, in approximating the integral of the internal forces by Gauss quadrature, the ${\textstyle L_{2}}$ orthogonality condition translates into orthogonality in the sense of Eq.(4.31). From the point of view of numerical implementation, however, it is preferable to cast the equilibrium condition and the subsequent developments in terms of the standard euclidean scalar product in ${\textstyle \mathbb {R} ^{n_{g}\cdot s}}$ –-working with the inner product defined by Eq.(4.31) is somehow a nuisance and complicates unnecessarily the involved algebra. This can be achieved by inserting the Cholesky decomposition of the weights matrix

 ${\displaystyle {\boldsymbol {W}}={\boldsymbol {W}}^{1/2}{\boldsymbol {W}}^{1/2},\;\;\;\;({\boldsymbol {W}}_{ii}^{1/2}{\mathrel {\mathop {:}}}={\sqrt {w_{i}}})}$
(4.32)

into Eq.(4.26):

 ${\displaystyle {\boldsymbol {{\mathcal {B}}^{*}}}^{T}{\boldsymbol {W}}{\boldsymbol {\mathcal {S}}}={({\boldsymbol {{\mathcal {B}}^{*}}}^{T}{\boldsymbol {W}}^{1/2})}({\boldsymbol {W}}^{1/2}{\boldsymbol {\mathcal {S}}})={\boldsymbol {0}}.}$
(4.33)

Defining now the weighted global stress vector and weighted matrix of strain modes as

 ${\displaystyle {\boldsymbol {\Sigma }}{\mathrel {\mathop {:}}}={\boldsymbol {W}}^{1/2}{\boldsymbol {\mathcal {S}}}={\begin{bmatrix}{\sqrt {w_{1}}}{\boldsymbol {\sigma }}(\mathbf {x} _{1};\cdot )\\{\sqrt {w_{2}}}{\boldsymbol {\sigma }}(\mathbf {x} _{2};\cdot )\\\vdots \\{\sqrt {w_{n_{g}}}}{\boldsymbol {\sigma }}(\mathbf {x} _{n_{g}};\cdot )\end{bmatrix}},}$
(4.34)

and

 ${\displaystyle {\mathbb {B} ^{*}}{\mathrel {\mathop {:}}}={\boldsymbol {W}}^{1/2}{\boldsymbol {{\mathcal {B}}^{*}}}={\begin{bmatrix}{\sqrt {w_{1}}}{{\boldsymbol {B}}^{*}}(\mathbf {x} _{1})\\{\sqrt {w_{2}}}{{\boldsymbol {B}}^{*}}(\mathbf {x} _{2})\\\vdots \\{\sqrt {w_{n_{g}}}}{{\boldsymbol {B}}^{*}}(\mathbf {x} _{n_{g}})\end{bmatrix}}}$
(4.35)

respectively, and inserting these definitions into Eq.(4.26), one finally arrives at:

 ${\displaystyle {\mathbb {B} ^{*}}^{T}{\boldsymbol {\Sigma }}={\boldsymbol {0}},}$
(4.36)

or equivalently:

 ${\displaystyle {\mathbb {B} _{i}^{*}}^{T}{\boldsymbol {\Sigma }}=0,\;\;\;\;\;\;\;i=1,2\ldots n_{u},}$
(4.37)

which shows that any statically admissible weighted stress vector is orthogonal, in the sense of the standard euclidean inner product, to the weighted strain modes ${\textstyle {\mathbb {B} _{i}^{*}}^{T}}$ (${\textstyle i=1,2\ldots n_{u}}$).

Comparing Eq.(4.36) with Eq.(4.17), it becomes clear that ${\textstyle {\mathbb {B} ^{*}}^{T}}$ plays the same role as operator ${\textstyle {\boldsymbol {G}}}$ in Eq.(4.17). In analogy with Eq.(4.18), thus, we can write

 ${\displaystyle \mathbb {R} ^{n_{g}\cdot s}={\mathcal {N}}({\mathbb {B} ^{*}}^{T})\oplus {\textrm {Range}}({\mathbb {B} ^{*}})}$
(4.38)

where ${\textstyle {\mathcal {N}}({\mathbb {B} ^{*}}^{T})}$ and ${\textstyle {\textrm {Range}}({\mathbb {B} ^{*}})}$ denote the null space and the range (or column space) of ${\textstyle {\mathbb {B} ^{*}}^{T}}$ and ${\textstyle {\mathbb {B} ^{*}}}$, respectively, and consequently decompose any ${\textstyle {\boldsymbol {\Sigma }}\in \mathbb {R} ^{n_{g}\cdot s}}$ as

 ${\displaystyle {\boldsymbol {\Sigma }}={\boldsymbol {\Sigma ^{ad}}}+{\boldsymbol {\Sigma ^{in}}}}$
(4.39)

with ${\textstyle {\boldsymbol {\Sigma ^{ad}}}\in {\mathcal {N}}({\mathbb {B} ^{*}}^{T})}$ and ${\textstyle {\boldsymbol {\Sigma ^{in}}}\in {\textrm {Range}}({\mathbb {B} ^{*}})}$. As in the continuous case (see Eq.(4.20)), the statically admissible component ${\textstyle {\boldsymbol {\Sigma ^{ad}}}}$ is now approximated by a linear combination of POD basis vectors obtained from converged stress snapshots (the methodology for obtaining these modes using the SVD is thoroughly explained in Appendix B.2):

 ${\displaystyle {\boldsymbol {\Sigma ^{ad}}}\approx {\boldsymbol {\Sigma }}^{*}=\sum _{i=1}^{n_{\sigma }}c_{i}^{ad}{\boldsymbol {\Psi }}_{i}={\boldsymbol {\Psi }}{\boldsymbol {c}}^{ad},}$
(4.40)

where

 ${\displaystyle {\boldsymbol {\Psi }}={\begin{bmatrix}{\boldsymbol {\Psi }}_{1}&{\boldsymbol {\Psi }}_{2}&\cdots &{\boldsymbol {\Psi }}_{n_{\sigma }}\end{bmatrix}}}$
(4.41)

denotes the (weighted) stress basis matrix and ${\textstyle {\boldsymbol {c}}^{ad}\in \mathbb {R} ^{n_{\sigma }}}$ stands for the vector of modal coefficients associated to such a basis matrix. Likewise, since the non-equilibrated component ${\textstyle {\boldsymbol {\Sigma ^{in}}}}$ pertains to the column space of ${\textstyle {\mathbb {B} ^{*}}}$, we can directly write

 ${\displaystyle {\boldsymbol {\Sigma ^{in}}}={\mathbb {B} ^{*}}{\boldsymbol {c}}^{in},}$
(4.42)

where ${\textstyle {\boldsymbol {c}}^{in}\in \mathbb {R} ^{n_{u}}}$. The low-dimensional (weighted) stress vector ${\textstyle {\boldsymbol {\Sigma }}^{ex*}}$ required in the proposed integration method is finally obtained as the sum of Eq.(4.42) and Eq.(4.40).

 ${\displaystyle {\boldsymbol {\Sigma }}\approx {\boldsymbol {\Sigma }}^{ex*}={\boldsymbol {\Psi }}{\boldsymbol {c}}^{ad}+{\mathbb {B} ^{*}}{\boldsymbol {c}}^{in},}$
(4.43)

or in a more compact format:

 ${\displaystyle {\boldsymbol {\Sigma }}^{ex*}={\boldsymbol {\Psi }}^{ex}{\boldsymbol {c}}.}$
(4.44)

where

 ${\displaystyle {\boldsymbol {\Psi }}^{ex}{\mathrel {\mathop {:}}}={\begin{bmatrix}{\boldsymbol {\Psi }}&{\mathbb {B} ^{*}}\end{bmatrix}},}$
(4.45)

and

 ${\displaystyle {\boldsymbol {c}}={\begin{bmatrix}{\boldsymbol {c}}^{ad}\\{\boldsymbol {c}}^{in}\end{bmatrix}}.}$
(4.46)

The matrix ${\textstyle {\boldsymbol {\Psi }}^{ex}\in \mathbb {R} ^{n_{g}\cdot s\times (n_{u}+n_{\sigma })}}$ defined by Eq.(4.45) will be hereafter called the expanded basis matrix for the (weighted) stresses, whereas ${\textstyle {\boldsymbol {c}}\in \mathbb {R} ^{n_{\sigma }+n_{u}}}$ will be correspondingly termed the expanded vector of modal coefficients. Inserting approximation (4.43) into Eq.(4.36), and considering that ${\textstyle {\mathbb {B} ^{*}}^{T}{\boldsymbol {\Psi }}={\boldsymbol {0}}}$ and that ${\textstyle {\mathbb {B} ^{*}}^{T}}$ is a full rank matrix, one finally arrives at the same equilibrium condition derived in the continuum case (see Eq. 4.24):

 ${\displaystyle {\boldsymbol {c}}^{in}({\boldsymbol {U}}^{*},{\boldsymbol {\epsilon _{M}}})={\boldsymbol {0}}.}$
(4.47)

Once the above equation is solved for ${\textstyle {\boldsymbol {U}}^{*}}$, the desired equilibrated stress vector ${\textstyle {\boldsymbol {\Sigma }}^{*}}$ is obtained by evaluating Eq.(4.40):

 ${\displaystyle {\boldsymbol {\Sigma }}^{*}={\boldsymbol {\Psi }}{\boldsymbol {c}}^{ad}({\boldsymbol {U}}^{*},{\boldsymbol {\epsilon _{M}}}).}$
(4.48)

(2) Indeed, functions ${\displaystyle {{\boldsymbol {B}}_{i}^{*}\!}\in L_{2}(\Omega )^{s}}$ (${\displaystyle i=1,2\ldots n_{u}}$) can be viewed as fluctuating strain modes, since they are the symmetric gradient of the displacement fluctuation modes, see Eq. 3.27.

## 4.4 Determination of modal coefficients

The next step in the development of the proposed integration scheme is to deduce closed-form expressions for the vectors of modal coefficients ${\textstyle {\boldsymbol {c}}^{ad}\in \mathbb {R} ^{n_{\sigma }}}$ and ${\textstyle {\boldsymbol {c}}^{in}\in \mathbb {R} ^{n_{u}}}$ in terms of the stress values computed at a set of ${\textstyle p={\mathcal {O}}(n_{u})}$ pre-specified sampling points (to be chosen among the set of Gauss points of the underlying finite element mesh). To this end, we need first to introduce some notation and terminology.

### 4.4.1 Gappy vectors

Let ${\textstyle {\mathcal {I}}=\{{\mathcal {I}}_{1},{\mathcal {I}}_{2}\ldots {\mathcal {I}}_{p}\}\subset \{1,2\cdots n_{g}\}}$ denote the set of indices of sampling points. Notationally, we write ${\textstyle {\hat {\boldsymbol {\Sigma }}}_{({\mathcal {I}})}\in \mathbb {R} ^{p\cdot s}}$ to designate the subvector of ${\textstyle {\boldsymbol {\Sigma }}}$ containing the rows associated to these sampling points; viz.:

 ${\displaystyle {\hat {\boldsymbol {\Sigma }}}_{({\mathcal {I}})}{\mathrel {\mathop {:}}}={\begin{bmatrix}{\sqrt {w_{{\mathcal {I}}_{1}}}}{\boldsymbol {\sigma }}(\mathbf {x} _{{\mathcal {I}}_{1}},\cdot )\\{\sqrt {w_{{\mathcal {I}}_{2}}}}{\boldsymbol {\sigma }}(\mathbf {x} _{{\mathcal {I}}_{2}},\cdot )\\\vdots \\{\sqrt {w_{{\mathcal {I}}_{p}}}}{\boldsymbol {\sigma }}(\mathbf {x} _{{\mathcal {I}}_{p}},\cdot )\end{bmatrix}}.}$
(4.49)

(When confusion is not apt to arise, the parenthetical subscript indicating the set of sampling indices will be dropped, and we shall simply write ${\textstyle {\hat {\boldsymbol {\Sigma }}}}$). It proves conceptually advantageous to regard this restricted or “gappy” –-a terminology that goes back to the work of Everson et al. [70]–- stress vector ${\textstyle {\hat {\boldsymbol {\Sigma }}}_{({\mathcal {I}})}}$ as the result of the application of a certain boolean operator ${\textstyle {\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}:\mathbb {R} ^{n_{g}\cdot s}\rightarrow \mathbb {R} ^{p\cdot s}}$ over the full vector ${\textstyle {\boldsymbol {\Sigma }}}$:

 ${\displaystyle {\hat {\boldsymbol {\Sigma }}}={\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}{\boldsymbol {\Sigma }}.}$
(4.50)

We call ${\textstyle {\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}}$ the selection operator associated to sampling indices ${\textstyle {\mathcal {I}}}$. This operator can be of course applied to any ${\textstyle {Y}\in \mathbb {R} ^{n_{g}\cdot s\times z}}$ (${\textstyle z\geq 1,z\in \mathbb {N} ^{}}$). For instance, the restricted matrix of weighted strain modes would be defined as:

 ${\displaystyle {\mathbb {\hat {B}} ^{*}}{\mathrel {\mathop {:}}}={\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}{\mathbb {B} ^{*}}={\begin{bmatrix}{\sqrt {w_{{\mathcal {I}}_{1}}}}{{\boldsymbol {B}}^{*}}(\mathbf {x} _{{\mathcal {I}}_{1}})\\{\sqrt {w_{{\mathcal {I}}_{2}}}}{{\boldsymbol {B}}^{*}}(\mathbf {x} _{{\mathcal {I}}_{2}})\\\vdots \\{\sqrt {w_{{\mathcal {I}}_{p}}}}{{\boldsymbol {B}}^{*}}(\mathbf {x} _{{\mathcal {I}}_{p})})\end{bmatrix}}}$
(4.51)

Furhtermore, it is straighforward to show that

 ${\displaystyle {\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}{\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}^{T}={\boldsymbol {I}},}$
(4.52)

(here ${\textstyle {\boldsymbol {I}}}$ is the ${\textstyle (n_{g}\cdot s)}$x${\textstyle (n_{g}\cdot s)}$ identity matrix) and that

 ${\displaystyle {\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}({A}{Y})=({\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}{A}{\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}^{T})({\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}{Y})}$
(4.53)

for any ${\textstyle {A}\in \mathbb {R} ^{n_{g}\cdot s\times n_{g}\cdot s}}$ and ${\textstyle {Y}\in \mathbb {R} ^{n_{g}\cdot s\times z}}$.

### 4.4.2 Least-squares fit

In the spirit of classical polynomial quadrature, such as Newton-Cotes formulae [35], the modal coefficients ${\textstyle {\boldsymbol {c}}^{ad}\in \mathbb {R} ^{n_{\sigma }}}$ and ${\textstyle {\boldsymbol {c}}^{in}\in \mathbb {R} ^{n_{u}}}$ are determined by fitting the low-dimensional approximation (4.43) to the weighted stresses calculated at the pre-specified sampling points. It should be noticed that, the variable subject to approximation –-the stress–- being a vector-valued function, the total number of discrete points to be fitted does not coincide with the number of spatial sampling points (${\textstyle p}$), but rather is equal to the product of such a number times the number of stress components (${\textstyle s}$). The well-posedness of the fitting problem, thus, demands that:

 ${\displaystyle p\cdot s\geq n_{\sigma }+n_{u}}$
(4.54)

i.e., the number of discrete points must be equal or greater than the number of parameters to be adjusted. For the equality to hold, both ${\textstyle n_{\sigma }+n_{u}}$ and ${\textstyle p}$ have to be multiple of ${\textstyle s}$; thus, an exact fit is in general not possible for arbitrary values of ${\textstyle n_{\sigma }}$ and ${\textstyle n_{u}}$, and recourse to an approximate fit is to be made. In this respect, we follow here the standard approach of using a least-squares, best-fit criterion, i.e., minimization of the squares of the deviations between “observed” (${\textstyle {\hat {\boldsymbol {\Sigma }}}}$) and fitted (${\textstyle {\boldsymbol {\hat {\Sigma }}}^{ex*}={\boldsymbol {\hat {\Psi }}}{\boldsymbol {a}}+{\mathbb {\hat {B}} ^{*}}{\boldsymbol {b}}}$) values (in our context, “observed” signifies “calculated through the pertinent constitutive equation”). This minimization problem can be stated as:

 ${\displaystyle {\boldsymbol {c}}={\begin{bmatrix}{\boldsymbol {c}}^{ad}\\{\boldsymbol {c}}^{in}\end{bmatrix}}={\textrm {arg}}\;{\underset {{\boldsymbol {a}}\in \mathbb {R} ^{n_{\sigma }},{\boldsymbol {b}}\in \mathbb {R} ^{n_{u}}}{\textrm {min}}}\,{\Vert {\hat {\boldsymbol {\Sigma }}}-\left({\boldsymbol {\hat {\Psi }}}{\boldsymbol {a}}+{\mathbb {\hat {B}} ^{*}}{\boldsymbol {b}}\right)\Vert }}$
(4.55)

where ${\textstyle \Vert \cdot \Vert }$ stands for the standard euclidean norm. Let ${\textstyle {\boldsymbol {{\hat {\Psi }}^{ex}}}={\boldsymbol {\mathcal {P}}}_{({\mathcal {I}})}{\boldsymbol {\Psi }}^{ex}=[{\boldsymbol {\hat {\Psi }}}\;\,{\mathbb {\hat {B}} ^{*}}]}$ be the gappy expanded basis matrix, and suppose that the sampling indices ${\textstyle {\mathcal {I}}}$ have been chosen so that ${\textstyle {\boldsymbol {{\hat {\Psi }}^{ex}}}}$ has full rank, i.e.:

 ${\displaystyle {\textrm {rank}}({\boldsymbol {{\hat {\Psi }}^{ex}}})={\textrm {rank}}([{\boldsymbol {\hat {\Psi }}}\;\,{\mathbb {\hat {B}} ^{*}}])=n_{\sigma }+n_{u}.}$
(4.56)

Then, it can be shown (see, for instance, Ref. [71]) that the solution of this standard, least-squares problem is provided by the following vector of coefficients:

 ${\displaystyle {\boldsymbol {c}}={\begin{bmatrix}{\boldsymbol {c}}^{ad}\\{\boldsymbol {c}}^{in}\end{bmatrix}}={\boldsymbol {\hat {\Psi }}}^{ex^{\dagger }}{\hat {\boldsymbol {\Sigma }}},}$
(4.57)

where

 ${\stackrel{\mathbf{^}}{\mathbf{\Psi }}}^{e{x}^{†}}:=\stackrel{}{\stackrel{}{\left({{\stackrel{\mathbf{^}}{\mathbf{\Psi }}}^{\mathbit{e}\mathbit{x}}}^{T}{\stackrel{\mathbf{^}}{\mathbf{\Psi }}}^{\mathbit{e}\mathbit{x}}{\right)}^{-1}}}$