Stochastic virtual element methods for uncertainty propagation of stochastic linear elasticity

This paper presents stochastic virtual element methods for propagating uncertainty in linear elastic stochastic problems. We first derive stochastic virtual element equations for 2D and 3D linear elastic problems that may involve uncertainties in material properties, external forces, etc. A stochastic virtual element space that couples the deterministic virtual element space and the stochastic space is constructed for this purpose and used to approximate the unknown stochastic solution. Two numerical frameworks are then developed to solve the derived stochastic virtual element equations, including a Polynomial Chaos approximation based approach and a weakly intrusive approximation based approach. In the PC based framework, the stochastic solution is approximated using the Polynomial Chaos basis and solved via an augmented deterministic virtual element equation that is generated by applying the stochastic Galerkin procedure to the original stochastic virtual element equation. In the weakly intrusive approximation based framework, the stochastic solution is approximated by a summation of a set of products of random variables and deterministic vectors, where the deterministic vectors are solved via converting the original stochastic problem to deterministic virtual element equations by the stochastic Galerkin approach, and the random variables are solved via converting the original stochastic problem to one-dimensional stochastic algebraic equations by the classical Galerkin procedure. This method avoids the curse of dimensionality of high-dimensional stochastic problems successfully since all random inputs are embedded into one-dimensional stochastic algebraic equations whose computational effort weakly depends on the stochastic dimension. Numerical results on 2D and 3D problems with low- and high-dimensional random inputs demonstrate the good performance of the proposed methods.


Introduction
Numerical techniques for solving complex partial differential equations are continuously developing at an incredible rate, e.g.finite difference methods, finite volume methods, finite element methods and spectral element methods etc [1].As a generalization of the finite element method, the virtual element method (VEM) has been proposed and received a lot of attention in the last decade [2,3,4,5].Compared to the classical finite element method, VEM can discretize 2D or 3D geometric domains utilizing arbitrary polygons or polyhedrons and is not limited to the regular elements used in the finite element method, which is thus highly flexible and mesh insensitive.In VEM, the shape functions can be non-polynomial.It does not require constructing explicit shape functions on elements since all numerical integration is transferred to edges rather than performed on the elements.Several limitations of finite element methods are also avoided, such as convex elements and element degradation caused by small edges and interior angles.Applications of VEM to various problems have been extensively studied, such as linear elastic problems [6,7,8], large deformation problems [9,10], contact problems [11,12,13], fracture and crack propagation problems [14,15,16], topology optimization [17,18] etc.Although extensive studies have been conducted, there is still a gap in using VEM to deal with problems with uncertainties.In many practical engineering problems, the inherent or epistemic uncertainty of systems are unavoidable.
Predicting uncertainty propagation on the physical models has become an important part of the analysis of systems [19], which leads to the development of dedicated numerical methods for uncertainty analysis.
In this paper, we focus on extending the deterministic VEM to stochastic VEMs (SVEMs) for the uncertainty analysis of 2D and 3D linear elastic stochastic problems that may involve random material properties and stochastic external forces, etc., which is currently still lacking.Our main contributions in this paper consist of two parts: the first contribution is to extend the deterministic virtual element discretization to stochastic cases and derive corresponding stochastic virtual ele-ment equations (SVEEs), and the second contribution is to present two numerical methods to solve the derived SVEEs efficiently and accurately.For the first contribution, we extend the deterministic virtual element space to a stochastic virtual element space that couples the classical virtual element space and the stochastic space, which can provide (stochastic) virtual element approximations for both deterministic and stochastic functions in the space.The constructed stochastic virtual element space can be considered as a deterministic virtual element space parameterized by random inputs.For each sample realization of the random input, it degenerates into a deterministic virtual element space and inherits all properties of classical virtual element spaces.Thus, we can simply approximate the stochastic solution using a linear combination of deterministic virtual basis functions with random coefficients (i.e. the unknown stochastic solution).Further, numerical techniques for calculating the gradients of virtual basis functions and the stabilization term in the deterministic VEM [3,7,8,18] can still be used to the stochastic discretization with slight modifications.In this way we can obtain SVEEs by assembling the stochastic stiffness matrix and the stochastic force vector with a complexity similar to the deterministic VEM.
Similar to VEM being a generalization of the finite element method, SVEM is also considered to be a generalization of the stochastic finite element method [20].Thus, numerical algorithms for solving the derived SVEEs can benefit from stochastic finite element solution algorithms, such as Monte Carlo simulation (MCS) and its improvements [21,22], spectral stochastic finite element methods [23,24], stochastic collocation methods [25,26], response surface and kriging methods [27,28], etc.For the second contribution, we develop two numerical methods to solve the derived SVEE, including a Polynomial Chaos expansion based SVEM (PC-SVEM) and a Weakly INtrusive approximation based SVEM (WIN-SVEM).The other methods mentioned above can also be extended to SVEMs in a similar way as in this paper.The PC-SVEM is a natural extension of the spectral stochastic finite element method [23].In this method, the stochastic solution is decomposed into a summation of a set of products of PC basis and deterministic vectors.By the use of stochastic Galerkin procedure, the original SVEE is transformed into an augmented deterministic equation whose size is much larger than the original SVEE.Also, the size increases dramatically as the degree of freedom of physical models, the stochastic dimension and the expansion order of PC basis increase, which leads to the curse of dimensionality when dealing with large-scale and/or high-dimensional stochastic problems.To address this issue, we further present a WIN-SVEM, which is an extension of our previous work for solving stochastic finite element equations [29,30].In this method, the stochastic solution is approximated by a summation of a set of products of random variables and deterministic vectors.Different from the PC-SVEM, both random variables and deterministic vectors are unknown a priori.To this end, we solve them using a dedicated iteration.The deterministic vectors are solved via a few number of deterministic equations that are obtained by a similar stochastic Galerkin process used to PC-SVEM.The random variables are solved via one-dimensional stochastic algebraic equations that are obtained by applying the classical Galerkin procedure to the original SVEE.In this way, all random inputs are embedded into these one-dimensional stochastic algebraic equations, and their efficient solutions are achieved using a non-intrusive sampling method with weak dimension dependence.
The proposed WIN-SVEM thus avoids the curse of dimensionality of high-dimensional stochastic problems successfully.
The paper is organized as follows: Section 2 presents the stochastic virtual element discretization for linear elastic stochastic problems and stochastic virtual element equations are then derived.In Section 3, the PC-SVEM is developed to solve the derived stochastic systems.Following that, the WIN-SVEM is proposed in Section 4 to solve the derived stochastic systems efficiently, with special emphasis on high-dimensional stochastic problems.2D and 3D numerical examples involving low-and high-dimensional random inputs are given in Section 5 to demonstrate the performance of the proposed methods.Conclusions and outlook follow in Section 6.

Stochastic elastic equations
Let (Θ, Ξ, P) be a suitable probability space , where Θ denotes the space of elementary events, Ξ is a σ-algebra defined on Θ and P is a probability measure.In this paper, we consider the following elastic stochastic equation where the deterministic domain Ω ⊂ R d with the boundary ∂Ω, the spatial dimension may be d = 2, 3 and the d-dimensional spatial coordinate is given by x denotes the divergence operator, σ (x, θ) is the stochastic stress tensor, the vector-valued dis- In this paper, we only consider linear elastic stochastic problems.The following linear stochastic strain tensor and linear elastic constitutive relation are adopted where C (x, θ) is a fourth order elastic tensor that may be related to stochastic material properties, e,g.stochastic Young's modulus and stochastic Poisson ratio.
To solve Eq. ( 1), let us consider its weak form written as follows: find a stochastic function where the functional space is defined as and H 1 (Ω) is the subspace of the space of square integrable scalar functions on Ω that contains both the function and its weak derivatives [31].The terms W (u (x, θ) , v (x) , θ) and F (v (x) , θ) are given by

Stochastic virtual element discretization
We adopt the stochastic virtual element discretization for the weak form Eq. (3).Specifically, the domain Ω is partitioned into n e non-overlapping polygonal elements Ω = n e e=1 Ω (e) , and each element Ω (e) , e = 1, . . ., n e includes n (e) vertices and m (e) edges.In this paper, we only consider the low-order virtual element, but the proposed method can be extended to higher-order virtual elements [32].We give the following approximate discretized virtual space V h Ω (e) ⊂ V of the element Ω (e) where P 1 represents the space of polynomials of degree up to 1, v h is a polynomial on each edge of Ω (e) and C 0 -continuity on the element Ω (e) , and ∇ • [C (θ) ε (v h )] vanishes on the element Ω (e) for all θ ∈ Θ.In this way, we couple the classical virtual element space and the random input θ.The function v h is not known on the element Ω (e) but explicitly known on the edge ∂Ω (e) .We now consider the weak form Eq. (3) on the discretized space: find u h (x, θ) ∈ V h : Ω × Θ → R d such that the following equation holds for P-almost surely θ ∈ Θ, where W h and F h are assembled by looping through all elements Ω (e) n e e=1 It is noted that W (e) h (u h (x, θ) , v h (x) , θ) and F (e) h (v h (x) , θ) cannot be evaluated in a similar way to the classical finite element method since the function v h is unknown on the element Ω (e) .To this end, a projection operator Π (e) : V h Ω (e) → P 1 Ω (e) d is defined similar to the deterministic VEM such that for ∀v h ∈ V h Ω (e) , p 1 ∈ P 1 Ω (e) d , θ ∈ Θ, holds.We can thus reformulate the term W (e) h (u h (x, θ) , v h (x) , θ) as where W (e) C (θ) and W (e) S (θ) are called the stochastic consistent term and the stochastic stabilization term, respectively.The last two terms are obtained according to Eq. ( 10) since Π (e) v h (x), . Furthermore, ε (e) Π (e) u h (x, θ) and ε (e) Π (e) v h (x) are (stochastic) constant strain tensors.Therefore, the stochastic consistent term W (e) C (θ) is evaluated via C (e) (x, θ) ε (e) Π (e) u h (x, θ) : ε (e) Π (e) v h (x) dx = a (e) C (e) (x, θ) ε (e) Π (e) u h (x, θ) : ε (e) Π (e) v h (x) , where a (e) is the area (for 2D polygonal element) or the volume (for 3D polygonal element) of the element Ω (e) .
Further, we let the stochastic solution u h (x, θ) and the function v h (x) on the element Ω (e) be approximated using a set of virtual basis functions where the solution vector u (e) i (θ) of the element T ∈ R d is the solution vector of the i-th vertex.The vector v (e) i ∈ R n (e) d has a similar expression but does not involve the random input θ.Applying the projection operator Π (e) to u h (x, θ) and v h (x) we have On the basis of this we rewrite the strain tensors ε (e) Π (e) u h (x, θ) and ε (e) Π (e) v h (x) as the following vector forms where the matrices e) are given by where the components ∂Π (e) ϕ i (x) calculations to edges and faces of the element Ω (e) based on Eq. ( 10), which is the same as the deterministic virtual element method, see [2,3,7,8] for details.In this way, it does not require knowing explicit representations of the functions {ϕ i (x)} n (e) i=1 and only needs to know their traces on edges.The Lagrangian linear basis functions similar to those used in the classical finite element method can be adopted for the purpose.Substituting Eq. ( 16) into Eq.( 13) we have where is the matrix form of the tensor C (x, θ) of the element Ω (e) .Two detailed representations of G (θ) for 2D and 3D problems can be found in the numerical example section.Hence, the stochastic element stiffness matrix corresponding to the stochastic consistent term W (e) C (θ) is given by Further, let us consider the calculation of the stochastic stabilization term W (e) S (θ) in Eq. ( 12), which can be achieved taking advantage of several numerical strategies [2,7,8].Here we adopt the approach presented in [3,18], which corresponds to where denotes the i-th coordinate value of the j-th vertex of the element Ω (e) .In the practical implementation, it is calculated as [18] W where the coefficient γ (e) S (θ) = 1 2 Tr G (e) (θ) , Tr (•) is the trace operator of matrices, I n (e) d ∈ R n (e) d×n (e) d is the identity matrix, and the deterministic matrix S (e) ∈ R n (e) d×n (e) d is given by where the vectors X (e) i and A (e) i are given by Hence, the stochastic element stiffness matrix corresponding to the stochastic stabilization term which is very close to that in the deterministic VEM, but the coefficient γ (e) S (θ) involves the random input θ.The calculation procedure of deterministic virtual element matrices can thus be inherited for k (e) S (θ).We only need to pay a little attention to the calculation of the random coefficient γ (e) S (θ).
In the last step, let us consider the stochastic term F (e) h (v h (x) , θ) calculated exactly using the one-point integration rule on the edges and the face (for 2D case) or the faces and the element (for 3D case) where n (e) Γ N is the number of edges (for 2D case) or faces (for 3D case) of the element Ω (e) , Γ (e) N, j is the j-th edge or face and includes n (e) Γ N, j vertices, and b (e) j is the length (for 2D element) or the area (for 3D element) of Γ (e) N, j .The vectors consisting of the values on each vertex are N, j N, j , where N, j and 0 otherwise.Thus, the stochastic element force vector is calculated as n (e) Assembling the above stochastic element stiffness matrices and stochastic element force vector we obtain the following SVEE where the global stochastic stiffness matrix K (θ) and the the global stochastic force vector F (θ) are assembled by where (•) represents the assembly operator for all stochastic element matrices and vectors, and the total degree of freedom (DoF) is given by n = n e d.
Further, the stochastic matrix G (x, θ) (or the tensor C (x, θ)) can be approximated using the following separated form in many cases where are deterministic matrices.For the non-separated stochastic matrix G (x, θ), the approaches for simulating random fields can be adopted to achieve Eq. ( 29)-like approximations for both Gaussian and non-Gaussian random inputs, e.g.Karhunen-Loève expansion and Polynomial Chaos expansion [33,34,35].In this way, we can reformulate the SVEE (27) as where the deterministic matrices {K i } m i=0 are assembled via since the stochastic element matrices k (e) C (θ) in Eq. ( 19) and k (e) S (θ) in Eq. ( 24) depend linearly on the matrix components {G i } m i=0 .

PC-SVEM: Polynomial Chaos based spectral stochastic virtual element method
The PC-based methods have been well developed and widely applied to solve a variety of stochastic problems.In this section, we present a PC-SVEM to solve the SVEE (27) (or the separated form Eq. ( 30)).In this method, the stochastic solution u (θ) is expanded using (generalized) PC basis as follows where d PC,i k i=1 are the corresponding deterministic vectors to be solved, {Γ i (θ)} k i=1 are the PC basis.In practice, we can choose different PC basis for different types of random inputs, such as the Hermite PC basis for Gaussian random variables and the Legendre PC basis for uniform random variables [23,24].The stochastic Galerkin approach is then used to transform SVEE (27) into the following deterministic equation [23] where P (θ) is the probability measurement of the random input θ.Further, the above equation can be rewritten as a compact form where the augmented deterministic matrix K PC ∈ R nk×nk and the augmented deterministic vectors d PC , F PC ∈ R nk are assembled by where the matrix and vector components K PC, ji ∈ R n×n and Further, if the separated form Eq. ( 30) is considered, the above calculation of the matrices K PC, ji is simplified as which only involves the numerical integration Θ ξ l (θ) Γ i (θ) Γ j (θ) dP (θ).For low-dimensional stochastic problems, the calculation is cheap enough benefiting from efficient numerical integration strategies on stochastic spaces [26].
It is noted that the total number of PC basis is k = (m+r)!m!r! , where (•)! is the factorial operator, r is the expansion order of the PC basis.Thus, similar to classical PC-based methods, the proposed PC-SVEM still suffers from the curse of dimensionality since the matrix/vector size nk in Eq. (34) increases sharply as the spatial DoF n of the physical model, the stochastic dimension m of the random input and the expansion order r of the PC basis increase.For instance, the size is about nk = 1 × 10 6 when n = 1 × 10 3 , m = 10 and r = 4, which requires extremely expensive computational effort.Although several methods are developed to alleviate the computational burden, e.g.dedicated iterative algorithms and sparse PC approximations [36,37,38].It is still challenging to solve very high-dimensional stochastic problems using the PC-SVEM.

A weakly intrusive stochastic virtual element method
To avoid the curse of dimensionality arising in the above PC-SVEM, we present a WIN-SVEM in this section, which can be considered as an extension of our previous work [29,30] on stochastic finite element methods to SVEM.To this end, we consider the stochastic solution u (θ) approximated by the following series expansion where {λ i (θ) ∈ R} k i=1 are scalar random variables, d WIN,i ∈ R n k i=1 are deterministic vectors, and k is the number of retained terms.It is noted that the number k and all pairs λ i (θ) , d WIN,i k i=1 are not known a priori.An iterative algorithm is presented to solve the pairs λ i (θ) , d WIN,i one by one.Specifically, we assume that the (k − 1)-th approximation u WIN,k−1 (θ) = k−1 i=1 λ i (θ) d WIN,i has been known and the goal is to solve the k-th pair λ k (θ) , d WIN,k .The original SVEE ( 27) can be rewritten as where the stochastic vector In this way, Eq. ( 39) only involves one unknown pair λ k (θ) , d WIN,k .However, different from that the random basis (i.e.PC basis) has been known in PC-SVEM, both the random variable λ k (θ) and the deterministic vector d WIN,k are unknown in this case.To avoid this issue, an alternating iteration is adopted to solve them.Specifically, if the random variable λ k (θ) has been known (or given an initial value), Eq. ( 39) is transformed into the following deterministic virtual element equation by taking advantage of the stochastic Galerkin procedure [23] similar to that in Eq. ( 33) which is equivalent to where the deterministic matrix K WIN,k = Θ K (θ) λ 2 k (θ) dP (θ) ∈ R n×n and the deterministic vector Existing numerical solvers can be adopted to solve it effi-ciently and accurately [39].Note that the size of Eq. ( 41) is the same as the original SVEE (27), which is different from the augmented size of PC-based derived equation (34) and can thus save a lot of computational effort.In practical implementations, we let the vector d WIN,k orthogonal to the obtained vectors d WIN,i k−1 i=1 to speed up the convergence, which is achieved by using the Gram-Schmidt orthonormalization where d WIN,i k−1 i=1 are normalized orthogonal vectors that meet d T WIN,i d WIN, j = δ i j , where δ i j is the Kronecker delta.
With the deterministic vector d WIN,k solved using Eq. ( 41), the random variable λ k (θ) is then recalculated taking advantage of the following classical Galerkin procedure Since the stochastic matrix K (θ) is positive definite and z T K (θ) z > 0, ∀z 0 ∈ R n , ∀θ ∈ Θ holds, the above equation can be rewritten as To avoid the curse of dimensionality arising in the high-dimensional problems, we adopt a nonintrusive sampling approach [30] to solve Eq. (44) instead of the PC-based approximation, which corresponds to where θ represents n s sample realizations of (θ), , respectively, and represents the element-wise division of two sample vectors, also known as Hadamard division operator.In this way, all random inputs are embedded into the sample realization vectors , which is insensitive to the stochastic dimension of random inputs.The curse of dimensionality can thus be avoided successfully, which will be discussed in detail in the next section.
We can calculate the k-th pair λ k (θ) , d WIN,k by performing the iterative process of Eq. (41) and Eq. ( 45) until reaching a specified precision.A similar iteration is also adopted to calculate other pairs λ k+1 (θ) , d WIN,k+1 , • • • until a good approximation of the stochastic solution is achieved.However, it is noted that the stochastic solution u WIN,k (θ) in Eq. ( 38) is approximated in a sequential way and it does not exactly fulfill the original SVEE (27).The approximation has low accuracy for some cases [29].We introduce a recalculation process to avoid this problem.
To this end, considered as a set of reduced basis functions and the random variable vector which requires to be solved repeatedly for n s sample realizations to get the final solution Λ θ (i) , i = 1, • • • , n s , but only very low computational effort is involved since the sizes of the reducedorder stochastic matrix D T WIN K (θ) D WIN ∈ R k×k and the reduced-order stochastic vector D T WIN F (θ) ∈ R k are greatly reduced compared to the original SVEE (27) in most cases.
Let us highlight the weak intrusiveness of the proposed method.On one hand, Eq. ( 38) is considered as a kind of intrusive approximation of the stochastic solution that is very similar to the PC-based intrusive approximation (32).However, on the other hand, the implementation for solving d WIN,k in Eq. ( 41) only involves deterministic calculations and the matrix K WIN,k keeps the same size and matrix properties as the original stochastic matrix K (θ), which is weakly intrusive.Also, Eq. (45) for calculating the random variable λ k (θ) is fully non-intrusive.In these senses, we implement the intrusive stochastic solution approximation only in weakly intrusive and fully non-intrusive ways.The method combines the high efficiency of intrusive methods and the weak dimensionality dependence of non-intrusive methods.It can solve high-dimensional stochastic problems efficiently and accurately.

High-dimensional stochastic problems
In this section, we will show that the proposed WIN-SVEM can be applied to high-dimensional stochastic problems without any modification.We explain this point from the perspective of the influence of high stochastic dimensions on solving Eq. ( 41) and Eq.(45).We only consider Eq. (30) in this section and a large number m is truncated in Eq. ( 30) to generate a high-dimensional stochastic problem.In this way, the deterministic matrix K WIN,k and the deterministic vector F WIN,k in Eq. ( 41) are calculated via where the probability integrals are approximated using the following non-intrusive sampling approach where are sample vectors of the random variables {ξ i (θ)} m i=0 , the operator represents the element-wise multiplication of sample vectors, E {•} is the expectation operator of the sample vector.Eq. (49) takes a total of k (m + 1) expectation operations, which is not sensitive to the stochastic dimension and has low computational effort even for very high stochastic dimensions.Note that although we only illustrate the high-dimensional input in the stochastic matrix K (θ), the above calculation also works efficiently for Θ F (θ) λ k (θ) dP (θ) if the stochastic vector F (θ) involves high-dimensional random inputs.
Further, the sample vectors in right side of Eq. ( 45) are calculated via which requires a total of k (m + 1) operations for , k and is also computationally cheap for high-dimensional stochastic problems.Therefore, both Eq.(41) and Eq. ( 45) are insensitive to the stochastic dimension.The proposed method can avoid the curse of dimensionality successfully.

Algorithm implementation
The above proposed WIN-SVEM for solving SVEEs is summarized in Algorithm 1, which includes two loop processes.The inner loop is from step 4 to step 9 and used to solve the kth pair λ k (θ) , d WIN,k .The outer loop from step 2 to step 13 is to approximate the stochastic solution using a set of pairs λ i (θ) , d WIN,i k i=1 .To execute the inner loop, a random sample vector λ (0) k θ ∈ R n s is initialized in step 3.In the numerical implementation, any nonzero vectors of size n s can be chosen as the initialization since the initial samples have little influence on the computational accuracy and efficiency of the proposed method.Following each inner loop, we only need to update the stochastic force vector F k+1 (θ) in step 10 and store the reduced-order matrix D WIN in step 11 in the outer loop.
Algorithm 1 WIN-SVEM for solving SVEEs  Orthonormalize d ( j) WIN,k using Eq. ( 42) 7: Update the random sample vector λ ( j) k θ ∈ R n s via Eq.(45) 8: Compute the locally iterative error d, j , j ← j + 1 9: Update the stochastic force vector Update the deterministic matrix Compute the locally iterative error u,k , k ← k + 1 13: end 14: Recalculate the random variable vector Λ (θ) ∈ R k via Eq.( 46) Two iterative criteria are involved in the above algorithm to check the convergence, i.e. d, j in step 8 for the inner loop and u,k in step 12 for the outer loop.The iterative error d, j is defined as which measures the difference between the vectors d ( j) WIN,k and d ( j−1) WIN,k and the calculation is stopped when d, j < d is met.Similarly, the iterative error u,k is defined as which measures the contribution of the k-th pair λ k (θ) , d WIN,k to the stochastic solution u WIN,k (θ).
However, Eq. ( 53) may be not a good error indicator in many cases [29] since the random variables {λ i (θ)} k i=1 are calculated in a sequential way and Θ λ 2 k (θ) dP (θ) may not keep decreasing.We avoid this problem by replacing {λ i (θ)} k i=1 in Eq. ( 53) with equivalent random variables . To this end, we calculate the autocorrelation function of the random variable vector which is decomposed by eigendecomposition into where Q ∈ R k×k is an orthonormal matrix and Z is a diagonal matrix consisting of descending eigenvalues of the matrix C ΛΛ .We construct an equivalent random variable vector Substituting the equivalent random variables λ i (θ) into Eq.( 53) we recalculate the iterative error u,k as where Z k is the k-th diagonal element of the matrix Z.It is noted that the above reformulation does not improve the approximation accuracy of the stochastic solution and only provides an equivalent representation.In this way, the iterative error u,k keeps decreasing as the retained item k increases.
More details of the basic implementation and the comparison between Eq. ( 53) and Eq. ( 56) can be found in [29].

Numerical examples
In this section, we test the proposed two methods with the aid of 2D and 3D numerical examples.For Algorithm 1, the convergence errors d in step 2 for the inner loop and u in step 4 for the outer loop are set as 1 × 10 −3 and 1 × 10 −6 , respectively.For both examples, n s = 1 × 10 4 random samples are used for performing MC simulations and generating reference solutions.In our cases, 1×10 4 samples are enough to achieve convergent probabilistic solutions.The same 1×10 4 random samples are also used in the proposed WIN-SFEM to eliminate the influence caused by sampling processes.Further, the examples are performed on one core of a desktop computer (sixteen cores, Intel Core i7, 2.50GHz).

Simulation of random inputs
In this example, we consider the SVEM-based plane stress analysis of a 2D model shown in Fig.
1 [40].The Voronoi mesh is adopted for the spatial discretization, including a total of 2018 vertices, 1000 Voronoi elements and 4036 DoFs.The model is fixed at the two red points as shown in Fig. 1.A stochastic force f (θ) = −1000 − 100ξ f (θ) (unit: N) is applied to the blue point along the y direction, where ξ f (θ) is a standard Gaussian random variable.Further, the material property matrix G (x, y, θ) is given by where the Poisson ratio ν = 0.3 and the Young's modulus E (x, y, θ) is modeled as a two-dimensional random field with the mean value E 0 (x, y) = 100 MPa and the covariance function [41] Cov EE (x 1 , y 1 ; where the standard deviation σ E = 10 MPa, and l x and l y are the correlation lengths in the x and y directions, respectively.The random field E (x, y, θ) is approximated via the following Karhunen-Loève expansion [23,34] where {ξ i (θ)} m i=1 are mutually independent standard Gaussian random variables and they are also independent of the random variable ξ f (θ).{κ i , E i (x, y)} m i=1 are eigenvalues and eigenvectors of the covariance function Cov EE (x 1 , y 1 ; x 2 , y 2 ).They are solved by the following Fredholm integral equation of the second kind which can be solved efficiently by taking advantage of existing eigenvalue solvers [42].To ensure E (x, y, θ) > 0, the samples θ (i) such that min x,y∈Ω E x, y, θ (i) ≤ 1 × 10 −3 are dropped out in numerical implementations.In this way, the deterministic matrices {G i (x, y)} m i=0 in Eq. ( 29) are given by In this example, we consider a low-dimensional case and a high-dimensional case via introducing different truncated items m in Eq. (59).For the low-dimensional case, we let the correlation lengths be l x = max(x) − min(x) and l y = max(y) − min(y).To achieve the truncated error κ m m i=1 κ i < 1 × 10 −3 , the number of truncated items is m = 6.Corresponding eigenvectors {E i (x, y)} 6 i=1 of the covariance function Cov EE (x 1 , y 1 ; x 2 , y 2 ) are shown in Fig. 2. It is noted that we can solve Eq. (60) using the discretized covariance matrix Cov EE ∈ R 2018×2018 of Cov EE (x 1 , y 1 ; x 2 , y 2 ), which is only dependent of the vertices and independent of elements.In this way, the eigenvectors in Fig. 2 are plotted on each vertex of the mesh.For the high-dimensional case, the correlation lengths are set as l x = 1 8 [max(x) − min(x)] and l y = 1 8 max(y) − min(y) , and m = 34 truncated items are retained to achieve κ 34 high-dimensional case, it is slower to converge to the specified truncated error and more truncated items are required to capture the local correlation property.

Low-dimensional case
In this section, we solve the low-dimensional case using the proposed PC-SVEM and WIN-SVEM.For the PC-SVEM, the second order Hermite PC basis of seven standard Gaussian random variables ({ξ i (θ)} 6 i=1 and ξ f (θ)) is adopted for {Γ i (θ)} 36 i=1 .The size of the augmented deterministic equation ( 34) is 145296, which is much larger than the 4036 DoFs of the original stochastic problem.We do not perform the numerical implementations for higher order PC basis in this example.
On one hand, the second order PC basis is enough to achieve a good approximation of the stochastic solution.On the other hand, the size of the augmented equation ( 34) is 484320 if the third order PC basis is adopted, which leads to too high computational burden in terms of storage and solution for a problem of such a small spatial scale.To show the computational accuracy of the proposed PC-SVEM and WIN-SVEM, we compare PDFs of the stochastic displacements u A,x (θ) and u A,y (θ) in the x and y directions of the point A (i.e. the blue point where the force f (θ) acts as shown in Fig. 1) obtained by PC-, WIN-and MCSbased SVEMs and their absolute errors in Fig. 6.It is seen from Fig. 6a and Fig. 6c that the PDFs of both stochastic displacements u A,x (θ) and u A,y (θ) obtained by PC-and WIN-SVEMs are in good accordance with those of MCS, which illustrates the good accuracy of the two proposed methods.
The comparison in logarithmic scales shown in Fig. 6b and Fig. 6d demonstrates that WIN-SVEM can achieve smaller absolute errors than PC-SVEM, especially for the stochastic displacement u A,y (θ) in the y direction.Further, WIN-SVEM can capture tails of the PDFs more accurately than PC-SVEM, which is very useful for many uncertainty quantification problems, such as the simulation of physical phenomena with long tailed probability distributions and the estimation of small failure probability in structural reliability analysis.Therefore, WIN-SVEM is recommended for problems with such requirements.Further, let us focus on the computational efficiency of the proposed methods.Computational times (unit: second) of the numerical execution of PC-SVEM, WIN-SVEM and MCS are listed in the second to fourth columns of Table 1, where the solving and recalculating times of WIN-SVEM are the computational times of step 2 to step 13 and the recalculation step 14 of Algorithm In this section, we solve the high-dimensional case with a total of 35 stochastic dimensions ({ξ i (θ)} 34 i=1 and ξ f (θ)).Only WIN-SVEM is adopted in this section.For PC-SVEM, even if we only use the second order PC basis, the size of the derived deterministic equation ( 34) is about 2.69 × 10 6 , which suffers from the curse of dimensionality.By using Algorithm 1, the iterative errors u,k of different retained terms for the high-dimensional case are shown in Fig. 7.The proposed WIN-SVEM still has good convergence for high-dimensional stochastic problems.k = 6 terms are retained in this case, which is slightly increased compared to the low-dimensional case.Computational times for this case are seen from the fifth and sixth columns of Table 1.WIN-SVEM has very low cost even for high-dimensional cases and is much cheaper than MCS.Compared to the low-dimensional case, the computational cost of the high-dimensional case does not increase dramatically as the stochastic dimension increases.In these senses, the proposed WIN-SVEM avoids the curse of dimensionality successfully.Further, PDFs of the stochastic displacements u A,x (θ) and u A,y (θ) in the x and y directions of the point A obtained by WIN-SVEM and MCS are compared in Fig. 8.For both two stochastic displacements, the computational accuracy of WIN-SVEM is still comparable to MCS.In this case, we consider the SVEM analysis of a 3D mechanical part as shown in Fig. 9a, where Dirichlet boundary conditions u x (θ) = u y (θ) = u z (θ) = 0 are imposed on the red surface and an external force f (x, y, z) = −500 N/mm 2 is applied along the x direction on the blue surface.
The model is discretized by use of the Voronoi mesh depicted in Fig. 9b, including a total of 28232 vertices, 4389 elements and 84696 DoFs.In this example, the stochastic material matrix G (x, y, z, θ) is given by where the Poisson ratio ν = 0.3 and the Young's modulus E (x, y, z, θ) is a three-dimensional random field with the covariance function where the standard deviation σ E = 41.8GPa, and the correlation lengths are given by l x = max(x)− min(x), l y = max(y)−min(y) and l z = max(z)−min(z).The random field E (x, y, z, θ) has a Eq. ( 59)like series expansion where the function E 0 (x, y, z) = 208 GPa, {ξ i (θ)} m i=1 are mutually independent uniform random variables on [0, 1].It is noted that the mean value of the random field E (x, y, z, θ) is E 0 (x, y, z) + 0.5 m i=1 √ κ i E i (x, y, z) instead of E 0 (x, y, z).Similarly, {κ i , E i (x, y)} m i=1 are eigenvalues and eigenvectors of the covariance function Cov EE (x 1 , y 1 , z 1 ; x 2 , y 2 , z 2 ) and can be solved by the Eq.(60)-like integral equation.The truncated number is set as m = 13 in this case to achieve the truncated error    stochastic problems.Regarding the computational accuracy, PDFs of the stochastic displacements u A,x (θ), u A,y (θ) and u A,z (θ) in the x, y and z directions of the point A (shown in Fig. 9b) obtained by WIN-SVEM and MCS are compared in Fig. 11.For the PDFs of stochastic displacements u A,x (θ) and u A,z (θ), WIN-SVEM is in very good agreement with MCS.The PDF of the stochastic displacement u A,y (θ) is slightly less accurate than those of u A,x (θ) and u A,z (θ), but acceptable accuracy is still achieved.We can simply retain more terms in the stochastic solution approximation to improve the accuracy if higher accuracy is required in practice.Further, the computational times of WIN-SVEM and MCS in this example are 213.45s(including the solving time 212.89s and the recalculating time 0.56s) and 2.74 × 10 4 s, respectively.The proposed WIN-SVEM is still much cheaper than MCS and a speedup of more than one hundred times is achieved.We also highlight that statistical properties of the stochastic solution are easily computed based on WIN-SVEM.Here we focus on the first and second order global statistical moments, that is, the mean value and the standard deviation dependent on spatial positions.According to Eq. ( 38), the mean value vector u is computed as which only involves the expectations of sample vectors λ i θ ∈ R 10 k i=1 and has very low computational effort.The mean value components u x , u y and u z of u in the x, y and z directions are seen from the first line of Fig. 12. Further, the standard deviation vector σ u is calculated via whose components σ u x , σ u y and σ u z in the x, y and z directions can be found in the second line of Fig. 12.

Conclusions
We presented two numerical approaches, PC-SVEM and WIN-SVEM, for solving stochastic systems derived from the stochastic virtual element discretization of 2D and 3D linear elastic stochastic problems.The deterministic virtual element method is first extended to SVEM that involves stochastic material properties and stochastic external forces, etc.Several key calculations of SVEM can be inherited from the deterministic virtual element method, e.g. the gradient computations of virtual basis functions in Eq. ( 17) and the stabilizing element stiffness matrix in Eq. (24) (except for the random coefficient γ (e) S (θ)).Numerical results demonstrate that both PC-SVEM and WIN-SVEM have comparable accuracy to MCS.However, PC-SVEM suffers from the curse of dimensionality and cannot be applied to high-dimensional stochastic problems.As a comparison, WIN-SVEM can efficiently solve both low-and high-dimensional stochastic problems without any modification, which has been verified by a numerical example of up to 35 stochastic dimensions.
Further, although only linear elastic stochastic problems are concerned in this paper, both the proposed PC-SVEM and WIN-SVEM can be applied to more general cases, which will be further investigated in subsequent research.
associated with stochastic external forces, and Γ N and Γ D are boundary segments associated with the Neumann boundary condition g s

Figure 1 :
Figure 1: Geometry and Voronoi mesh of the 2D model.

Figure 4 :
Figure 4: Iterative errors of different retained terms in the low-dimensional case.

Figure 5 :
Figure 5: Components of the stochastic solution: the deterministic vectors d u x ,i 4 i=1 in the x direction (the first line), the deterministic vectors d u y ,i 4 i=1 in the y direction (the second line) and PDFs of the random variables {λ i (θ)} 4 i=1 (the third line), respectively.
Comparison of PDFs of the stochastic displacement u A,x (θ) in the x direction of the point A obtained by PC-, WIN-SVEMs and MCS, respectively.
Absolute errors of PDFs of the stochastic displacement u A,x (θ) in the x direction of the point A obtained by PC-and WIN-SVEMs to the MCS reference PDF.
Comparison of PDFs of the stochastic displacement u A,y (θ) in the y direction of the point A obtained by PC-, WIN-SVEMs and MCS, respectively.
Absolute errors of PDFs of the stochastic displacement u A,y (θ) in the y direction of the point A obtained by PC-and WIN-SVEMs to the MCS reference PDF.

Figure 6 :
Figure 6: PDFs of the stochastic displacements u A,x (θ) and u A,y (θ) in the x and y directions of the point A obtained by PC-, WIN-SVEMs and MCS and their absolute errors.

Figure 7 :
Figure 7: Iterative errors of different retained terms in the high-dimensional case.
Comparison of PDFs of the stochastic displacement u A,x (θ) in the x direction of the point A obtained by WIN-SVEM and MCS, respectively.
Comparison of PDFs of the stochastic displacement u A,y (θ) in the y direction of the point A obtained by WIN-SVEM and MCS, respectively.

Figure 8 :
Figure 8: PDFs of the stochastic displacements u A,x (θ) (left) and u A,y (θ) (right) in the x and y directions of the point A obtained by WIN-SVEM and MCS, respectively.

Figure 9 :
Figure 9: Geometry of the 3D mechanical part (left) and its Voronoi mesh (right).
Comparison of PDFs of the stochastic displacement u A,y (θ) in the y direction of the point A obtained by WIN-SVEM and MCS, respectively.
Comparison of PDFs of the stochastic displacement u A,z (θ) in the z direction of the point A obtained by WIN-SVEM and MCS, respectively.

Figure 11 :
Figure 11: PDFs of the stochastic displacements u A,x (θ) (top left), u A,y (θ) (top right) and u A,z (θ) (bottom) in the x, y and z directions of the point A obtained by WIN-SVEM and MCS, respectively.

Figure 12 :
Figure 12: The mean functions u x , u y and u z (the first line) in the x, y and z directions, and the standard derivation functions σ u x , σ u y and σ u z (the second line) in the x, y and z directions.

Table 1 :
1, respectively.The cost of the recalculation process of WIN-SVEM is almost negligible since only deterministic equations with size 4 are solved for different sample realizations.It is found that both PC-SVEM and WIN-SVEM are much cheaper than MCS.However, since a very large deterministic equation needs to be solved in PC-SVEM, it is more computationally intensive than WIN-SVEM.Computational costs of the stochastic dimensions 7 and 35.