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.


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: In the first formulation, (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), it is assumed that K is symmetric and positive definite on V . A matrix K ∈ R n×n is positive definite on a subset V ⊂ R n if v T Kv > 0 for all nonzero v ∈ V . Moreover, a symmetric positive definite matrix K constitutes a scalar product u, v = u T Kv on V , a property that will be used later. In Fig. 1 three examples of networks are shown. The network in Fig. 1(a) 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 Ku = F , some degrees of freedom have to be prescribed, resulting in the restricted solution space V . The network in Fig. 1(b) can represent a conductive medium, governed by the discrete Poisson equation. The temperature at each node is contained in u. The network in Fig. 1(c) 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.
(a) (b) (c) Figure 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]).
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 Find u ∈ V ms : v T Ku ms = v T F, ∀v ∈ V ms .
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 Ku = F is not necessarily solvable without prescribing some degrees of freedom. Consider fixed constraints with zero displacement (non-zero displacement is treated in Section 3.4) and let N D ⊂ {1, . . . , n} be the set of indices corresponding to the fixed degrees of freedom. Let N = {1, . . . , n} \ N D . Denote by V ⊂ R n the restricted solution space defined by The variational formulation of the network displacement problem reads: 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 macroscale. See Fig. 2(a) 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.
(b) A bilinear basis function Λi : (c) The interpolated bilinear basis function λi ∈ VH ⊂ R n . 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. 2(a). 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 the coarse nodal basis functions of the grid T . For a quadrilateration, bilinear basis functions are suitable, illustrated in Fig. 2(b).
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: i ∈ M D if ∃j ∈ N D : Λ i (p j ) = 0 and i ≡ j (mod d). (4) 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.
, 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. 2(c) for an illustration of the interpolated vector λ i of the coarse nodal basis function Λ i . The coarse space is defined as the span of the interpolated non-fixed basis functions, that is  [4]). With C H = B T H , an equivalent definition of the coarse space V H is as the range of the map B H C H : V → V H , that is The detail space is defined as the null space of the restriction matrix: With C H = B T H it holds that v ∈ W if the bilinear weighted average of v is zero for each interpolated bilinear basis function λ i , i.e. λ T i v = 0, ∀i ∈ M. The coarse space and the detail space constitute a splitting of V such that each v ∈ V can be uniquely decomposed as v = v H + w where v H ∈ V H and w ∈ W . 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.
To prove the second part, 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 :

Proposition 1. If B H has linearly independent columns and C
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 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.
The multiscale solution, u ms ∈ V ms , to the original problem (3), is defined as the solution to the problem Proposition 3. Let K be symmetric and positive definite on V , and F ∈ V . Then there exists a unique solution to problem (5).
where u ms is the solution to the multiscale problem (5), solves the original problem (3).
Proof. Same arguments as used in Proposition 3 show that there exists a unique such u f . According to Proposition 2, v ∈ V can be decomposed as v = v ms + w, 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 (5), 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 Proposition 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 {λ i − φ i } i∈M constitute a basis for V ms .
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 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 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 (5) 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 Sec. 6, the modifications φ i decay fast away from its coarse node. Therefore the problems (6), of calculating φ i , can be localized with preserved convergence rates. The localization is accomplished by solving each problem (6) 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  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. 6(a). Another choice is to use a fixed number of layers of coarse elements surrounding the considered element E, as depicted in Fig. 6(b). 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 tõ 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 (6). Secondly, the attained modified basis is used to solve the global multiscale problem (5). In this section, elementwise assembling and localization, described in Sec. 3.2, is employed.
First, a useful matrix notation is introduced. Consider a matrix A ∈ R a×b . Let A ⊂ {1, 2, . . . , a} and B ⊂ {1, 2, . . . , b} be subsets of the matrix row and column indices respectively. A new matrix A(A, B) ∈ R |A|×|B| is extracted from A by only considering rows corresponding to indices in A and columns corresponding to indices in B, that is A(A, B) = (a ij ) (i,j)∈A×B . 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 The full modifications are thereafter calculated according tõ and used to assemble the modified basis matrix Finally, the following localized version of the global problem (7) 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 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 (10) 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 (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 (11). 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 (13) into (12) and using (14) 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 quasi-uniform finite element mesh with mesh parameter H ≈ m −1/2 . The following two norms on the space V are introduced: Since K is symmetric and positive definite on V the smallest eigenvalue of K, denoted ν, fulfils It is assumed that ν 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, the L 2 -norm and a Poincaré inequality. A bilinear weighted interpolant π H : which is expected to hold for slowly varying F . 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 This error is viewed as acceptable 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 7. 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 Ku 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 all test functions v f ∈ W . For since v T f Ku 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 Sec. 6 it is shown, through numerical validation, that this relation is true for more complicated network models. The theoretical analysis of this result for the general network model proposed in Sec. 5 is complicated and postponed to future work. Assuming this result gives the following error bound.
Theorem 1. Assuming equation (17), the following error bound is true for any load vector F ∈ V : Proof. First consider the case F ∈ V H . The vectors u ms ∈ V ms andũ ms ∈Ṽ ms solve SinceṼ ms ⊂ V , Galerkin orthogonality gives Moreover, since F ∈ V H it holds that u = u ms . Using this together with the assumption (17) with w = u = u ms leads to |||u ms −ũ ms ||| ≤ min v∈Ṽms |||u ms − v||| ≤ Ce −cρ |||u 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 −Ku + 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 ij and assume that the edge has a width w ij . All edges is assumed to have a uniform thickness z in the direction into the plane. The direction vector, d a ij , 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 ij and given by In Fig. 7 some of the introduced notation is depicted. The length change ∆L ij of an edge, which is not exact but approximated by taking the dot product of the displacement difference onto the direction vector, is illustrated.  Figure 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 ij of an edge (i, j) is calculated.

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 ij . Consider an edge (i, j) as shown in Fig. 7(b). 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

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. 7(a). When the nodes are displaced, an angular change ∆θ ijl 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. Figure 8: The initial position of the edge pair (i, j, l) is shown with dashed lines. After displacement, forces act at the nodes of the edge pair to restore the angular change between the edges. The angle deviation of each edge is denoted δ ji and δ jl , respectively.

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. 7(a). The forces acting at the outer nodes a ∈ {i, l} will be and at the central node

Assembling of elasticity matrix
The governing equation −Ku + F = 0 contains all node equilibrium equations as a matrix system. 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 ij , K II ijl , K III ijl ∈ 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 1, 2a}, {1, . . . , n})u = F I a (i, j), a ∈ {i, j}, 1, 2a}, {1, . . . , n})u = F III a (i, j, l), a ∈ {i, j, l}.
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.  9(a). 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 1 andF correspond to the modification of K and F by explicitly including the boundary conditions.
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. 9(b) 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 Section 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. 6(a). 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 (19) are seen in Fig. 11.
for the two problem types and their three different setups are shown in Fig. 12 and Fig. 13. In some setups the errors can be compared with the so-called FEM-error, corresponding to solving the multiscale problem (8) 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.      (19). 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.      (20) with correction as in Section 13. The coarse mesh size H and its square H 2 are included in the plots to clarify convergence rates.

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 (17) 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 [14,19,12] to study virtual paper sheets and macroscale mechanical properties such as tensile strength, tensile stiffness, bending resistance, z-strength and fracture propagation.