A multiscale virtual element method for the analysis of heterogeneous media

We introduce a novel heterogeneous multiscale method for the elastic analysis of two‐dimensional domains with a complex microstructure. To this end, the multiscale finite element method is revisited and originally upgraded by introducing virtual element discretizations at the microscale, hence allowing for generalized polygonal and nonconvex elements. The microscale is upscaled through the numerical evaluation of a set of multiscale basis functions. The solution of the equilibrium equations is performed at the coarse scale at a reduced computational cost. We discuss the computation of the multiscale basis functions and corresponding virtual projection operators. The performance of the method in terms of accuracy and computational efficiency is evaluated through a set of numerical examples.


INTRODUCTION
Advances in automated manufacturing and, in particular, additive manufacturing have led to the widespread application of components possessing complex and fit-for-purpose material layouts in the construction, aerospace, and automotive industries. 1 Additively manufactured functionally graded composites and foams can be tailored to increased mechanical properties when compared with traditional layered composites or metals, for example, higher strength to weight ratios and higher damping to weight ratios. 2 However, the corresponding manufacturing processes can be extensive, are prone to errors, and necessitate several design iterations before a desirable layout is finally produced. This motivates the development of computational methods that can lead to augmenting desirable mechanical traits while reducing undesirable ones while still at the design stage.
Yet, the flexibility provided by manufacturing poses a series of challenges vis-a-vis the numerical simulation of structural components with exotic material layouts. The distribution of material heterogeneities at the microscale or mesoscale significantly affects the overall mechanical response of the component. Hence, both from a physical and a computational perspective, the problem of assessing the mechanical performance of a component becomes a multiscale one. Numerical simulation of physics across multiple length scales can, in principle, be done with the standard finite element Sreekumar and Triantafyllou are equally Contributing. approach. 3 However, this would necessitate the use of extremely fine mesh discretizations to resolve the corresponding heterogeneities, hence leading to high-computational costs. 4 Computational homogenization techniques have been developed to efficiently address this issue (see . These are based on the definition of a representative volume element (RVE) 9 that forms the basis for evaluating the effective constitutive behavior of a heterogeneous domain. 10,11 The structural problem is then solved at a macroscopic scale using standard discretization methods, for example, finite elements. However, computational homogenization methods rely on the assumptions of periodicity and scale separation, which are often not the case in highly heterogeneous domains. Increasing the size of the RVE 12,13 to accommodate more information about the heterogeneities will lead to increased computational costs, effectively negating the advantages of upscaling.
Multiscale finite element methods (MsFEM) 14 offer a very interesting alternative to account for such heterogeneous problems. MsFEM use the information about the morphology at the microscopic or mesoscopic scale to construct equivalent numerical problems subject to certain kinematical constraints. Solutions to these equivalent problems yield multiscale basis functions, which are used to map mesoscale information to the macroscopic scale. These basis functions are dependent on the geometry and material properties of the underlying constituents. An inherent drawback of this method, that prevents its application to structural engineering problems, is its inability to account for the Poisson effect. The enhanced multiscale finite element method (EMsFEM) 15 has been developed to overcome this problem by introducing additional coupling terms when computing the multiscale basis functions. However, the EMsFEM can only treat quadrilateral meshes.
Accounting for complex morphologies of the heterogeneities encompassed by an RVE using traditional element geometries such as quadrilateral and triangular elements may necessitate quite fine mesoscopic mesh discretizations, thus driving up the cost of computing multiscale basis functions. Optimization of the underlying mesh could thus prove critical to improving the performance of the method. Mesh optimization requires, as a prerequisite, numerical methods that allow for more flexible mesh generation capabilities.
Polygonal and polyhedral finite element methods (PFEM) [16][17][18][19][20][21][22] are increasingly being applied to several fields of computational mechanics and are well suited in modeling complex microstructure morphologies. To this point, these formulations have been applied to topology optimization, 23,24 computational fracture and damage modeling, [25][26][27][28] contact problems, 29,30 and fluid mechanics. 31 However, PFEM are limited by the availability of suitable basis functions. The typically employed analytical Wachspress shape functions are valid only over simple and convex polyhedra. 32 Hence, for the case of complex and, in general, nonconvex polygons and polyhedrons, numerical approaches have to be employed that, however, increase the computational toll of the method 33 as in the case of, for example, numerically evaluated harmonic 17,34 and maximum-entropy shape functions. [35][36][37] This issue is further augmented in nonlinear problems. 38,39 The virtual element method (VEM) emerged as a versatile alternative for solving partial differential equations using flexible element geometries to address the aforementioned issues. [40][41][42][43][44][45][46] The development of the VEM can chronologically be traced back to finite volume methods 47,48 and more recently the mimetic finite difference method (MFD). [49][50][51][52][53][54] More specifically, the VEM emerged as a variational reformulation of the low-order MFD method 51 and its higher order generalization. 55 MFDs, when extended to finite element methods (FEMs), were able to model traditional trial/test function spaces over generalized element geometries without having explicit representations for these basis functions inside the element interior. The aforementioned approaches aim to improve the accuracy of the standard FEM, when applied to polyhedral meshes, by enriching the space of trial functions with possibly nonpolynomial expressions. Contrary to the MFD, the VEM attempts to preserve the polynomial accuracy over simplexes. 40 This allows the use of complex polyhedral and nonconvex element shapes with more general continuity requirements such as (div) conformity. 56 In this work, we derive a novel heterogeneous multiscale method, termed the enhanced multiscale virtual element method (EMsVEM), for the analysis of highly heterogeneous domains across multiple length scales. To achieve this, we employ the VEM to resolve heterogeneities at the fine scale and derive appropriate multiscale basis functions to project the VEM onto a coarse finite element mesh. The solution of the governing equations is then performed at the coarse mesh at a significantly reduced computational cost. We proceed to examine and discuss the advantages and drawbacks that arise from employing virtual elements at the fine scale in terms of accuracy and computational efficiency.
The remainder of the article is organized as follows: in Section 2, the formulation of the VEM for linear elasticity is discussed. The derivations pertaining to the EMsVEM are provided in Section 3. The method is validated against the standard FEM and analytical solutions in Section 4 and the merits and bottlenecks of the EMsVEM are discussed. Concluding remarks are provided in Section 5.

Problem statement
The case of a domain Ω ⊂ R 2 is considered herein subjected to plane stress or plane strain conditions. The domain is subjected to a body force b and is prescribed a displacement and traction u and t, respectively over the boundaries Γ u and Γ t , respectively. The displacement and traction boundaries are defined to be nonintersecting, that is, Γ u ∩ Γ t = Ø. The corresponding equilibrium, constitutive, and compatibility equations are defined in Equations (1a), (1b), and (1c), respectively, as subjected to generalized inhomogeneous displacement and traction boundaries defined as follows: Here, C is the material constitutive tensor, is the stress tensor, and is the infinitesimal strain tensor. The displacement field u is being solved for. The symbols ∇, ∇ S , and div(⋅) denote the gradient, symmetric gradient, and divergence operators, respectively. Equation (1) is recast into their corresponding weak form, considering the following function spaces where [ 1 (Ω)] 2 represents the standard Hilbert function space of dimension 2. Defining an appropriate weighting function v ∈  0 , the weak formulation of the governing equations defined in Equation (1) is expressed in the following compact form where a(⋅, ⋅) and F(⋅) are bilinear and linear functionals, respectively, defined as In this work, the virtual element approach is employed to discretize Equation (5).

Virtual element discretization
The aim of the VEM is to generalize the standard finite element domain decomposition to polygonal (or polyhedral in three-dimensional) elements with any number of edges (not restricted to merely quadrilateral, triangular, hexahedral, or tetrahedral meshes). This includes also the case of nonconvex elements. Let Ω be decomposed as shown in Figure 1B into n el nonoverlapping two-dimensional polygonal elements contained in a set  h * . The corner vertices of each element where N i v denotes the number of corner vertices. The edges are denoted by e i j for j = 1, 2, … N i e , where N i e = N i v is the total number of edges. Each edge e i j connects vertices i j and i j+1 . For a kth-order VEM, the edge e i j is assumed to contain k − 1 internal edge-points. The case of a second-order polygonal element is shown in Figure 2. In the following,  i will be written as  for brevity.
Similar to the standard FEM, the approximations u h and v h to the actual displacement and weighting field, respectively, are defined as where  h and  h0 are termed "virtual spaces." *The parameter h is interpreted as the maximum diameter of all elements contained in  h .
Hence, Equation (4) can be established in the following discretized form where the bilinear forms a(⋅, ⋅) and F(⋅) in Equation (8) are evaluated from the piecewise element bilinear operators a  (⋅, ⋅) and F  (⋅), respectively, as The global virtual spaces  h and  h0 and local element virtual spaces   h and   h0 are part of standard virtual element literature 57 (see explicit definitions provided in Appendix A for the sake of completeness). Remark 1. The spaces employed within the VEM are examples of generalized finite element spaces that allow the trial and test fields to assume nonpolynomial forms over an element. This is seen as a weaker construction of the more conventional finite element space, where the field variables are restricted to only members of a polynomial space.
The displacement approximation, that is, u h is implicitly interpolated at the virtual element degrees of freedom (DoFs) that include The values of the DOFs are denoted by the operator dof(⋅) that assumes the following expression where [M k−2 ()] 2 denotes a space of k − 2-order vector monomials of dimension 2. The weighting function v h is also interpolated at the same DoFs.
Explicit definitions for these vector monomial spaces are provided in Appendix B for the sake of completeness. The DoFs of the space have been illustrated over an arbitrary polygonal element for k = 1, 2, 3 in Figure 3. Hence, for the case of a two-dimensional domain where the displacements assume two DoFs per node, the dimension of the virtual element space is Since the VEM functions do not have an explicit expression over the domain, essential operations such as ∇(⋅) and consequently, (⋅) (Equation (1c)) cannot directly be performed. Conversely, a strain-specific projection operator Π k ∶   h () → [P k ()] 2 is defined according to the following optimality criterion where [P k ()] 2 is a space of kth-order two-dimensional vector polynomials. This criterion ensures that the error arising from the polynomial approximation is energetically orthogonal to the approximating subspace [P k ()] 2 . It follows that the strain energy is computed exactly, despite the simplifying assumption introduced by the projection. In VEM literature, this property is called k-consistency.
Remark 2. Introducing the strain specific Π k projection onto a polynomial space is a deviation from the classical Galerkin framework. Approximating a nonpolynomial function using a polynomial basis falls into a class of "variational crimes" 40 and introduces additional error into the formulation. This error is a result of performing numerical integration for nonpolynomial functions using polynomial quadrature rules. However, one can derive a priori bounds and estimators to quantify and control this error (see Reference 58).
The approximating subspace [P k ()] 2 is spanned by vector monomials belonging to [M k ()] 2 . However, these monomials also contain zero deformation modes, that is, rigid body modes that contribute zero strain energy to the formulation described in Equation (12). To avoid spurious results arising from ill-conditioned or singular matrices, such modes are excluded when computing Π k . Instead, they are treated in a stability term. This is discussed in Section 3.2.
Following this reasoning, Equation (12) is eventually established in the following form where K () belongs to the kernel of rigid body motions and [M k ()] 2 ∖K () is the kernel of the strain operator (⋅).
The contents of these sets can be determined using the additive kinematic decomposition rules mentioned in Reference 59. The monomial spaces used for the VEM formulation are provided in Appendix B for the sake of completeness. The exact procedure followed for computing Π k and the associated local stiffness matrix a  (u h , v h ) are discussed within a multiscale context in Section 3.2.

Overview
To this point, the EMsFEM has been developed to treat regular heterogeneous domains as shown in Figure 4A, where each rectangular coarse element clusters its own representative portion of the underlying fine mesh. The fine mesh is designed to resolve all mesoscale heterogeneities. These are upscaled to the coarse level, where the solution is obtained at reduced computational cost. In this work, we treat the most general case of complex geometrical domains at the microstructure as shown in Figure 4B. To achieve this, we employ the VEM to accurately and efficiently resolve the heterogeneities at the fine scale. Similar to the EMsFEM, the heterogeneous domain is coarsely discretized into  M( ) , = 1 … n M el , coarse elements, where n M el is the number of coarse elements. Each coarse element clusters its own set of  m(i) , i = 1 … n m el microelements, where n m el is the number of microelements in the  M( ) .
For each type of coarse element, a set of multiscale basis functions is evaluated using the standard VEM. These basis functions are then employed to map the fine scale onto the coarse scale, where the solution of the governing equations is performed. The EMsVEM procedure is schematically depicted in Figure 5.
The multiscale basis functions required for the upscaling procedure are evaluated through the solution of a homogeneous version of Equation (8) over the th coarse element domain  M( ) , that is, These homogeneous equations are subjected to kinematical constraints that account for heterogeneities and the deformability of the boundary. They are imposed over the RVE in the form of linear or periodic boundaries. The procedure for such enforcements are discussed in Reference 2 and will not be detailed here. The choice of the coarse element boundary conditions plays an important role vis-a-vis the accuracy of the method. Resonance errors may occur when the coarse-element length scale approaches the fine-element scale, that is, a large number of coarse elements are employed to mesh the domain. Such errors can be overcome by an oversampling strategy. 60

Virtual fine-scale stiffness matrices
The bilinear operator used in the coarse element subproblem, defined in Equation (14), is decomposed in a manner analogous to Equation (9): Contrary to the standard EMsFEM, in this work, we employ virtual elements at the fine scale. Consequently, we adopt a VEM discretization for the elementwise bilinear functional shown in Equation (15): where Π m k denotes the projector discussed in Section 2.2.
Remark 3. To account for the potential nonpolynomial expressions, standard Lagrange polynomials are not a suitable basis. Rather, a canonical basis is utilized that is implicitly approximated in terms of a linear combination of conventional scaled vector monomials. This is convenient because such monomials have simple analytic expressions. The problem of giving an explicit form to the canonical basis is now reduced to computing the coefficients involved in the linear combination. To accomplish this, functions belonging to virtual spaces need to be projected onto a polynomial space Expanding Equation (16) and exploiting the symmetry of the bilinear functional and the energetic orthogonality condition in Equation (13), the following relation is obtained where the first term consists entirely of polynomial expressions and can be exactly computed analytically (k = 1 methods) or numerically (k ≥ 2 methods) through the DOFs defined in Equation (10). This is called the consistency term and involves the elliptic projection of the virtual element functions, that is, Π m k u h and Π m k v h , respectively. The second term on the right-hand side of the expression contains unknown nonpolynomial integrands u h and v h . These terms are not computable as elements of a virtual element space lack an explicit definition over the element interior.
The noncomputable term is replaced with a computable bilinear form   m (⋅, ⋅) that can account for strain energy associated with higher order deformation modes while ensuring coercivity; this is called the stability term. Hence, the final VEM approximation of the strain energy is expressed as There exists several ways to define the approximated stability operator   m (⋅, ⋅). One is referred to References 61 and 62 for a more comprehensive discussion on the stability term. In this work, we adopt the D-recipe stabilization 63 that also ensures correct scaling in relation to the consistency term; the corresponding definition is provided in Appendix E for the sake of completeness.
Remark 4. The consistency term on its own is rank deficient. The stability term is included to restore coercivity. In the context of quadrilateral or hexahedral elements, the stability term can be interpreted as the hourglass stiffness. The term can assume any bilinear form that satisfies the fundamental properties of coercivity and stability and reduces to zero over polynomial subspaces. The easy-to-compute stability term introduces additional error into the formulation. However, this approximation is chosen such that the error convergence rates are still optimal. It is for this reason that one should not expect the VEM formulation to reduce to the familiar FEM approach in the case of, for example, quadrilateral elements. In the select case of the three-noded triangular element, the polynomial space used is complete to degree one. In this case, the stability term reduces to zero over polynomial subspaces and the consistency term coincides with the classical three-noded FEM stiffness matrix.
The fine scale stiffness matrix K el, m(i) is now defined on the basis of Equation (18): where, K C is the consistency stiffness matrix and K S is the stability stiffness matrix whose expressions are provided in Appendices D and E, respectively. As discussed in Section 2.2, the strains (⋅) cannot be directly evaluated on u h and v h as these are implicitly defined. The projector Π m k provides the functions with an explicit form by approximating them with a monomial expansion. The projector Π m k is evaluated using the following orthogonality condition at the microscale The aforementioned monomial expansion for the term Π m k u h is expressed as follows: where The quantity n k is and denotes the number of vector-valued monomials in Substituting Equation (21) into Equation (20) and simplifying, the following expression is derived The procedure to evaluate the fine scale matricesG and B is provided in Appendix C. Equation (23) is conveniently recast in matrix form as where m k is a vector comprising the coefficients s i . Hence, the fine-scale projection operator is eventually derived as The microelement stiffness matrix in Equation (19) is used to assemble the coarse element-specific RVE stiffness matrix K m using a standard direct stiffness approach, that is, where A is the assembly operator. The coarse element-specific RVE load vector f m is also similarly assembled from its microelement contributions as where f el, m(i) is the nodal force vector defined at the fine scale. This is established on the basis of Equations (5b) and (9b) as where Γ t denotes the traction boundaries associated with the element under consideration. The body force and traction terms can be computed by defining appropriate nodal quadrature rules over the microelement interior and edge, respectively, as detailed in Reference 64.

Constructing multiscale basis functions
In the multiscale finite element framework utilized herein, the microdisplacement components of the fine mesh are mapped to the macrodisplacement nodal components of the corresponding coarse element according to Equation (29).
where u mx,i and u my,i , i = 1, … , n m are the displacement components of the ith micronode, n m is the number of micronodes within the coarse element, and n M = 4 is the number of macronodes of the coarse element. Moreover, u Mx,J and u My,J are the macrodisplacement components defined at the Jth macronode of the coarse element; N iJxx , N iJyy , N iJxy , and N iJyx correspond to the microbasis interpolation functions. Equations (29) hold if and only if the microbasis interpolation functions satisfy the following conditions At the fine element scale, this mapping assumes the following form where u m(i) is the displacement vector for the ith fine element in the th RVE. The array N m(i) denotes the multiscale basis functions mapping the associated coarse-nodal displacements in u M( ) to the current fine element. Collecting Equation (31) for all microelements within the coarse element, the following equation is established where u m comprises the nodal displacement of the fine mesh contained within coarse element . A set of interpolation functions satisfying Equation (30) (hence Equation (29) also) can be established by solving the following boundary value problem where K m is the RVE stiffness matrix that is assembled according to Equation (26), u S is the vector of boundary micro-nodal displacements, and u are the the prescribed displacements obtained from imposing linear, periodic or oversampling boundary conditions over the coarse-element/RVE boundary (see References 2,60). Each solution derived for a prescribed set of boundary conditions u IJ corresponds to a column of the multiscale basis function matrix N m . Algorithm 1 summarizes the procedural steps required for the evaluation of the microbasis functions.
Remark 5. To retain consistency with the literature, each coarse element is termed herein an RVE. This, however, should not be confused with the typical definition of the RVE as employed within a scale separated homogenization theory.
As an example, Equation (31) assumes the following form for the case of the RVE shown in Figure 6 and the microelement i = 6, where the vector of nodal displacements for fine element #6 is and the corresponding multiscale basis function matrix is expressed as 1 N xy7,1 N xx8,1 N xy8,1 N xx13,1 N xy13,1 N xx12,1 N xy12,1  N xy7,1 N yy7,1 N xy8,1 N yy8,1 N xy13,1 N yy13,1 N xy12,1 N yy12 respectively. For the RVE shown in Figure 6, u m is a 34 × 1 vector. Furthermore, N m is the coarse-element matrix of multiscale basis functions; for the case of Figure 6 this is a 34 × 8 matrix. Each column j of N m corresponds to the deformed configuration of the RVE when the jth coarse DoF is equal to 1 whereas all other coarse DOFs are equal to zero.

Governing multiscale equilibrium equations
The elementwise equilibrium equations at the fine scale of the RVE under consideration are where K el, m(i) and f el, m(i) denote the fine-element stiffness matrix and load vector defined in Equations (19) and (28), respectively. Substituting Equation (31) into Equation (38) results in the following expression: Premultiplying this by N T m(i) , the following relation is obtained: where K el M( ),m(i) corresponds to fine-element stiffness matrix, that is however, defined at the coarse nodes and f el M( ),m(i) to the corresponding vector of nodal forces respectively. The reduced order coarse element equilibrium equation can be established in the following form where K el M( ) and f el M( ) are the coarse element stiffness matrix and load vector, respectively. These quantities are a priori unknown but can be derived on the basis of strain energy equivalence between the coarse element and its underlying fine scale mesh.
The overall strain energy E int of the coarse element can be, in principle, established as where and correspond to the strain and stress fields defined over the RVE. However, the RVE strain energy can also be considered to be additively decomposed into the contributions of the associated fine-scale elements, that is, Comparing Equation (43) with Equation (44), the following expression is established Substituting Equation (31) into Equation (45) results in the following expression:

Equation (46) holds if and only if
Similarly, the following expression must hold for the RVE reduced order nodal load vector, that is, The reduced order RVE stiffness and nodal load matrices defined in Equations (47) and (48) can be assembled using a direct stiffness method to eventually derive the reduced order global equilibrium equation where and A denotes a standard direct stiffness assembly operator.

Solution at the coarse scale
Having computed K el, m(i) , the multiscale basis functions can be evaluated by assembling K m and solving Equation (33). The global system of equations at the coarse scale has already been established in Equation (49): where K M and f M denote the 2n M × 2n M global coarse stiffness matrix and 2n M × 1 global coarse load vector, respectively. The global coarse stiffness is assembled from its local coarse stiffness contributions as shown in Equation (50). The global coarse load vector is evaluated from the local coarse element load contributions. Finally, the EMsVEM solution is performed for the n M coarse-scale nodes u M . The microdisplacements can be evaluated from the solution of the reduced order solution of Equation (51) using the following down-scaling procedure. The coarse elementwise displacements are first extracted from u M and stored in the vector of macroelement nodal displacements u M( ) , = 1 … n M el . The displacements associated with the ith fine element in the th coarse element/RVE is computed using Equation (31). The strains and stresses associated with these fine-scale displacements can be computed as follows: where B is the term provided in Equation (23). These strains and stresses are uniform over the ith fine-element domain.
To allow for compatibility with traditional postprocessing routines, the element is decomposed into subtriangles. The evaluated stresses and strains are then associated with desired quadrature integration points. The process flow of the EMsVEM is graphically shown in Figure 7.

APPLICATIONS
The EMsVEM is verified by comparing against the standard VEM and analytical solutions. A first-order VEM (k = 1), is used in all cases. In the following, the four different element-types illustrated in Figure 8 are used for the verification. The Centroidal Voronoi Tessellations (CVT) shown in Figure 8B were generated by Lloyd's algorithm 65 with the generated seeds being forced to coincide with the associated centroid of each polygon. The Random Voronoi elements (RAND) shown in Figure 8D were created by a random set of seeds. The polygonal meshes were generated using PolyMesher. 66 To investigate the fidelity of the proposed method, the  2 norm and the  1 semi-norm are employed for the displacement and stress/strain approximations, respectively.

Square plate under tension
The case of the homogeneous square plate shown in Figure 9A is considered herein. The plate is fully clamped at the bottom and subjected to a traction vector t = A mesh discretization of 40 × 40 quadrilateral plane stress elements with full integration is employed for the FEM and VEM as shown. The EMsFEM and EMsVEM solutions are derived considering a coarse mesh consisting of 5 × 5 quadrilateral coarse elements. This is illustrated by the coarse grid in Figure 9B. Each coarse element contains 8 × 8 quadrilateral fine elements. The case of rectangular elements only is considered in this case for comparisons between the virtual element and finite element based methods to be meaningful. Linear boundary conditions are used to derive the multiscale basis functions for this example.
The  2 -norm of the errors between the FEM and the EMsFEM is 5.7 × 10 −3 . The corresponding norm for the VEM to EMsVEM comparison is practically the same, that is, 5.5 × 10 −3 . A convergence study is also performed by retaining the same number of coarse elements as shown in Figure 9B, and increasing the number of fine-scale elements.
The evolution of the  2 -norm as a function of the number of elements is provided in Figure 10, where the EMsFEM and EMsVEM behave in an identical fashion. The reference slope is also provided for the comparison of the error convergence rates; near optimal rates are observed for both methods.

Cantilever beam subjected to parabolic traction
The cantilever structure shown in Figure 11 is fully clamped at the left and is subjected to a parabolic traction at its free end. The domain has a length L = 8 m, height D = 4 m, and thickness t = 1m. The material has a Young's modulus E = 10 7 N/m 2 and a Poisson's ratio = 0.3. The parabolic traction at the free end assumes the following form where p = −1000 N is the total load applied and I = tD 3 ∕12.
The displacements are evaluated analytically on the basis of plane strain assumptions 67 according to Equations (55a) and (55b) for the horizontal and vertical component, respectively where E = E 1− 2 and = 1− . Three different geometries of microstructure elements are considered; these are summarized in Table 1. The problem is also solved with the EMsFEM for the case of quadrilateral microelements.
The convergence behavior of the EMsFEM and the EMsVEM with respect to different element-types specified in Table 1 is investigated in the form of relative  2 and  1 error convergence plots shown in Figures 12 and 13. The evolution of the errors are studied as a function of the number of coarse elements, for five microstructure configurations. At low discretizations ( Figures 12A,B and 13A,B), the EMsFEM solution provides the best accuracy. However, for finer microstructure configurations ( Figures 12E and 13E), the CVT , RAND, and QUAD element meshes used by the EMsVEM offer accuracies approaching the QUAD EMsFEM method, over all coarse-element mesh discretizations.
The convergence rates are provided in Tables 2 and 3. These are found to nearly coincide with the expected theoretical slopes of -2 in the  2 and -1 in the  1 relative error norms, respectively. The theoretical convergence rates are derived in Reference 68. One can conclude that near-optimal convergence rates are obtained by the method over all coarse-element discretizations. This suggests that the EMsVEM is a viable alternative to the EMsFEM when a flexible mesh generation is required to account for fine-scale heterogeneities.
The primary advantage offered by the EMsVEM over the EMsFEM lies in its ability to handle any kind of microstructure configuration. This is illustrated by computing displacements with the EMsVEM for a 10 × 10 coarse-element discretization with an arbitrarily chosen microstructure definition described in Table 4. The contour plots of the resulting total displacements are shown in Figure 14. Figure 14A,B illustrate the method's ability to sufficiently handle widely varying microstructural configurations within a single problem. The corresponding y-displacements obtained for the free end of the neutral axis of the structure are compared against the analytical solution in Table 5.

Cantilever beam with a periodic microstructure
In this example, a cantilever beam with periodically repeating circular inclusions is considered. A 30 × 6 coarse grid is created over this domain as shown in Figure 15A. The microstructure enclosed within each coarse element contains a This example is provided to establish the EMsVEM as a useful tool in the analysis of composites for driving down computational costs. To prove this claim, we draw attention to the relaxed conformity requirements on polygonal meshes when compared with conventional quadrilateral finite element meshes. This flexibility is exploited to minimize the number of nodes involved in the microstructural discretization while still retaining satisfactory accuracy. In particular, the microstructure is discretized using three approaches, that is, uniform quadrilateral elements, uniform polygonal CVT elements, and an adaptively refined mesh as shown in Figure 16. Mesh 1, schematically depicted in Figure 16A, comprises 10 116 uniform quadrilateral elements and is treated as a reference solution. Mesh 2, contains 5000 uniform CVT elements and is shown in Figure 16B. Mesh 3 ( Figure 16C) is adapted from Reference 69, wherein efficient polygonal discretizations are exploited for performing a nonlinear analysis on fiber composites. The properties of each mesh are summarized in Table 6. Periodic boundary conditions are used to derive the multiscale basis functions, for all cases.
The  2 and  1 error norms obtained by Mesh 2 and Mesh 3 are summarized in Table 7. It is observed that both meshes achieve comparable accuracy despite the fact that Mesh 2 has a considerably larger number of DoFs. The computational time for all cases, averaged over five runs is also shown in Table 7. This result illustrates appreciable benefits attainable through the EMsVEM by using more flexible mesh generation capabilities.
Contour plots of the total displacements are shown in Figure 17A,B for the EMsVEM-P model with Mesh 3 and the standard VEM, respectively. The corresponding von-Mises stresses VM are provided in Figure 18A,B for the EMsVEM and the VEM, respectively. Both methods provide practically identical results, as also manifested by the  2 and  1 norms shown in Table 7.

F I G U R E 14
Total displacement contours for the arbitrarily chosen microstructure defined in Table 4. Units are in m

Cantilever beam with a highly heterogeneous material distribution subjected to parabolic traction
The case of the cantilever beam examined in Section 4.2 and shown in Figure 11 is considered here also. In this case, we investigate the effect of the material heterogeneity on the performance of the proposed EMsVEM both in terms of accuracy and computational efficiency. The results obtained from the EMsVEM are compared against the standard VEM. Comparisons are also provided against the standard FEM and the EMsFEM for the case of quadrilateral elements. To investigate the effect of the assumed boundary conditions for the evaluation of the multiscale basis functions, two variants are considered, that is, linear and periodic boundary conditions. The abbreviations of the multiscale methods used in this example are shown in Table 8.
Three cases are considered vis-a-vis the geometry of the fine scale, that is, with heterogeneities having (a) QUAD, (b) CVT, and (c) RAND shapes. Furthermore, two discretizations are examined per case, to assess the the effect of the scale separation on the accuracy of the EMsVEM. The models employed along with their corresponding coarse and fine-scale discretizations are summarized in Table 9.
Each microelement in both discretization schemes is randomly assigned a Young's modulus generated by a uniform distribution. The lower and upper bounds of the Young's modulus are considered to be E l = 1 GPa and E u = 100 GPa, respectively. The material distribution associated with discretization scheme B is illustrated in Figure 19. This material distribution within the coarse element is assumed to repeat periodically over the entire domain. The Poisson's ratio is 0.3 in all cases.  The profile of the vertical displacements along the neutral axis (y = 0) of the structure are shown for all the cases in Figure 20. For discretization scheme A, shown in Figure 20A,C,E, a significant deviation in the displacements obtained through FEM and VEM and the multiscale solutions is observed. The methods using periodic boundaries, that is, EMsFEM-P and EMsVEM-P are found to approximate the complete solutions better than EMsFEM-L and EMsVEM-L.
On the other hand, in the case of the finer discretization scheme B, as shown in Figure 20B,D,F, all methods provide practically identical results. The relative error of the EMsVEM to the VEM solution for linear and periodic boundary conditions is shown in Table 10. These are practically identical to the relative errors of the EMsFEM to FEM solution which are also shown in Table 10 for completeness.
The horizontal displacements (u x ) obtained for the microstructure using RAND elements and the discretization scheme B are indicatively shown in Figure 21A,B for the VEM and EMsVEM-P, respectively. The Von-Mises stress distribution ( VM ) obtained for all cases using the discretization scheme B, are shown in Figure 22.

Discussion on computational gains
To assess the computational effectiveness of the EMsVEM, we evaluate the time required for the assembly and inversion of the global state matrices for all the methods employed. The average times over five runs are shown in Figure 23.
There is an appreciable reduction in the time taken for assembly and matrix inversion when using multiscale methods, as evidenced by Figure 23A,B. This is to be expected as the number of nodes involved in the multiscale assembly and solution procedures are significantly reduced when compared to the standard FEM and VEM. Furthermore, the time required in the EMsVEM is significantly lower when compared to the EMsFEM. This is attributed to the fact that there are no iterative evaluations over quadrature integration points done in the case of first-order EMsVEM as opposed to the EMsFEM. In Figure 23C, postprocessing times are also compared for all methods. This is to account for the effect potential overheads might have on the computational efficiency of the methods, especially within While it is clear from Figure 23D that postprocessing is a significant factor for the multiscale methods, Figure 23C reveals that this is indeed lower than in the corresponding fine scale implementations. This is due to the fact that down scaling is performed per coarse element, hence implicitly vectorizing the corresponding strain and stress computation loops.

CONCLUSIONS
A novel EMsVEM is developed for the analysis of heterogeneous elastic structures. The novelty of the EMsVEM rests on the utilization of the VEM at the fine scale to enable flexible meshing that can efficiently capture any type of domain heterogeneity. Numerically computed multiscale basis functions upscale the heterogeneities allowing the solution procedure to be performed at the coarse scale. Subsequently, the solution is downscaled to capture local fine scale variations. The EMsVEM efficiently models complex geometries at the microscale while retaining the accuracy of the standard FEM and its EMsFEM variant. The proposed method is validated against analytical solutions and the standard FEM and VEM. Our convergence studies demonstrate that the EMsVEM is well behaved even in cases of significant heterogeneities in the fine scale.
Based on the computational advantages of the VEM, stemming from its flexible element geometries and the relaxed requirement on domain integration, the EMsVEM is shown to significantly reduce the complexity of the computational problem compared to the standard FEM and the EMsFEM. The computational upscaling procedure proposed is shown to significantly reduce the computational costs compared to the standard VEM.
The ability of the method to handle adaptively refined virtual element meshes at the fine scale is also demonstrated. According to our observations, intelligently chosen polygonal fine scale discretizations drive down computational costs without significantly sacrificing accuracy.
Similar to the standard EMsFEM, the boundary conditions used to derive the multiscale basis functions play a critical role in the accuracy of the method. Periodic boundary conditions are shown to offer a better account of the heterogeneous behavior when significant heterogeneities are present.
A limitation of the proposed method is that it is confined to the utilization of quadrilateral coarse element geometries. Extending the method to account for the most general case of polygonal coarse element representations is currently under development.

C. EVALUATION OF THE VEM MATRICESG AND B
Only the expressions necessary for the derivations provided in the main manuscript are provided herein. For a more extensive discussion (see Reference 57). The weak form defined in Equation (23) is rewritten here for convenience: This is written in a more explicit form as follows: (C2) In this, the left-hand side termG can be evaluated in a straightforward manner since it contains only vector valued monomials. The right hand side term is expanded according to Equation (C3) as The domain integral B i in Equation (C3) vanishes for k = 1. When k ≥ 2, further manipulation of the integral is necessary. The term div( m (m j )) is expanded out as follows: div( m (m j )) = where m ∈ [M k−2 ( m )] 2 . The coefficients of the linear expansion d j directly depend on the material properties and can be obtained through inspection (see Reference 57). Substituting Equation (C4) into the expression for B i in Equation (C3) yields the following form where dof i in Equation (C5) is the value of u h at the ith DoF. Here, the definition of the dof operator at interior nodes (see Equation (10)) is exploited. The boundary integral B b in Equation (C3) can be evaluated directly using a Gauss-Lobatto quadrature rule according to Equation (C6) where l i is the length of the edge under consideration, w j is the Gauss-Lobatto weight of the boundary node located on that edge, N V (x j ) is the matrix of the canonical basis functions evaluated at the edge node j, N P is the matrix of C (m j ), and N e denotes the outwards looking normal edge vector written in matrix form as [ n x 0 n y 0 n y n x ] .
Hence, matrix B in Equation (C3) eventually becomes and Equation (C2) can be expressed in the convenient matrix representatioñ whereG is the left hand term in Equation (C2) and can be readily computed. For k = 1 methods, this can be fairly trivial as the integrand is comprised of constants. For higher order methods, numerical integration over subtriangulated domains might be necessary. It can also be verified using the relationG = BD as shown in. 59 m k is the vector representation of the projection operator. One solves for this: m k =G −1 B. (C10)

D. CONSISTENCY TERM
The consistency term K C shown in Equation (19) is defined as where the matrixG and m k are provided in Appendix C and n k is the number of vector valued monomials defined in Equation (22).

E. STABILIZATION TERM
To avoid hourglass modes and to ensure the coercivity of the bilinear approximation and positive definiteness of the stiffness matrix in Equation (19), an additional contribution is required that stems from the so-called stabilization term. The stability term classically employed in several works was originally proposed in. 40,56 The associated stability stiffness matrix contribution K (1) S is evaluated according to Equation (E1) as where the factor | m | has been included to ensure scaling of the stability energy with the element size. However, it has been reported in Reference 68 that this choice of the stability term may lead to nonoptimal error convergence rates in the case of three dimensional higher order methods. An alternative strategy, the D-recipe stabilization, was proposed in Reference 63. This has been adapted into an engineering context 64 and reads as follows: where is a modified stability parameter, that is, where = 1 and * = | m | tr(C) tr(D T D) .