Wavelet based reduced order models for microstructural analyses

This paper proposes a novel method to accurately and efficiently reduce a microstructural mechanical model using a wavelet based discretisation. The model enriches a standard reduced order modelling (ROM) approach with a wavelet representation. Although the ROM approach reduces the dimensionality of the system of equations, the computational complexity of the integration of the weak form remains problematic. Using a sparse wavelet representation of the required integrands, the computational cost of the assembly of the system of equations is reduced significantly. This wavelet-reduced order model (W-ROM) is applied to the mechanical equilibrium of a microstructural volume as used in a computational homogenisation framework. The reduction technique however is not limited to micro-scale models and can also be applied to macroscopic problems to reduce the computational costs of the integration. For the sake of clarity, the W-ROM will be demonstrated using a one-dimensional example, providing full insight in the underlying steps taken.


Introduction
An increasing amount of engineering applications rely on materials with complex microstructures to achieve desired material properties that are tailored to the application [25]. Computational homogenisation provides an accurate and efficient framework to investigate the macroscopic properties arising from the microstructure using the Hill-Mandel conditions [22] to establish the scale transition. An overview of recent advances in computational homogenisation is presented by Geers et al. [16].
internal variables describing the constitutive behaviour of the material on each integration point in the model imposes strong requirements on the available memory.
Following Hernández et al. [20], the hyper-reduction methods can be divided into two classes. On the one hand the nodal vector approximations that approximate the integral by introducing a global basis to represent the integrand. This allows the basis to be evaluated only once (offline), thereby reducing costs of the subsequent evaluations of the integral. This basis is weighted using coefficients that minimise the interpolation error between the nonlinear integrand and its modal approximation on a set of sampled points in the least-squares sense. Examples of this class are DEIM and MPE [2,5]. On the other hand the integral (quadrature) approaches approximate the integral using a reduced set of integration points which have empirically determined weights in the methods ECM or ECSW [12,20]. A comparison of the Empirical Interpolation Method, in particular High-Performance Reduced Order Modelling [21], and the Empirical Cubature Method [20] for solving micromechanical equilibrium problems is presented in [38].
Besides Reduced Order Models, the computational costs have also been reduced using numerically efficient solvers, e.g. the Fast Fourier Transform [30] or wavelet bases to reduce the number of equations by projecting a fully discretised multi-scale problem onto a lower dimensional space, e.g. [3,9,17,19]. Some alternative wavelet-based reduction methods can be found in literature. The wavelet-based MOR [14] has been proposed to provide an alternative subspace, for those problems where there is no time to build a POD basis or when a global POD basis cannot adequately represent the local behaviour. However, this alternative approach does not attain significant compression ratios that are characteristic for POD. In this work, we will therefore still rely on a POD for the first reduction.
In this paper, the dimensionality of the problem is first tackled using the classical Reduced Order Modelling approach. Next, a novel integration scheme is proposed to limit the number of function evaluations by adaptively selecting the quadrature points using a sparse wavelet representation of the integrand. Due to the hierarchical nature of the wavelet bases, local refinement comes naturally with a wavelet representation. This approach can be considered as a form of multi-resolution analysis (MRA) used to perform the local refinement similar to Meyer [27] and Mallat [24]. The multi-resolution analysis provides control over the approximation accuracy of local phenomena [3] using a pre-defined tolerance.
A combined approach of ROM and MRA reduces the computational costs of both the assembly and solution of the microstructural models. Furthermore, it requires only a tolerance to indicate the required accuracy of the approximation, allowing to bound introduced errors. The introduction of an a priori determined reduced integration scheme for the high-dimensional parameter space is thereby omitted, greatly reducing the dimensionality of the snapshot space.
This paper is outlined as follows. First, the mechanical problem and the corresponding standard Reduced Order Model are introduced in Sects. 2 and 3 respectively. In Sect. 4, a one-dimensional mechanical model is introduced, aiming for a comprehensive analysis of the proposed method. Evidently, a 1D problem may not be as convincing in reproducing a state not captured in the snapshot space as a 2D problem might. Therefore, more convincing examples will be provided in 2D in forthcoming work. The section shortly outlines Reduced Order Modelling, after which the wavelet reduced model using MRA is introduced. As no assumptions are made with respect to the spatial dimensionality of the problem and multi-dimensional adaptive wavelet transforms are available in literature [32], the method can be readily extended to two or three dimensional problems. The presentation in a 1D context however, allows the reader to fully comprehend the underlying principles and added value of the method. The extension to a wavelet reduced multi-dimensional model will be discussed in a following publication. Section 5 presents the results of the wavelet reduced analysis of two onedimensional elasto-plastic microstructural models, whereby the accuracy and the reduction in computational costs are discussed. The paper closes with the conclusions in Sect. 6.

Notation
First the mechanical model is described in general context using a tensorial notation. Vectors and tensors are typeset using a boldface symbol, e.g. w, ε, σ etc. Spaces are noted using a calligraphic font S and columns and matrices are denoted by single and double underlined symbols, e.g. column vector c and matrix M.
Single and double contractions between vector or tensor valued quantities are denoted with a dot or a double dot respectively, e.g. a · b = s and ∇q(x) : σ = s where s represents a scalar quantity. Furthermore all basis functions N (x) and R(x) are denoted in uppercase. The corresponding coefficients that The grid point of the corresponding scaling function on a dyadic wavelet grid is denoted by x j i in the wavelet coefficients, the notation d j i [ f (x)] is introduced where i and j are the index and the level of the wavelet function. The integrands occurring in the weak form of the linear momentum balance are denoted using Greek symbols, e.g. the force and stiffness integrands are denoted by ϕ and κ respectively. A short-hand notation is employed for time discretised parameters, such as the history parameters. The time-step on which the parameter is sampled is denoted by a superscript t for the current time step and t + Δt for the next time step, e.g. ξ t = ξ (t) and ξ t+Δt = ξ (t + Δt).

Mechanical model
To outline the principles of W-ROM, this work considers a two-scale model similar to those presented in [13,16,29]. The model is comprised of a macro-mechanical model and a micro-mechanical model representing the topological and material information on the corresponding macro-and microscale. The subscript M and m are introduced to distinguish between macro-and micro-scale quantities respectively. The loading in the macro-scale model is transferred to the microscale model via periodic boundary conditions and the stress state in the macro-scale model is given by the volume averaged stress in the micro-scale model, conforming to the Hill-Mandel condition [23]. This procedure is schematically depicted in Fig. 1.
The reduced order approach is illustrated using a microscale mechanical model depicted in Fig. 2, where a Lagrangian description is employed with x the position vector of a material point x(t) at time t, V and S denote the microstructural domain and its boundary respectively and O is the origin. Neglecting the inertia and body-forces, the linear momentum balance for this problem is given by where σ is the Cauchy stress tensor and ξ are the history parameters describing the local material state. The model is formulated in a small strain framework. The local micro-scale strain tensor ε consists of a contribution of the macroscopic strain ε M and the symmetric gradient of the micro-fluctuations w(x). This leads to the following definition of the micro-scale strain ε(ε M , x, t)

Weak formulation
After multiplication of Eq. (1) with weighting function q(x) and integration by parts, the standard Bubnov-Galerkin weak form is obtained where V is the domain of the microstructural model with boundary S, on which an outward pointing unit normal n is defined. The traction t is defined as t = n·σ . The internal and external forces are denoted by f int (w) and f ext respectively. The micro-fluctuation field w(x) is constrained using periodic boundary conditions at the boundary S and the macroscopic strain ε M (t) results from the macro-scale kinematics. Using Hill-Mandel, the macroscopic stress is obtained by volume averaging the microstructural stress σ (ε, ξ), i.e.
where |V| denotes the volume of the microstructure.

Spatial discretisation
The microstructural model is discretised using a standard Lagrangian finite element basis. The weighting and trial functions q and w in the weak form of the linear momentum balance (3) are approximated using the discretised weighting and trial functions q h and w h respectively.
where N (x) is a set of Lagrangian interpolation functions. After substitution of the discretised weighting and trial functions, the finite element problem can be solved for different macroscopic strains ε M (t) using the Newton-Raphson procedure.

Reduced order modelling
To reduce the dimensionality of the finite element discretisation, the microstructural model is reformulated using a set of reduced basis functions. In the traditional FE discretisation, the spatial accuracy of the model is proportional to h − p for p > 0, where h and p are the element size and order respectively. The number of physical deformation modes present in the microstructure is often much lower than the number of FE shape functions required to accurately capture the deformation field. Therefore, the number of degrees of freedom can be reduced using a set of global shape functions which are sufficiently detailed to capture the local phenomena occurring in the microstructural model. The number of degrees of freedom required for the discretisation is thereby no longer proportional to the spatial accuracy of the basis.

Proper orthogonal decomposition
The applied Reduced Order Modelling technique makes use of the Proper Orthogonal Decomposition (POD) by Pearson [33] and Schmidt [36] to extract essential physical modes from a set of solutions constructed using the original discretisation. Snapshots of the discretised micro-fluctuation coefficients w are collected in a snapshot-matrix X. The snapshot matrix is then decomposed into proper orthogonal modes and corresponding eigenvalues λ using the POD. The normalised eigenvalues represent the contribution of the corresponding mode to the snapshot matrix and are used to truncate the reduced basis to V , consisting of the first n m modes with a contribution above the predetermined tolerance δ RB . This process is schematically depicted in Fig. 3.

Construction of the reduced basis
The snapshot matrix X is obtained by collecting snapshots of the micro-fluctuation coefficients w under n s different macroscopic strains ε M as columns in a snapshot matrix X = w(ε 1 M ), w(ε 2 M ), . . . , w(ε n s M ) , after which the POD yields the modal coefficients v j and corresponding eigenvalues λ j . The reduced basis functions R(x) are then constructed out of a linear combination of the Lagrangian basis The weighting and trial functions are then discretised using the reduced basis functions only.
Substitution of the reduced weighting and trial functions (8) into the weak form (3) yields the ROM [41]: From the equations it is clear that the integration can not be performed a priori since the internal forcef int (W ) is nonlinearly dependent on the reduced micro-fluctuations. The integration scheme requires the stresses σ at every integration point in the microstructure. Therefore the integration procedure is still proportional to the number of elements used to discretise the microstructural problem and the reduction in memory usage and floating point operations (FLOPs) required to assemble and solve the microstructural problem is only marginal. This problem can be resolved using, for example, hyperreduction techniques such as HP-ROM [21] or ECM [20]. In the first approach, not only the displacement-field but also the stress-field is projected onto a reduced basis. The modal contributions are determined in a least-squares sense by sampling the stresses in a reduced set of integration points. The latter approach reduces the integration scheme itself directly by choosing a subset of integration points and optimising the quadrature weights based on snapshots of the integrands. Both schemes rely on snapshots of either the stress-field or the complete integrand to reduce the integration.

Wavelet-Reduced Order Model
The Wavelet-Reduced Order Modelling technique will be demonstrated on a one-dimensional microstructural model. The 1D setting provides a transparent view of all underlying principles and implementation aspects, which assists the reader in getting a full comprehensive view of the method proposed. To obtain the 1D model, all previous equations are simplified to scalar expressions, e.g x becomes x, ε becomes ε, σ becomes σ etc. The homogenisation of the 1D micro-structural problem leads to a homogenised force instead of a homogenised stress. This model with length and spatially varying cross-sectional area A(x) is schematically depicted in Fig. 4. An elasto-plastic model is used to describe the material behaviour. The microstructure is loaded with a macroscopic strain ε M = Δ with Δ the increase in length of the microstructure.
To reduce the computational integration costs, the integrands in the reduced order model are approximated using wavelets. In this case, the approximation is performed using Deslauriers-Dubuc interpolating wavelets [7,8,10]. This wavelet family was chosen because it yields a compact interpolating scheme. Furthermore the basis functions are smooth and therefore compatible with the smooth fields often present in mechanical problems.
The W-ROM relies on several wavelet techniques, such as wavelet synthesis, wavelet analysis, MRA and data compression. These concepts are briefly explained in Appendix A. The interested reader is also referred to the book by Goedecker [18].

Wavelet representation of the integrand
The Reduced Order Model requires an integration step to determine the reduced internal forcesf int and the corresponding tangent stiffnessK . To approximate these integrals efficiently, the integrands are projected on the wavelet grid using MRA. The MRA automatically determines the sampling points for interpolation, and the integration is performed on the wavelet approximation of the reduced integrand. Periodic Node The integrands ϕ(ε, ξ ) and κ(ε, ξ ) for the internal force and stiffness are defined in (12) and (13) respectively.
For the present analysis, periodicity of the microstructural model is used to approximate the function values at the boundary of the domain. Note, however, that wavelets are not limited to periodic boundary conditions. For the treatment of general boundary conditions using a wavelet basis, the reader is referred to Vasilyev et al. [40].
In order to approximate the integrand ϕ(ε, ξ ) of the internal force using wavelets, the gradients of the reduced basis functions ∂ x R(x) are required on the points in the dyadic wavelet grid. The gradients are found by sampling the finite element discretised gradient at every dyadic grid point in the wavelet basis using an inverse mapping of the physical to the isoparametric coordinate. Both discretisations are schematically depicted in Fig. 5.
In this work, second order Lagrangian finite elements are used to obtain the full-order finite element solution yielding a piecewise linear approximation of the strain field. This finite element basis is sufficiently smooth to be sampled directly using the Deslauriers-Dubuc interpolating wavelets. When discontinuous strain fields are required to capture the homogenised behaviour accurately, a discontinuous wavelet transform [39] could be used to interpolate the stress and strain fields.

Multi-resolution wavelet approximation
The integrands ϕ(ε, ξ ) and κ(ε, ξ ) are projected onto Deslauriers-Dubuc interpolating wavelet basis. The basis is constructed out of multiple grid levels j = 0, 1, . . . , j max . Each level contains a set of scaling functions φ j i (x) and wavelets ψ j i (x) with index i = 0, 1, . . . , 2 j n and grid level j = 0, 1, . . . , j max . The scaling functions and wavelets are derived from the scaling function φ(x) and mother wavelet ψ(x) by scaling their width by a factor 1/2 at each increasing level and translating them with a step of the grid size Δx j corresponding to the level j. The scaling function and wavelet coefficients are found by projecting the integrands onto the corresponding scaling function φ respectively. Note that the Deslauriers-Dubuc scaling functions are interpolating. Therefore, the scaling function coefficient is given directly by the function values sampled on the wavelet grid point, i.e. s Using the reduced basis gradient ∂ x R(x j i ) sampled on the dyadic wavelet grid, the local strain in each grid point is computed using ε are discretised using the same multi-resolution sparse wavelet grid. This allows for the wavelet interpolation of the history parameters on newly added grid points avoiding the need to store the history coefficients of every dyadic wavelet grid point. Using the strain ε j i and the history ξ j i , the stress σ are computed on the wavelet grid. This enables a MRA of the internal force integrand ϕ(ε, ξ ) associated with the reduced micro-fluctuations W on the wavelet-grid. Dur-ing the MRA, the analysis of the grid points is conducted on each level in the hierarchy. The analysis of subsequent levels by evaluating and wavelet transforming the neighbouring points of current level grid points enables the error control, which other hyper-reduction methods ignore.
Note that for an accurate approximation of the internal force and tangent stiffness integrands the initial wavelet grid needs to be sufficiently fine such that the adaptive refinement scheme detects regions requiring local refinement. In principle the initial grid should be as coarse as possible, yet still resolving the presence of fine-scale features, e.g. in a heterogeneous microstructure there need to be grid points in the neighbourhood of the microstructural features to pick up their influence on the initial grid, even if this is still very inaccurate.

Integration of the wavelet representation
The internal force and tangent stiffness matrix required for the Newton-Raphson procedure result from integrating the wavelet representations of the internal force inte- The MRA approximation of an integrand approximationf (x) is defined using a set of coarse scaling coefficients and a sparse set of wavelet coefficients, respectively. The fieldf (x) is then given by: To evaluate the integral of the weak form the unit-integral property of the Deslauriers-Dubuc interpolating wavelets is used [37]. The integral of the mother scaling function φ(θ) constructed on a dyadic wavelet grid with sufficient grid points (≥ 2m − 1), coordinate θ and a grid spacing Δθ = 1 is given by [4].
Using the biorthogonal refinement relations [18] the mother wavelet function is expressed as a function of the scaling function.
For Deslauriers-Dubuc interpolating wavelets this simplifies to the following expression.
Mapping the coordinate θ onto the physical coordinate x, the integrals of the Deslauriers-Dubuc scaling functions and wavelets on an equidistant dyadic grid, with Δx j the grid spacing on level j, are given by: Substitution yields the following simple summation of the coefficients times the grid size at the coefficient level.
Combining the wavelet approximation of the integrands ϕ and κ and the integration and substitution into a Newton-Raphson procedure to solve the linear momentum balance leads to the W-ROM outlined in Algorithm 1. Note that the integration grid is rebuilt within every iteration. When the history coefficients from the previous increment or iteration are required, the wavelet representation thereof is stored. Upon refinement, the history variables need to be interpolated between the available grid points when the history is required on newly added grid points. The Deslauriers-Dubuc basis functions used to discretise the integrands are also used to interpolate the history variables.

Require:
The external force f ext , the macroscopic strain ε M , the Newton-Raphson and MRA tolerances δ NR and δ w , the sampled strain modes ∂ x R(x j i ) and the wavelet representation of the history variables s 0 Compute stress, history ξ t+1 and tangent stiffness Refine intermediate grid points Use wavelet synthesis to retrieve interpolated values σ Compute stress, history ξ t+1 and tangent stiffness ϕ Store the new history parameters

Numerical examples
To evaluate the performance of the W-ROM, a onedimensional microstructure with elasto-plastic material behaviour is modelled. The required integration points are monitored to assess the computational costs. The accuracy of the resulting macroscopic forces f M (ε M ) of the full order model (FOM), the ROM and the W-ROM are compared.
The accuracy with respect to the standard ROM is evaluated by comparing the micro-fluctuation coefficients W . To demonstrate the flexibility of the wavelet integration a second microstructure consisting of two domains with different material parameters is modelled. The spatially dependent cross-sectional area given by Eq. (21) is used for both models. An average cross-sectional area of A avg = 0.1 mm 2 and an area fluctuation of ΔA = 0.02 mm 2 are used for both examples. A wavelength of L = 0.5 mm is selected for a microstructure with length = 1 mm. The cross-sectional area is then given by: To homogenise the one-dimensional problem, the following relation between the macroscopic strain ε M and the macroscopic force f M is introduced.
The full-order model is discretised with 1000 quadratic Lagrangian elements and the weak form is integrated using 2 Gauss-Legendre quadrature points per element. The dyadic wavelet grid has 9 points at level j = 0 such that the 4th degree Deslauriers-Dubuc scaling functions and wavelets are fully supported within the domain. Note however that a finer initial grid may be used when this is required. Periodic boundary conditions are employed for both the finite element and the wavelet discretisations.

Example I: Micro-structural model with one material
The first example consists of a unit-cell, schematically depicted in Fig. 6, with an elasto-plastic material having a Young's modulus of E = 1.0 GPa, a yield-stress of σ y = 0.02 GPa, and a hardening coefficient of K = 0.4 GPa.
The model is loaded with a macroscopic strain ε M increasing from 0 to 0.1 in 10 equidistant increments. Snapshots of the resulting micro-fluctuation fields w(x) are stored in the snapshot-matrix X to extract the modes constituting the reduced basis R(x). Using the POD, the modes and their corresponding eigenvalues required to formulate the reduced basis R(x) are retrieved. The modes and singular values are plotted in Fig. 7. In this example, 3 modes are used to capture the micro-fluctuation field. After the derivation of a reduced basis, it is projected onto the dyadic wavelet grid to complete the W-ROM.
The accuracy of the reduced models is assessed by comparing the reduced micro-fluctuation coefficients and the resulting integrated macroscopic force obtained using the ROM and W-ROM (using Algorithm 1) with the results obtained using the FOM. The relation between computational cost and accuracy of the W-ROM is investigated by adopting different tolerances and comparing the corresponding sparse wavelet grids.

Accuracy
To assess the accuracy, the macroscopic strain ε M is increased from 0 to 0.1 in 20 equal increments. Note that the macroscopic strain is applied in twice the number of increments, i.e. every odd strain increment lies in the middle of two strains used to generate the snapshots X such that the ROM needs to interpolate between the available modes. The resulting micro-scale stress and strain fields obtained with the FOM, Micro-scale forces (inc. 20) Fig. 8 Microscopic stress, strain and force fields for the FOM, ROM and W-ROM (δ w = 10 −2 ). A section of the force curve is enlarged 5× to visualise the small deviations in the force equilibrium ROM, and W-ROM are plotted in Fig. 8 for load increments 5, 7 and 20 (using a tolerance δ w = 10 −2 ).
Due to the limited number of snapshots and modes used, the ROM is not able to represent the force equilibrium exactly close to the plastic deformation front. Small deviations of the force equilibrium are visible between the ROM and FOM results. Since the W-ROM relies on the same basis, these small deviations are also present here. In order to obtain a more accurate representation of the strain field close to the plastic deformation front, more snapshots are required to obtain the missing modes to construct the ROM. Moreover, if the ROM problem would be solved with a wavelet approximation as well, this concern would be completely alleviated. This will be explored in future work.
Note, however, that the W-ROM approximates the reduced order model accurately using only a reduced set of integration points. The strains of the W-ROM are compared to the lossless integrated result of the ROM. Since both models make use of the same basis (in the wavelet representation the strains are sampled from the reduced basis directly), it suffices to compare the reduced coefficients. The evolution of the coefficients over the load increments is plotted in Fig. 9, where the integration error in the coefficients of the W-ROM model is defined with respect to the ROM model as follows.
To investigate the influence of the wavelet representation on the macroscopic force f M , the macroscopic force-strain curve for the FOM, ROM and W-ROM is shown in Fig. 10. The error f M in the approximation of the macroscopic force f M using the ROM and W-ROM model with respect to the FOM model is defined as where f * M is the macroscopic force approximated with either the fully integrated ROM or W-ROM model.
Even though the tolerance δ w = 10 −4 only applies to the approximation of the internal force, the error in the wavelet representation of the macroscopic force remains of approximately the same order. The ROM captures the macroscopic force f M exactly up to numerical precision in the linear regime and in the post-yielding regime up to the given tolerance. In this regime, the wavelet representation proves sufficiently accurate to maintain the ROM accuracy. The influence of the applied tolerance on the macroscopic force f M is investigated by approximating the internal forces in the W-ROM model using different tolerances δ w = 10 −2 , 10 −3 , . . . , 10 −6 . The resulting errors of the ROM and W-ROM with respect to the FOM case are plotted in Fig. 11.
A stricter tolerance consistently leads to a better force approximation. The errors for the tolerances δ w = 10 −5 and δ w = 10 −6 are plotted in dashed lines since the wavelet representation reached the maximum level (here set to j max = 12) and accordingly the reduction in error stagnates. A larger j max would allow to capture more details, through which the approximation of the internal force integrand ϕ(x 0 ) becomes even more accurate resulting in a more accurate force integral.
This becomes clear by considering the sparse wavelet grids for the various tolerances δ w , as depicted in Fig. 12. At tolerance δ w = 10 −4 the maximum level is reached and the grid starts to fill in. The accuracy of the approximation can be enhanced using a higher maximum level j max , or by increasing the order of the bases. Sparse grid δ w = 10 −5 Fig. 12 Sparse dyadic grids used for the wavelet representation using various tolerances δ w

Reduction
The sparsity of the wavelet representation allows the approximation of the micro-scale stress and internal force integrand fields using a limited number of sample points. Depending on the tolerance, the algorithm automatically identifies points at locations where refinement is needed, thereby eliminating the need for an a priori determined set of integration points [20] or modal reconstruction of the micro-scale stress field [21]. The wavelet basis enables an approximation up until a preset tolerance picking up local phenomena in high detail while using few sample points in regions without significant fluctuations. When the maximum level is reached, the wavelet coefficient gives an indication of the order of magnitude of the remaining approximation error. The number of points required for an approximation with tolerance δ w is depicted in Fig. 13. For a tolerance of δ w = 10 −2 in the internal force integrand only 76 integration points are required, compared to the original FEM problem with 2000 integration points. The method gains on the full order model up to a tolerance δ w = 10 −5 , which is remarkable considering there are only limited regions without strong fluctuations in the force integrand field in the one-dimensional example model.
Note, that for a tolerance δ w ≤ 10 −5 the wavelet representation of the integrand requires more integration points than the original Gauss scheme to project the integrand onto the wavelet basis. The accuracy of the integrand in the FOM how- When using a tolerance δ w lower than the accuracy of the FOM the wavelet representation will approximate the artefacts of the original Lagrangian discretisation of the finite element problem. This locally requires extra sampling points which do not yield extra accuracy of the physical solution, but merely in the discretised solution. If a higher accuracy is desired a finer finite element mesh discretisation for the FOM needs to be employed first.
Evidently, this emphasises again that the tolerance used in the wavelet reduction/approximation should be balanced with respect to the underlying discretisation error of the FOM. Note that this problem will no longer exist if wavelets Micro-scale forces (inc. 20) Fig. 15 The evolution of the micro-scale strain (left-hand side), stress (middle) force (right-hand side) approximated with a FOM, ROM and W-ROM model for increments 5, 7 and 20. A section of the force curve is enlarged × 5 to visualise the small deviations in the force equilibrium are used to solve the BVP for the FOM/ROM as will be done in future work.

Example II: Two-phase microstructural composite rod
Next, a somewhat more complex microstructure is considered. Similar to the first case, this problem consists of an expanding plastically deformed zone. However, this microstructure consists of two phases of elasto-plastic material which yield at different stress-levels with a different hardening rate, resulting in a discontinuous micro-scale strain field. The microstructure is schematically depicted in Fig. 14.
The resulting micro-scale stresses and strains are plotted in Fig. 15 for the load increments 5, 7 and 20. At increment 5 the onset of plastic deformation is visible on the right hand side of the microstructure, at increment 7 the right phase shows a significant plastically deformed zone and at the final increment 20 both phases reveal plasticity. When both phases are deforming plastically, a strain discontinuity originates due to the difference in yield stress and hardening parameters. Here one can see a small approximation error in the ROM and the W-ROM stress, strain and force field. This error is however localised due to the (automatically) refined wavelet grid near the discontinuity and hardly influences the integrated macroscopic stress. As stated before, this problem can be solved by adopting a wavelet family which allows for discontinuities [39]. Figure 16 shows the resulting micro-fluctuations resolved with ROM and W-ROM. Comparing the W-ROM method relative to the fully integrated ROM method, the errors in the reduced micro-fluctuation parameters W are of the order O 10 −4 , which is much smaller than the imposed tolerance of δ w = 10 −3 .
When looking at the macroscopic force-strain curve in Fig. 17 of the FOM, ROM and W-ROM models, a similar trend is observed. In the plastic regime, the error in the ROM and W-ROM with respect to the FOM are of order O 10 −4 using approximately 300 integration points. This is considerably smaller than the imposed tolerance.
The wavelet grid used for the W-ROM approximation presented in Fig. 18 shows that the discontinuities at the boundaries and centre of the domain are adequately picked up by the automatic refinement strategy until the maximum wavelet level j max is used. The W-ROM model uses ∼ 300 Due to the C 0 kinks originating at the plastic front, the refinement algorithm locally uses a lot of sampling points. In this work it is shown that the local high accuracy can be obtained. If a high reduction factor is favoured, the user has several options to fine tune the compression: (i) choosing the maximum wavelet level lower; (ii) a less strict tolerance; (iii) using higher order bases. This implies that (i) the minimum width of the wavelets increases or (ii) the approximation does not refine that much in narrow regions around kinks and discontinuities; (iii) fine scale details are approximated using less sampling points. Options (i) and (ii) will give a locally less accurate microscale solution. Note however, that the error defined locally is quite strict when one is interested in the homogenised stresses or forces, since the error made in the homogenised result is often lower than the local error.

Conclusions
This paper presents a novel hyper-reduced method by introducing a wavelet basis and a MRA to integrate the weak form up to a pre-set tolerance. The innovative aspects of the proposed method are: -The wavelet reduced integration does not require an offline calibration step using a second set of (stressbased) snapshots to construct the reduced integration scheme. This reduces the input parameters required to construct the integration scheme to the strain modes only (projected on the dyadic wavelet grid) and a tolerance to control the level of approximation. -The reduced integration uses the wavelet-based MRA to enable an adaptive and local refinement of the dyadic grid used to integrate the stress, internal forces and stiffness. Conversely, the adaptive dyadic grid will also coarsen when smooth functions are to be approximated. -The error in the integration is controlled using a pre-set tolerance. This is verified by quantifying the microfluctuation error of the W-ROM relative to the FOM. The resulting error in micro-fluctuation coefficients is typically bounded by the imposed tolerance.
The reduction of the number of stress evaluations to perform the integration with respect to the original Gauss integration scheme used in the ROM is demonstrated and the compression is shown to be inversely proportional to the tolerance.
Hyper-reduced models such as Empirical Interpolation Methods and Empirical Cubature Methods rely on offline determined integration points, and sometimes a basis for the integrand. This allows for a higher compression ratio, but requires a detailed sampling of the high-dimensional snapshot space including the history parameters for path dependent materials, to ensure that no physical modes are missing since the integration scheme is constructed a priori. Many hyper-reduction techniques combine the sampling of the non-linear terms, such as the stress field σ or the internal force integrand ϕ, with the sampling of the kinematics. Therefore no additional computations are required for the construction of the snapshots of the nonlinear term. The accuracy of the reduced integration of history dependent material models however, requires extensive sampling of both the kinetics and the kinematics present in the parameter space.
The Wavelet-Reduced Order Model, on the other hand, adaptively determines the points required for accurate integration in-situ and only requires prior accurate sampling of the kinematics. The sampled kinematics are required to construct the symmetric gradient of the reduced basis ∇ 2 R(x) and is therefore less sensitive to the sampled snapshots, and algorithmic choices to select integration points.
Note that there is no theoretical limitation to expand W-ROM to two or three dimensional microstructural models. This 1D model provides a transparent view of all underlying principles and implementation aspects of the method proposed. Expanding the method to multiple dimensions will however require multi-dimensional adaptive wavelet transforms as presented by Paolucci et al. in [32] and more elaborate storage solutions for the internal variables and modes by approximating each variable independently on a separate wavelet grid to allow fast wavelet transforms while still limiting the memory usage.
To approximate the integrands, a sequence of approximation spaces is used, V j containing the scaling functions φ j i (x). The scaling function spaces in the sequence are subspaces of their successors, i.e. V j ⊂ V j+1 . The wavelet spaces W j complement the scaling function spaces V j such that the combination of both spaces V j ⊕ W j are spanning the (more detailed) space V j+1 . This hierarchical relation between the wavelets and scaling functions on different levels is schematically depicted in Fig. 20.
In this work, the wavelets are translated along grid points on level j spaced Δx j apart. At each level the wavelet is scaled down using a factor 1 2 in width for each subsequent level forming a dyadic wavelet grid. The basis consisting of coarse scaling functions and increasingly fine wavelet functions living on multiple levels is formed. In each subsequent level j, the wavelet functions ψ To perform a MRA on a function f (x), without considering the computational efficiency, the function can be projected onto the fine level scaling functions φ j max i (x). The scaling function representation can be decomposed into a combination of scaling functions and wavelets at the lower level j max − 1. The projection of the function f (x) on the scaling functions φ j max −1 i (x) can be analysed again using the wavelet analysis to obtain a representation in terms of level j = j max − 2 scaling functions φ In the multi-resolution analysis the wavelet transform is applied hierarchically to decompose the different frequency components in the function f (x). The multi-resolution analysis is schematically shown for the Haar wavelet family in Fig. 21.
On a dyadic wavelet grid the wavelet analysis and synthesis can be performed discretely using a set of filter coefficients hierarchically relating the scaling function coefficients and wavelet coefficients. The level j scaling function and wavelet coefficients are related to the level j +1 scaling function coefficients via the analysis filter coefficients denoted byh k and g k . The synthesis of the level j scaling function and wavelet coefficients to the level j + 1 scaling function coefficients is performed using the synthesis filter coefficients h k and g k . For a more detailed explanation of multi-resolution analysis the reader is referred to the books [6,26].

A.1.1 Wavelet synthesis
Starting from a coarse representation of the function f (x) comprising of 0th level scaling functions φ 0 i (x), an increas- Table 1 The wavelet coefficients used to generate the 4th degree interpolating Deslauriers-Dubuc wavelet and scaling functions ingly detailed reconstruction of the function f (x) is obtained by adding higher level wavelets. This reconstruction relies on the wavelet synthesis relations to interpolate intermediate points on finer grid levels j + 1. The synthesis of the discretised field is performed using the synthesis relations for bi-orthogonal wavelets (25). Here, m is the degree of the current filter, and h k and g k are the synthesis filter coefficients listed in Table 1 used to perform the synthesis. The analysis filter coefficients are denoted bỹ h k andg k . The 4th degree Deslauriers-Dubuc level j = 0 scaling functions φ 0 i (x) and wavelets ψ 0 i (x) are depicted in Fig. 22. When using these Deslauriers-Dubuc interpolating wavelets, the wavelet coefficients d j+1 i on the fine level are given directly by the difference between the level j wavelet interpolation and the sampled function value f (x j i ) 1 . This 1 Note that this property is specific for the Deslauriers-Dubuc interpolating wavelet family.
is shown for wavelet synthesis on a dyadic grid using the 4th degree (m = 4) bi-orthogonal Deslauriers-Dubuc scaling function coefficients h k and wavelet coefficients g k . A schematic overview of the synthesis is depicted in Fig. 23.
The interpolated (or approximated) values on the intermediate grid points can be found by interpolating wavelets assuming that the wavelet coefficients d For Deslauriers-Dubuc wavelets the filter coefficient g k is described using the Kronecker-delta function −δ 1 k . Therefore the only value d j i−k contributing to approximationf j+1 2i+1 is k = 0. This coefficient represents the error between the interpolated function valuef Note that relation (27) is specific for the Deslauriers-Dubuc interpolating wavelet. For other wavelets the general synthe-sis relations for bi-orthogonal wavelets (25) can be applied to find the wavelet coefficients.

A.2 Data compression
Using the Deslauriers-Dubuc wavelet family allows the wavelet representation of the function f (x) to be constructed from the coarse level up. By hierarchically applying equation (27) over each level, the correcting wavelet coefficients d j i are computed. Only the points for which the error |d j i |/| f avg | ≥ δ w , in which δ w is the tolerance on the approximation error relative to the absolute average of the function values on the initial grid | f avg |, are marked for further refinement. When a point x j i is marked for further refinement the surrounding grid points on level j + 1, x j+1 2i and x j+1 2i+1 are approximated in the next level, providing a systematic refinement scheme. The MRA approximationf (x/ ) is proven to converge to the function f (x/ ) as the number of levels tends to infinity [6].
The process of truncating wavelets is schematically depicted for the Deslauriers-Dubuc wavelet family in Fig. 24. In this example several wavelets have coefficients d j i that are below the imposed tolerance at level j = 3. These wavelets are no longer considered for further refinement in higher lev- els. At level j = 4 all coefficients are below the imposed tolerance δ w . The wavelet representation is sufficiently accurate and all contributions of higher level wavelets can be neglected.