Numerical upscaling of discrete network models

In this paper a numerical multiscale method for discrete networks is presented. The method gives an accurate coarse scale representation of the full network by solving sub-network problems. The method is used to solve problems with highly varying connectivity or random network structure, showing optimal order convergence rates with respect to the mesh size of the coarse representation. Moreover, a network model for paper-based materials is presented. The numerical multiscale method is applied to solve problems governed by the presented network model.


Introduction
Network structures are used to model a wide variety of phenomena, such as flow in porous media, traffic flows, elasticity of materials, body deformation in computer graphics, molecular dynamics, and fiber materials. In these applications, the microscale behaviour determines the macroscale properties of the system. Often a full microscale model is difficult or impossible to work with because of the vast computational complexity. Therefore, there is an interest in constructing coarser, but still accurate, representations of the entire system. Such a procedure is sometimes referred to as upscaling or homogenization. In this work a numerical upscaling method for discrete networks is presented.
There exist several numerical upscaling methods for partial differential equations (PDE) based on the idea of homogenization, such as the heterogeneous multiscale method (HMM) [20], the multiscale finite element method (MsFEM) [8], and the more recent works [3,17]. The upscaling approach presented in this paper is based on the localized orthogonal decomposition method (LOD) [4,15], which in turn is inspired by the variational multiscale method (VMM) [9]. Multiscale methods applied to network problems are for instance investigated by Ewing [5] and Ilev et al. [11] who study the heat conductivity of network materials and develop an upscaling method by solving the heat equation locally over small sub-domains. These local solutions are used to compute an effective global thermal conductivity tensor. Della Rossa et al. [2] investigate network models of traffic flows and derive a governing PDE for the macroscale by formulating traffic flow equations for single network nodes and interpreting the relations as finite difference approximations. The macroscale parameters are resolved using a two-scale averaging technique. Chu et al. [1] develop a multiscale method for networks representing flows in a porous medium. The medium is modelled as a network where nodes represent pores and edges represent throats. The conductance of each throat is assumed to be given by Hagen-Poiseuille equation, and using mass conservation equations for the flow through the network, a model for the microscale is attained.
The numerical upscaling method proposed in this work is developed for general unstructured networks. The network is supposed to represent the microscale, and the macroscale is represented by a finite element mesh which is coarse in comparison to the fine scale network. The coarse grid does not have to be related to the network in any way except that both cover the same computational domain, and therefore the method can be applied to arbitrary network geometries. The coarse FEM grid is used to define a macroscale solution space spanned by basis functions defined at each coarse grid node as in standard FEM. The upscaling idea is to modify the coarse basis functions to account for the microscale features of the network. This is accomplished by solving local sub-network problems at each coarse basis function. The modified basis functions are thereafter used to solve a global low-dimensional system resulting in an accurate macroscale solution. The method leads to modified basis functions that decay exponentially, and hence localization of the local sub-network problems can be utilized, reducing the computational cost considerably while preserving optimal convergence rates.
Moreover, this paper includes a two-dimensional network model, which can be used to model paper-based materials in form of fiber networks. The macroscale mechanical properties of paper-based materials are of great interest. Paper is a heterogeneous material built up of fibers bonded together into a network structure. The mechanical properties of paper depend primarily on the properties of the fibers and the bonds between them. In [12,14,19], computational fluid dynamics and advanced contact modeling are used to simulate the paper forming process. One future aim is to utilize that framework together with the proposed multiscale method to create virtual fiber networks and investigate the macroscale mechanical properties. A network representation including fibers and bonds is a suitable methodology to study the mechanical properties of paper [10,13,18]. Moreover, the varying properties of single fibers and bonds, as well as an interest for fracture propagation simulations, call for an upscaling approach. The presented network model is based on forces arising at the nodes when the network is displaced, acting to restore the initial configuration. The network model is similar to lattices models like [16,21] where edges are represented by springs. Moreover, angle springs between pair of edges are included. A novelty of the network model in this work is a third type of force phenomenon resulting in an effect similar to the Poisson effect. Force equilibrium equations at each node result in a matrix equation which can be very large. For a regular network, the model converges to the linear elasticity equation when the length of the network edges tends to zero. The numerical upscaling method is applied to the network model and numerical examples are solved to demonstrate the convergence rates of the method. The examples show how the proposed numerical upscaling method resolves fine scale features which the standard FEM cannot.
The outline of this text is as follows. In Sect. 2, the general problem formulation is stated. Thereafter, in Sect. 3, the theory of the numerical upscaling method is presented. Sect. 4 contains error analysis, and in Sect. 5, the two-dimensional network model is described. In Sect. 6, numerical examples are presented, showing the convergence rates of the proposed method. Lastly, in Sect. 7, conclusions and future work are discussed.

Problem formulation
Consider a problem modelled by a network with properties governed by a connectivity matrix K ∈ R n×n . The matrix K can for instance be the discrete Poisson operator describing heat conduction, the finite difference discretization of the linear elasticity operator, or represent a more complex model, such as of the mechanics of a fiber network. Let F ∈ R n denote the load vector and let the solution vector be denoted u, belonging to a vector space V ⊂ R n . The network problem can be stated in two equivalent ways, either: or: Fig. 1 Three examples of networks: a regular square network (a), a regular square network with randomly perturbated nodes (b), and a fiber network (c) (generated as in [12]) In the first formulation, (2.1),K andF denotes modifications of K and F by explicitly including the restriction of u to the space V , for instance by holding some nodes fixed.
To ensure existence and uniqueness of the second formulation, (2.2), it is assumed that K is symmetric and positive definite on Moreover, a symmetric positive definite matrix K constitutes a scalar product u, v = u T K v on V , a property that will be used later. In Fig. 1 three examples of networks are shown. The network in Fig. 1a exemplifies a finite difference grid for the unit square, with K as the resulting discretization of the linear elasticity operator. This problem setup can be used to find the node displacements u under applied node forces F. To attain a solvable system K u = F, some degrees of freedom have to be prescribed, resulting in the restricted solution space V . The network in Fig. 1b can represent a conductive medium, governed by the discrete Poisson equation. The temperature at each node is contained in u. The network in Fig. 1c is a fiber network building up a paper sheet. The fibers are modelled as chains of edges connected at nodes with bonds between fibers at common network nodes.
The objective of this paper is to develop a numerical upscaling method for networks, circumventing the computational issues arising when materials of macrosize are considered. The idea is to reduce the size of the system by introducing a subspace V ms ⊂ V , as a coarse representation of the network. This space is called the multiscale space and it should fulfil the condition that dim V ms is much lower than dim V . The multiscale solution is attained from the problem The aim is to construct a multiscale space such that an error u − u ms is small. To achieve this, a FE-type coarse space is first introduced, which does not have the desired approximation properties. This coarse space is then modified by solving local sub-networks problems, resulting in the desired multiscale space. In the following section such a numerical homogenization method is presented.

Numerical homogenization of networks
Consider a network with N nodes and properties governed by a symmetric and positive semi-definite matrix K ∈ R n×n , where n = d · N is the number of degrees of freedom of the network, and d denotes the number of degrees of freedom at each node. For instance, for an elastic network, where node displacements are to be solved, the number of degrees of freedom at each node will be two or three, depending on if the network is two-or three-dimensional. In the following presentation the space is assumed to be two-dimensional, but the method works analogously for three dimensions. Denote by p i ∈ R 2 the position of the node corresponding to degree of freedom i = 1, . . . , n. Note that groups of d degrees of freedom correspond to the same position.
Let the solution vector be denoted u ∈ R n . The ordering of nodes and their degrees of freedom is arranged such that if d = 2, u(1) and u (2) correspond to the first and second degree of freedom of node 1, u(3) and u(4) correspond to the first and second degree of freedom of node 2, and so on, with analogous ordering if d is larger. Here u(i) denotes the i:th component of vector u. Let F ∈ R n denote the load vector. The system K u = F is not necessarily solvable without prescribing some degrees of freedom. Consider fixed constraints with zero displacement (non-zero displacement is treated in Sect. 3.4) The variational formulation of the network displacement problem reads: (3.1) For the problem to be solvable it is assumed that K in addition to being symmetric, also is positive definite on the restricted solution space V .

Coarse grid representation
The overall idea of the upscaling method is to introduce a coarse grid, representing the network at the macroscale. See Fig. 2a for an illustration of a network with a coarse grid representation. At each coarse node, d number of basis functions are defined similarly as in the finite element method. These basis functions span a low dimensional solution space which gives an insufficient description of the fine scale features. To include the fine scale information, the basis functions are modified by solving local sub-network systems. Thereafter the modified basis functions are used to solve a global system, smaller than the full system including all nodes, resulting in an upscaled approximation of the original problem. In what follows, the details of this procedure are described. Let the coarse grid be denoted T , containing M coarse nodes and let m = d · M be the degrees of freedom of the coarse grid. One choice of coarse grid is a quadrilateration as in Fig. 2a. It is assumed that the coarse grid constitutes a good approximation of the computational domain of the network, and that each coarse element contains at least one network node and that N > M. Let Λ i : R 2 → R, i = 1, . . . , m, denote   Let M D ⊂ {1, . . . , m} be the set of indices corresponding to fixed coarse degrees of freedom. The fixation of coarse grid nodes is determined from the set of fixed network nodes, N D . Consider a coarse node with basis function Λ i , describing for instance the x-displacement of that node. If there exists a network node with fixed degree of freedom j such that p j lies in the support of Λ i , and j also describes x-displacement, then the coarse degree of freedom i should be fixed. This is illustrated in Fig. 3. The fixation condition is equivalently stated as: The modulo operation is used to check if two degrees of freedom, here i and j, correspond to the same displacement direction (e.g. x or y for d = 2). For networks with more complex boundary geometry, the coarse grid has to be refined at the boundary to attain a proper representation of the fixed boundary conditions.
Let M = {1, . . . , m}\M D denote the set of nonprescribed coarse degrees of freedom. The positions of the coarse nodes, {P i } m i=1 , defined similarly as the positions of the network nodes, are a subset of R 2 , likewise as the nodes of the network. However, these two subsets do not have to be related, but as already noted, it is assumed that each coarse element contains at least one network node.
Next, two vector spaces are introduced, the coarse space V H , and the detail space W . The coarse space is defined from the coarse basis functions in the following way. Let λ i ∈ R n , i = 1, . . . , m, be the interpolation of the coarse nodal basis functions to the network nodes given by See Fig. 2c for an illustration of the interpolated vector λ i of the coarse nodal basis The coarse space is defined as the span of the interpolated non-fixed basis functions, that is The detail space is defined as the null space of the restriction matrix: The coarse space and the detail space constitute a splitting of V such that each v ∈ V can be uniquely decomposed as Before proving this fact, a lemma is stated showing the relation between the spaces R m H , R n and V H , and the maps in-between, illustrated in Fig. 4.

Lemma 3.1 If B H has linearly independent columns and C H = B T H , then for each
To prove the second part,

Proposition 3.1 If B H has linearly independent columns and C H
Using the detail space W , together with the connectivity matrix K , the multiscale space V ms is defined as the K -orthogonal complement of W : The spaces W and V ms constitute another splitting of V implying that every v ∈ V can be decomposed uniquely as v = v ms + w where v ms ∈ V ms and w ∈ W .

Proposition 3.2 Assume B H has linearly independent columns and C H = B T
H . If K is symmetric and positive definite on V , then V = V ms ⊕ W uniquely.
which has a unique solution since K is symmetric and positive definite on V . Define v ms = v H − z and The multiscale solution, u ms ∈ V ms , to the original problem (3.1), is defined as the solution to the problem where v ms ∈ V ms and w ∈ W . Using orthogonality, and that u f and u ms are solutions to their respective problem, it can be derived that To solve the multiscale problem (3.3), it is convenient to construct a basis for the multiscale space V ms , which can be used to simplify the problem to a matrix equation. A basis for V ms is constructed using the vectors λ i , by defining modification vectors φ i ∈ V , i ∈ M , as solutions to the problems (3.4)

Proposition 3.5 If K is symmetric and positive definite on V , B H has linearly independent columns and C H = B T H , then the vectors {λ
Proof The problem to find φ i such that w T K φ i = w T K λ i , ∀w ∈ W , has a unique solution φ i ∈ W since K is symmetric and positive definite on V . By construction, it is also true that λ i −φ i ∈ V ms . To prove linear independence, consider a i (λ i −φ i ) = 0 and apply B H C H to both sides. Using that C H φ i = 0 gives B H C H a i λ i = 0 = B H C H 0, and by the second part of Lemma 3.1 it follows that a i λ i = 0 implying a i = 0.
Using the above constructed basis for V ms , {λ i − φ i } i∈M , to assemble the matrix reduces the variational form of the multiscale problem (3.3) to the equivalent matrix problem where U ms ∈ R m H and u ms = B ms U ms . At this point, a multiscale space with low dimension compared to the solution space V has been constructed as was desired in the problem formulation. Moreover, it has been shown that the problem can be solved through a matrix equation after constructing a basis for the multiscale space. However, the problems of solving the modified basis functions have the same size as the original problem. It turns out that this can be circumvented, utilizing that the modifications φ i decay exponentially, implying that the problems can be localized. This is presented in the following section.

Localization
The described method requires systems to be solved which are as large as K . However, as will be demonstrated by numerical examples in Sect. 6, the modifications φ i decay fast away from its coarse node. Therefore the problems (3.4), of calculating φ i , can be localized with preserved convergence rates. The localization is accomplished by solving each problem (3.4) on a restricted domain, called patch.
As in FEM, it is natural to assemble the stiffness matrix elementwise. In this work, it is suitable to assemble the connectivity matrix K over each coarse element E, such that K = E∈T K E . The local element connectivity matrices K E : V → V are assembled for each element E ∈ T by only considering edges in each element. See Fig. 5 for an illustration. For unstructured networks, edges may intersect the elements. Such a situation is resolved by temporarily dividing the edges at intersection points. Using the decomposition of the connectivity matrix into local element matrices, the modifications φ i can analogously be assembled as E∈T Proposition 3. 6 The sum of the elementwise modifications, E∈T φ E i , solves the original problem (3.4).
Proof It follows that  With elementwise assembly, it is convenient to use patches centred at each element. For an element E ∈ T , let its patch be denoted ω E ⊂ R 2 . One suitable choice of patch geometry is a circle with center coinciding with the element center, as illustrated in Fig. 6a. Another choice is to use a fixed number of layers of coarse elements surrounding the considered element E, as depicted in Fig. 6b. Let the patch size be described by the parameter ρ such that ρ H is the radius of the patch, where H denotes the coarse element size. In the two examples shown in Fig. 6, the value of ρ corresponds to 1.5. Let N E ⊂ N , denote the degrees of freedom of network nodes that are in the patch ω E . Similarly, let M E ⊂ M , denote the degrees of freedom of coarse nodes that are in the patch of element E. For each element E, its localized subspace of the detail space is defined according toW The localized modificationφ i is attained by solving the problems for each element E, and taking the sumφ i = E∈Tφ E i . The resulting localized multiscale space is denotedṼ ms , and is defined as

Algebraic formulation
In this section the algebraic formulation of the numerical multiscale method is presented. The multiscale method consists of two main steps. First, the basis for the multiscale space V ms is constructed by calculating each φ i from (3.4). Secondly, the attained modified basis is used to solve the global multiscale problem (3.3). In this section, elementwise assembling and localization, described in Sect. 3.2, is employed.
First, a useful matrix notation is introduced. Consider a matrix A ∈ R a×b . Let Using the introduced matrix notation, the following matrices and vectors are defined for E ∈ T and i ∈ M : The multiscale method can now be formulated as solving several matrix systems. For each element E ∈ T , and each coarse degree of freedom i ∈ M , the following system is solved Finally, the following localized version of the global problem (3.5) is solved: where the fine multiscale solution is calculated asũ ms =B MŨms . Summarized, the algebraic formulation of the numerical multiscale method is:

Non-zero fixed boundary conditions
Consider a displacement problem where F = 0 and the set of prescribed degrees of freedom, N D , includes degrees with non-zero displacement. Let (3.8) Assume u = u 0 + g H where u 0 ∈ V and g H for simplicity is a linear combination of λ i , i = 1, . . . , m. The problem (3.8) is then equivalent to The problem is in this way transformed back to the previous formulation and can be solved with the same method.
Proposition (3.4) states that the exact solution of the displacement problem is u = u ms + u f . For the case with non-zero prescribed displacements and g H ∈ span({λ i } m i=1 ), which was described above, the correction vector u f turns out to be a linear combination of the vectors φ i , i = 1, . . . , m, as the following derivation shows. Consider the non-homogeneous displacement problem (3.9). The correction is attained from where g H was assumed to be a linear combination of λ i , i.e.: The modification vectors φ i , i = 1, . . . , m, is attained from the problems Note that the problems are solved for all i = 1, . . . , m, compared to before, when only i ∈ M were considered. This is necessary to construct the correction u f as will be seen next. Inserting (3.11) into (3.10) and using (3.12) gives implying that the correction is the sum For general fixed displacements g H , not necessarily in V H , see the work [6].

Error analysis
This paper concerns a quite general network model described by a connectivity matrix K and a right hand side load F. In this section, some error bounds are shown. Because of the generality of K , assumptions are needed. It is assumed that the coarse grid is a quasi-uniform finite element mesh with mesh parameter H ≈ m −1/2 . The following two norms on the space V are introduced: , whereh i is the average length of all edges connected to the node corresponding to degree of freedom i. It is assumed that K is symmetric and positive definite on V , and that the smallest eigenvalue of K , denoted ν, fulfils where ν is bounded from below by a constant independent of n. For the finite element method posed on a quasi-uniform mesh these definitions and assumptions correspond to the energy norm, a weighted l 2 -norm and a Poincaré inequality. A bilinear weighted which is expected to hold for v ∈ V . Introducing the notationF, with componentsF i = F ih −2 i , the following variant of the Cauchy-Schwarz inequality can be derived using the original version, The weighted load vectorF is O(1) in regard to n for distributed surface loads.
The main source of error in the proposed method is the localization of the multiscale basis functions to patches. To isolate this contribution the vector π H F ∈ V H is introduced. Consider the modified model problem: Given the assumptions, it is noted that which together with (4.2) implies This error is viewed as acceptable since the constant C is independent of data variation in the connectivity matrix, and without loss of generality it is assumed that F ∈ V H . First it is shown that the multiscale solution u ms is equal to u if F ∈ V H .

Proposition 4.1 Let u and u ms be defined according to
If F ∈ V H , then it holds that u = u ms .
Proof Since F ∈ V H , there exists anF ∈ V such that B H C HF = F. To show v T K u ms = v T F, ∀v ∈ V , it is noted that it holds for all v ∈ V ms from the definition of u ms . Due to the splitting V = V ms ⊕ W it is enough to show that it also is true for since v T f K u ms = 0 by orthogonality and B T H v f = C H v f = 0 by the definition of the space W . Therefore u ms solves the same equation as u, and due to uniqueness u ms = u.
The error committed by localization is difficult to study without stronger assumptions on the connectivity matrix K . It has however been analysed for several concrete cases, for instance when K is the stiffness matrix arising from discretizing the Poisson equation [15] and the elasticity equations [7] with the finite element method. For these cases it is true that where C and c are constants independent of H . In Sect. 6 it is shown, through numerical validation, that this relation is true for the more complicated network model considered there. The theoretical analysis of this result for the general network model proposed in Sect. 5 is complicated and postponed to future work. Assuming this result gives the following error bound.
Theorem 4.1 Assuming Eq. (4.5), the following error bound is true for any load vector F ∈ V : Proof First consider the case F ∈ V H . The vectors u ∈ V, u ms ∈ V ms andũ ms ∈Ṽ ms solve Moreover, since F ∈ V H it holds that u = u ms . Using this together with the assumption (4.5) with w = u = u ms leads to |||u ms −ũ ms ||| ≤ min v∈Ṽ ms |||u ms − v||| ≤ Ce −cρ |||u ms |||.
This together with the fact that For the general case F ∈ V , letû,û ms andû ms denote the different solutions to the problems in (4.6), but with load vectors π H F. Usingû =û ms and the inequalities (4.4) and (4.7) finally result in |||u −ũ ms ||| = u −û +û ms −û ms +û ms −ũ ms

Network model for paper-based materials
In this section, a two-dimensional elasticity network model is presented, which can be used to model fiber networks in paper-based materials. An elasticity network consists of nodes and edges. When the nodes are displaced, internal forces act to restore the displacements. These forces act at two types of elements, either on edges or on edge pairs (two edges connected at a joint node). Three types of internal forces are included in this model, one type is related to edges, and two types are related to edge pairs. The model is two-dimensional, static and assumes small deformations. The network mechanics are governed by force equilibrium equations assembled at each node. The general form of the equation for each node i reads where F External i ∈ R 2 are all externally applied forces. The internal force F Internal i ∈ R 2 is, as mentioned, a sum of three contributions: The first force contribution is related to the edges of the network and acts to compensate for changes in length of the edges. The second force contribution acts on edge pairs to compensate changes in the angle between the edges of each edge pair. The third force contribution is included to model the Poisson effect. It acts on edge pairs by introducing a resistance to changes in the total length of the two edges of the pair. In the following sections, the three forces are described. Preparatory some nomenclature is introduced. Let the network consist of N nodes. Let (i, j) denote the edge connecting node i to j and let E denote the set of all edges. Note that (i, j) = ( j, i). Edge pairs are denoted by (i, j, l) where j is the central node. Denote by P the set of all edge pairs. Note that (i, j, l) = (l, j, i). Each node i has two degrees of freedom, the x-directed displacement and the y-directed displacement, contained in the vector δ i ∈ R 2 .
All force equilibrium equations can be assembled into a system of the form −K u + F = 0 where K ∈ R n×n is called the elasticity matrix, u ∈ R n is the node displacements, and F ∈ R n is the external forces. The elasticity matrix K is attained by summation of matrices assembled at edges and edge pairs. The node displacement vector u is arranged according to and the elasticity matrix is assembled such that the force equilibrium equation for node i is at row 2i − 1 for the x-component, and at row 2i for the y-component. Let the length of edge (i, j) be denoted L i j and assume that the edge has a width w i j . All edges is assumed to have a uniform thickness z in the direction into the plane. The direction vector, d a i j , of an edge (i, j), with respect to node a ∈ {i, j}, is defined as The length change of edge (i, j) is denoted ΔL i j and given by In Fig. 7 some of the introduced notation is depicted. The length change ΔL i j of an edge, which is not exact but approximated by taking the dot product of the displacement difference onto the direction vector, is illustrated.

Extension of edges
The first force contribution acts at edges due to their internal resistance to length change. When the nodes of an edge are displaced so that the projection of the difference of the node displacements onto the initial edge direction is nonzero, anti-parallel forces arise at the nodes of the edge to restore the length. The tendency of an edge to restore its length is described by the elastic modulus k i j .
Consider an edge (i, j) as shown in Fig. 7b. When the nodes are displaced, the edge (i, j) will give rise to two forces F I a (i, j), a ∈ {i, j}, acting on node a according to (a) (b) Fig. 7 Sketches showing the notation used for the network nodes, edges and edge pairs. To the right it is illustrated how the approximate length change ΔL i j of an edge (i, j) is calculated

Angular deviations of edge pairs
The second force contribution acts at edge pairs from their internal tendency to resist change of the angle between their two edges. When a change in angle occurs, two torques arise at the connecting node acting on one edge each to restore the change. By transforming these torques to force couples the effect can be converted into the force equilibrium equations. Consider an edge pair (i, j, l) as depicted in Fig. 7a. When the nodes are displaced, an angular change Δθ i jl occurs, giving rise to two torques acting on edge (i, j) and ( j, l) respectively, at the position of node j. The angular change is a sum of two contributions according to Each term, δ ji and δ jl , is the angle deviation of respective edge from its initial orientation, as can be seen in Fig. 8. By using the assumption α ≈ tan α, the angles δθ ja , a ∈ {i, l} can be calculated according to where the edge normals are calculated according to  Transforming the torques to force couples gives the resulting three forces F II a (i, j, l), a ∈ {i, j, l} from edge pair (i, j, l), acting on node a, according to These forces are illustrated in Fig. 8.

Poisson effect of edge pairs
The third force contribution results in an effect similar to the Poisson effect and acts at edge pairs. The idea is to add forces that work to keep the total length of the two edges of the pair constant. Hence, when one edge changes length, two kind of forces occur, on one hand forces acting to restore the length of the specific edge, on the other hand forces acting to change the length of the other edge in the pair. Consider an edge pair (i, j, l), as shown in Fig. 7a. The forces acting at the outer nodes a ∈ {i, l} will be and at the central node

Assembly of elasticity matrix
The governing equation −K u + F = 0 contains all node equilibrium equations as a matrix system. HereK is the modification of K by explicitly including boundary conditions. Since each edge and each edge pair leads to separate force contributions, the total elasticity matrix K can be assembled from separate element matrices for each of the three force contributions. Let K I i j , K II i jl , K III i jl ∈ R n×n denote the matrices assembled from the first, second and third force contribution respectively, at different elements [edges (i, j) or edge pairs (i, j, l)]. These matrices are sparse and the only nonzero elements are defined by the relations The elasticity matrix K is assembled according to The matrix K is symmetric and semi-positive definite. With proper fixation of nodes, the matrix will be positive definite on the restricted solution space. Moreover, for a regular network with uniform coefficients k, κ, η and γ , the presented model is equivalent to the finite difference discretization of the two-dimensional linear elasticity equations.

Numerical results
In this section, two network problems are solved using the proposed multiscale method and the fiber network model. Both problems are similar, with different boundary conditions and load vectors. Consider a unit square network with nodes and edges in a regular grid pattern, as shown in Fig. 9a. Let the number of nodes be (r + 1) 2 .
In the examples r = 2 7 = 128 will be used. Let l i and μ i represent the standard Lamé parameters, with one value for each node i. Set the parameters of the network model to , are the mean values of the two nodes of edge (i, j). Using the above parameters, the elasticity matrix K is assembled. Let h = 1/r denote the length of each network edge.
The first network problem, called the fixed boundary problem, is where u(i) = 0 for all network nodes on the unit square boundary, and F(i) = 1 h 2 . As mentioned earlier,K andF correspond to the modification of K 1 and F by explicitly including the boundary conditions. Regular square network.

(a) (b) (c)
Network with coarse scale grid representation.
Coarse scale basis function.  The second problem, called the displaced boundary problem, is where u(i) = 0 for all network nodes with x = 0, and u(i) = 0.1 for all x-directed degrees of freedom with x = 1.
To solve the two problems using the proposed numerical multiscale method, a coarse FEM grid is introduced. The grid is similar to the network but with (R + 1) 2 nodes. The basis functions Λ i are chosen as classic bilinear. See Fig. 9b for an illustration of the network and coarse FEM grid. Let H = 1/R be the width of the coarse elements. With the described network geometry, the fixed boundary conditions correspond to fixation of coarse nodes at the boundary.
Each problem is solved with three different setups, first a basic setup with l i = μ i = 1. Secondly, by using a realization of random coefficients l i and μ i sampled in [0.1, 10] with uniform distribution. Thirdly, each node i is displaced [δx i , δy i ] where δx i and δ y i is randomly sampled in [−0.4h, 0.4h] with uniform distribution. For nodes with initial coordinate x = 0 or x = 1, it is enforced that δx i = 0, and similarly for nodes with y = 0 or y = 1, δ y i = 0. In Fig. 10, a network with such random structure is shown for r = 32. The benefits of the multiscale method will be demonstrated by utilizing localization as described in Sect. 3.2, by introducing patches that the local modification problems are solved over. Which degrees of freedom that are included in each patch set N E , are chosen based on the radius ρ H . A degree of freedom i ∈ N will be in N E , if | p i − c E | ≤ ρ H , where c E ∈ R 2 is the center of element E. Hence patches will have the circular form as was illustrated in Fig. 6a. The rapid decay of the modified basis {λ i − φ i } i∈M is demonstrated by solving the modificationφ i with different ρ and computing the relative error |||φ i −φ i |||/|||φ i |||. The degree of freedom i is chosen as one of the central nodes but the trend is similar for all nodes. The resulting errors for the first problem (6.1) are seen in Fig. 11.
Moreover, the validity of the assumed Poincaré inequality (4.1) is investigated by solving the generalized eigenvalue problem and calculating the smallest eigenvalue. Here D denotes the diagonal matrix with valuesh 2 i , i = 1, . . . , n on the diagonal. This is done for the three different setups and r ∈ {2 2 , 2 3 , 2 4 , 2 5 , 2 6 }. In Fig. 12, it is seen that the resulting eigenvalues are independent of n. The values are normalized with the value for r = 2 6 .
Next, the two problems (6.1) and (6.2) are solved using r = 2 7 = 128 for different coarse grids with R = 2, 4, 8, 16, 32. The problems are solved using localization with patch radius ρ H = C H log 2 (H −1 ) = k/2 k where k = 1, 2, 3, 4, 5. For the first    13 Relative errors for the first problem (6.1). The coarse mesh size H and its square H 2 are included in the plots to clarify convergence rates. For the basic setup and the setup with random coefficients, also the so-called FEM-error is included problem (6.1), the constant is chosen to C = 1, and for the second problem (6.2) it is C = 1.5. The second problem, with non-zero fixed displacement, is solved with correction as described in Sect. 3.4. The resulting relative errors for the two problem types and their three different setups are shown in Figs. 13 and 14. In some setups the errors can be compared with the so-called FEM-error, corresponding to solving the multiscale problem (3.6) with non-modified basis B H instead of the modified multiscale basisB ms . It can be seen that this FEM-solution behaves poorly for the setup with random coefficients.

Conclusion and future work
In this paper a numerical multiscale method for discrete networks is proposed. For a set of different numerical examples, the convergence rates of the proposed method are examined. For regular networks with low connectivity variation the method resembles the convergence rates of the ordinary FEM. For networks with randomly varying connectivity it is shown that the multiscale method performs better than ordinary FEM. The method is moreover used to solve network problems with random structure, indicating error convergence rates at least linear in energy norm and quadratic in l 2 -norm.
A challenging theoretical problem for the future is to extend the result in (4.5) beyond finite element based discretization to more general networks, including the presented network model. Further on, other fiber network models, for example beam based, as well as three-dimensional, should be considered. Thereby the presented numerical upscaling method will be applied to realistic macroscale paper networks,  investigating the problem size which can be studied and the computational efficiency of the proposed method. Moreover, the proposed method will be used together with the paper forming simulation framework presented in [12,14,19] to study virtual paper sheets and macroscale mechanical properties such as tensile strength, tensile stiffness, bending resistance, z-strength and fracture propagation.