{'fraction_non_alphanumeric': 0.0388403772266853, 'fraction_numerical': 0.013133077191756898, 'mean_word_length': 4.151313422094279, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Local multiscale model reduction using discontinuous Galerkin coupling for elasticity problems

In this paper, we consider the constrained energy minimizing generalized multiscale finite element method (CEM-GMsFEM) with discontinuous Galerkin (DG) coupling for the linear elasticity equations in highly heterogeneous and high contrast media. We will introduce the construction of a DG version of the CEM-GMsFEM, such as auxiliary basis functions and offline basis functions. The DG version of the method offers some advantages such as flexibility in coarse grid construction and sparsity of resulting discrete systems. Moreover, to our best knowledge, this is the first time where the proof of the convergence of the CEM-GMsFEM in the DG form is given. Some numerical examples will be presented to illustrate the performance of the method.


Introduction
Many problems in science and engineering exhibit high heterogeneity. One example is the seismic meta material for which the Young's modulus of Foamed plate is about 1.6×10 5 Pa, and it is about 2.07×10 11 Pa for steel [28]. In order to solve problems in highly heterogeneous media, one usually needs to discretize the governing equation with high resolution mesh in order to obtain reliable solutions, and this will inevitably results in large-scale systems which are hard and sometime even impossible to solve with current computing resources. Therefore, some model reduction techniques are necessary. There are some existing model reduction techniques for handling highly heterogeneous media or material. One way is to homogenize the detailed heterogeneous media based on some rules so that one can solve the equation with low resolution grid [2,9,16,24]. Another way is to solve the equation in a coarse grid using multiscale basis functions. These functions are usually constructed by solving some carefully designed local problems, and they include heterogeneous information of the media [1,11,12].
There are in literature various multiscale model reduction techniques based on different types of multiscale methods for highly heterogeneous problems. Some examples include the multiscale finite element method [4,18], variational multiscale method [19,19], multiscale mortar methods [3,27], localized orthogonal decomposition method [22], multiscale finite volume method [10,20,21]. Among these methods, the multiscale finite element method (MsFEM) and its generalization the Generalized Multiscale Finite element method (GMsFEM) have achieved success for various types of heterogeneous problems [12,25,26]. The basic idea of both the MsFEM and the GMsFEM is to solve the equation in a relative low resolution grid with appropriately constructed multiscale basis functions instead of polynomial basis functions as in traditional finite element methods. These basis functions are constrained forms trusted by solving some local heterogeneous problems with specific boundary conditions and they include small scale medium information. Thus these basis functions allow us to perform accurate but inexpensive coarsegrid simulations. Motivated by the localized orthogonal decomposition method (LOD) [13,17,22,23] the constraint energy minimizing generalized multiscale finite element method (CEM-GMsFEM) [8] is developed. The CEM-GMsFEM basis functions are constructed based on local spectral problems and an energy minimization principle. The resulting basis functions have exponential decay property for high contrast media. Moreover, the convergence depends only on the coarse mesh size and is independent of the scale and contrast of the media. The method has shown success for many problems, for example [5,25]. The original CEM-GMsFEM is based on continuous Galerkin coupling. It is extended to the discontinuous Galerkin (DG) formulation in [6]. The DG coupling allows the basis functions to be discontinuous and thus gives more flexibility for complex applications.
In this paper, we consider the CEM-GMsFEM with discontinuous Galerkin (DG) coupling for the linear elasticity equations in highly heterogeneous media. We use the classical internal penalty discontinuous Galerkin (IPDG) method as the coarse grid scheme to couple the multiscale basis functions. The construction of multiscale basis functions consists of two major steps. In particular, we will construct auxiliary basis functions by solving local spectral problems on each coarse element, which is followed by solving energy minimization problems in oversampled domains. This will yield the final multiscale basis functions. We propose two variants of energy minimization problems. The first one is based on solving constraint energy minimization problems while the second one is the relax version by solving unconstrained energy minimization problems [14,15]. In comparison with previous convergence methods, the form of the CEM-GMsFEM under DG coupling proposed in this paper provides better convergence for elastic problems. By choosing an appropriate number of oversampling and basis functions, an accurate approximation can be achieved even with very large contrasts. Rigorous analysis and numerical examples for both versions will be provided.
This paper will be organized as follows. In Section 2, we present the isotropic elasticity problem in heterogeneous media and the IPDG method. In Section 3, CEM-GMsFEM is presented. In Section 4, we provide some stability and convergence results. In Section 5, we give some numerical results to illustrate the performance of this method. Finally, we conclude the paper in section 6.

Preliminaries
In this section, we present the isotropic elasticity problem in heterogeneous media and the DG method. We will also introduce the mesh required for our method, and some norms for our analysis.

Isotropic elasticity problem in heterogeneous media
In this paper, we consider the following elastic problem in the domain Ω ⊂ R 2 , div (σ (u)) + f = 0, (1a) where u is the displacement field, and f is the given source term. In the above system, σ = σ (u) represents the stress tensor, = (u) represents the strain tensor, C = C i,j,p,q (x) is the fourth-order tensor representing material properties where i, j, p, q = 1 or 3 for two-dimensional space. Using Voigt notation, the elastic tensor can be expressed in terms of the following coefficient matrix In addition, for isotropic elastomer, its constitutive equation is σ(u) = λtr ( (u)) I + 2µ (u), and λ and µ are Lame constants satisfying λ + µ > 0, µ > 0. So, we have C 11 = λ + 2µ, C 33 = λ + 2µ, C 13 = λ, C 55 = µ. We recall that, in (1c), grad u = ∂ui ∂xj 1≤i,j≤2 , and

Interior Penalty Discontinuous Galerkin (IPDG) method
We first introduce some concepts of grids as illustrated in Figure 1. Let T H be a coarse grid for the domain Ω. An element K in T H is called a coarse element, where K is a triangle or quadrilateral in the two-dimensional case. We call H the size of the coarse grid, and we denote the diameter of the cell K as h K . We let N c be the total number of coarse nodes, and N be the number of coarse elements. We also denote the collection of all coarse grid edges by E H . For any real number s, the discrete Sobolev space H s T H is defined as follows: The corresponding discrete Sobolev norm is defined as follows: where · s,K is the norm on H s (K). Besides, we denote the seminorm of the discrete gradient as In addition, we define the fine grid as T h , and h > 0 is the size of the fine grid. We note that the fine grid is a refinement of the coarse grid T H . In Figure 1, we illustrate an example of coarse and fine grids when rectangular partitions are adopted. In the rest of the paper, we assume this type of rectangular partition is used. Figure 1: Illustration of the grids. We use black to denote the fine grid and red to denote the coarse grid. A coarse element K is shown.
Next we give the formulation of the DG method. For each coarse block K i , we denote the restriction of the Sobolev space H 1 (Ω) on K i by V h (K i ). That is, for the discrete space, V h (K i ) is the space of piecewise bilinear functions with respect to the fine grid in K i . Then the global DG space is We remark that functions in V h are continuous within coarse blocks, but discontinuous across the coarse grid edges in general. The formulation of IPDG method then reads: find u h ∈ V h such that where the bilinear form a DG (u, v) is defined as where and η can take the value of −1, 0 or 1, where we pick η = 1 in this paper that corresponds to the classical Symmetric Interior Penalty Galerkin (SIPG) method. Here γ > 0 is the penalty parameter. Besides, η = 0 is for the Incomplete Interior Penalty Galerkin (IIPG) method and η = −1 is for the Non-symmetric Interior Penalty Galerkin (NIPG) method. In the above bilinear form, [[v]] is the vector jump and [[v]] is the matrix jump to be defined as follows. For coarse grid edges in the interior of the domain, we define where + and − respectively represent the values on two adjacent cells K + and K − sharing the coarse grid edge, and n + and n − are the unit outward normal vector on the boundary of K + and K − . In addition, we define the average of a tensor σ as where σ + , σ − are two tensors defined on K + and K − sharing a common interior coarse edge. For the coarse grid edges lying on the boundary of the domain, we define where n is the normal vector pointing outside of the domain. For each coarse element K i , we define the oversampling domain K i,p by extending K i by p coarse grid layers. In addition, for the i-th coarse grid node, we let ω i be the coarse neighborhood defined by the union of all coarse elements having the i-th node. Figure 2 shows an example of an oversampling domain for p = 1 and a coarse neighborhood.

Construction of multiscale basis functions
In this section, we will present the construction of the multiscale space V cem that will be used in our multiscale method. The main steps include generating multiscale auxiliary functions and then computing multiscale basis functions for global coupling.

Multiscale auxiliary functions
We first show how to construct auxiliary functions. For the j-th coarse block K j , recall that V h (K j ) is the space of piecewise bilinear elements defined on T h in K j . We compute λ j i , ψ j i ∈ R × V h (K j ) by solving the following eigenvalue problem in which is a set of partition of unity functions for the coarse grid, k 2 = λ + 2µ and k 1 = We remark that the choice of the bilinear form (15) is motivated by our convergence analysis. The term |∇χ i | 2 is merely a scaling factor, and is not crucial in our analysis. The most important part in (15) is the inclusion of k 2 . Such choice gives small and contrast dependent eigenvalue for each high contrast channel network in a coarse element. See [14,15] for more details. We will arrange the eigenvalues of the above problem (13) in non-decreasing order. The local auxiliary space W aux (K j ) is formed by the first G j eigenfunctions corresponding to the first G j smallest eigenvalues, that is, The global multiscale auxiliary space W aux is the sum of the local auxiliary spaces as follows, We next define the notion of ψ j i -orthogonality in the space Finally, a projection operator π : V h → W aux is defined by

Multiscale basis functions
In this section, we construct the multiscale trial basis functions. Let ψ j i ∈ W aux be the i-th eigenfunction of the spectral problem (13) for the coarse element K j . We will construct a global basis function φ j i,glo and a local basis function φ j i,ms . The global basis function has support in the whole computational domain. The global basis functions span an approximation space that gives good approximation property. We will show that the global basis functions have an exponentially decay property, which states that the global basis functions have small values outside an oversampling region. The properties suggest the construction of local basis functions, which can be computed efficiently by solving local problems and are more suitable for numerical computations.
We first consider the global version of the basis functions. For each multiscale auxiliary function ψ j i ∈ W aux , there is a multiscale trial basis function φ j i,glo that satisfies the following constraint energy minimization problem where the minimum is sought in the space V h . We note that φ j i,glo is ψ j i -orthogonal. By using the Lagrange multiplier, the above minimization problem can be rewritten as : Then we define the global multiscale space Now, we define our local multiscale trial basis functions. This construction is motivated by an exponential decay property of the basis function φ j i,glo . Therefore, local multiscale basis functions are constructed on some oversampling regions. We introduce the subspace V h (K j,p ) defined by the union of all where the zero Dirichlet boundary condition on ∂K j,p is weakly imposed. Similarly, we use the Lagrange multiplier to write the above problem as: find φ j i,ms ∈ V h (K j,p ) and δ j i,ms ∈ W aux (K j,p ) such that Then we get the localized multiscale trial space, i.e.

Relaxed minimization
We recall that the basis functions in the previous section solve constraint energy minimization problems (23). These constraints can be relaxed. In particular, we consider the following unconstrained energy minimization problems for the constructions of basis functions. Given a function ψ j i ∈ W aux (K j,p ), we compute φ j i,ms ∈ V h (K j,p ) by  This problem is equivalent to the following variational formulation In Figure 3 and Figure 4, we show examples of multiscale basis functions, from which we can see clearly their decay behaviors. Furthermore, the global multiscale basis function φ j i,glo ∈ V h is defined as This problem is equivalent to the following variational formulation Finally, the multiscale solution u off ms is defined as the solution of the following problem, find u off

Convergence Analysis
In this section, the stability and convergence analysis of our method are carried out. We will first give an outline of our analysis. In Lemma 4.1, we will prove the coercivity and continuity of the bilinear form of a DG in our scheme (30). In Lemma 4.3, we will prove that the global basis functions have an exponential decay property. In particular, we will show that the difference between the global basis function and the corresponding local basis function is small if the oversampling region is sufficiently large. In Lemma 4.4, we prove a stable decomposition property. Theorem 4.1 is the main result, which states that the error of approximation in the DG norm is bounded above by O(H) if the size of the oversampling region is proportional to the logarithm of the coarse grid size H and the logarithm of the contrast ratio of the coefficient. In addition, Lemma 4.5 is an analog of the Lemma 4.3 for the relaxed method. Finally, Theorem 4.2 is the main convergence result for the relaxed method. For our analysis, we define the following DG norm The following constants are independent of the mesh size and the number of oversampling layers unless otherwise indicated. We assume that there exist positive constants 0 < c m0 ≤ c m1 such that for a.e.
x ∈ Ω, C (x) is a positive definite matrix with where λ min (C (x)) and λ max (C (x)) are the minimum and the maximum eigenvalues of C (x). We also define the b-norm as · b := b (·, ·). Moreover, for λ j i in (13), we let In the following, we will show the continuity and coercivity of the IPDG bilinear form. For a coarse element K, let n ∂K be the unit outward normal vector on ∂K. Moreover, we define V h (∂K) as the restriction of V h (K) to ∂K. Given u ∈ V h (K), the normal flux σ (u) n ∂K belongs to V h (∂K) and is defined by the following formula: wherev ∈ K is an extension of v to K defined by solving with the boundary conditionv = v on ∂K. Using the results in [7], there is a uniform constant A > 0 such that We will employ the following coercivity and continuity results on the IPDG bilinear form in our study, assuming the penalty value γ is suitably large.
Lemma 4.1 Assume that the penalty parameter γ satisfies γ > Ac 2 m1 c −2 m0 , where the constant A > 0 is defined in (33). Then the bilinear form a DG defined in (6) is continuous and coercive, that is and and by the definition of the DG-norm, we have To estimate the second term in the right hand side of (36), we note that Then by the Cauchy-Schwarz inequality, we obtain We hence proved (35).
(34) can be proved in a similar way. By the definition of a DG and the properties of inequalities, we obtain Notice that The first term S 1 is easily estimated as follows: By the arithmetic-geometric mean inequality, the third term S 3 is also easily estimated as follows: For S 2 , again using the Cauchy-Schwarz inequality and (33), we can estimate it as follows: Then collecting above estimations, if γ ≥ (34) has been proved.
To conclude, we must prove that the global basis functions are indeed localizable. For each coarse block K, B is a bubble function satisfies We take B = Π j χ j where the product is taken over all vertices j on ∂K and define the constant In the following, we will prove a technical result in Lemma 4.2. This is a surjectivity result of the projection operator π. It says that for any function ψ in the auxiliary space W aux , we can find a function φ in the DG space V h such that the DG norm of the function φ is controlled by the norm of the function ψ. Also, the support of the function φ is contained in the support of the function ψ.
Lemma 4.2 Given ψ ∈ W aux , we could find φ ∈ V h , such that We write M = C T C π (1 + c 2 ) , where C T is the maximum number of vertices over all coarse elements.
Proof Let ψ ∈ W aux (K j ) . Then the solution of (23) is found by the following steps. First, we find Note that the existence of (23) is equivalent to find where β is a constant to be determined. Note that ψ is supported in K j . Let B is a bubble function and φ = Bψ, then we obtain Finally, using the (13), we note that This proves that (23) has a unique solution. We see that φ and ψ also satisfy (37) respectively. Hence we are done to prove this lemma.
Next we will prove the convergence of the method.

Convergence Analysis
In this part, we prove the convergence of our DG coupling method using the multiscale basis functions defined in (23). First of all, in Lemma 4.3 we will prove a localization result of the global multiscale basis functions defined in (20). In the result shown in (39), we see that the difference between the local basis function and the global basis function is less than a constant M 1 , which decays exponentially as the number of oversampling size p increases.

Lemma 4.3
We consider an oversampling area K j,p . That is, K j,p is obtained by extending K j by p coarse grid layers. Let ψ j i ∈ W aux be the auxiliary multiscale basis function, φ j i,glo be the solution of (20), and φ j i,ms be the solution of (23). Then we have where M 1 = 4M c 1 + 1 Proof By using the solution of the minimization problems (20) and (23), we can get and where (ψ ) j i , (ψ ) j i,glo ∈ W aux . Substracting the above two equations and restricting w ∈ V h (K j,p ) , we have For ψ j i ∈ W aux , using Lemma 4.2, there exists a function φ j i ∈ V h such that So, we obtain We will now estimate φ j i,ms − φ j i,glo DG , that is, the difference between global multiscale basis functions φ j i,glo and localized multiscale basis functions φ j i,ms . We will observe that the global multiscale basis functions decay and have tiny values outside of a sufficiently large oversampled area. In our proof, we shall employ a cutoff function. Define a function Then supp π χ p+1,p j g ⊂ K j,p+1 \ K j,p .
Remember that π is a projection, restricting ψ j i ≤ 1, then we note that Combining (44) and (45), that is Hence, let l = ψ + χ p+1,p j g, then we have Subsequently, we will estimate (47) in three parts, first estimating 1 − χ p+1,p j g DG and ψ aDG(Kj,p+1\Kj,p) separately, then further estimating the first part's results in the second part, and finally obtaining our final estimate of the value between φ j i,ms and φ j i,glo in the third part.
Then by the definition of g, we obtain Recall the multiscale basis function φ j i,ms satisfies the constraint energy minimization problem a DG φ j i,ms , ς 2 g − γ 1 = 0, therefore, For K ⊂ K j,p \ K j,p−1 , π (g) = 0, then we have then there holds π ς 2 g b(Kj,p\Kj,p−1) ≤ c By the above inequality and the relation of Ω \ K j,p−1 = (Ω \ K j,p ) ∪ (K j,p \ K j,p−1 ), we have Part3. Using (51) in (50), we note that By using (51) in (59), we obtain (60) By the definition of g, we have .
This completes the proof.
We will next prove a stable decomposition property. In particular, we will now estimate Then we obtian Hence, by the definition of ϑ, we note that In addition, by the projection of π and the definition of the cutoff function, we have This completes the proof of the lemma. Now, we are ready to establish our main theorem, which estimates the error between the solution u and the multiscale solution u off ms . To proof the following theorem, the approximate solution u gl ∈ V cem obtained in the global multiscale space V cem is defined by Theorem 4.1 Let u be the solution of (6) and u off ms be the solution of (30). Then we have where u gl is the multiscale solution using corresponding global basis in (66).
H and χ i are bilinear partition of unity, then there exists a constant For the second term in the right hand of (69), recall the Lemma 4.2 and Lemma 4.3, there exists a constant c 3 > 0 such that Next we estimate the first term u − u gl DG in (69). By using the orthogonality property and the eigenvalues in (13), we have and Then we apply Lemma 4.2 to the function This completes the first part of the theorem. Next, we estimate u gl b , Dividing both sides of the above inequality by u gl b , we obtian Substituting (74) into (73), we have We further demonstrate that our method achieves linear convergence in coarse grid size. To obtain our desired goal, using the fact that |∇χ i | = O H −1 and k 1 = Nc i=1 k 2 |∇χ i | 2 , and assuming χ i are bilinear partition of unity, (75) can be rewritten as Finally we only need to proof c

Convergence Analysis of Relaxed Method
In this section, we analyze the convergence of relaxed version of the CEM-GMsFEM. We first consider the basis functions defined in Section 3.3 and estimate the difference between φ j i,ms and φ j i,glo .
Lemma 4.5 Let us consider such an oversampling area K j,p . That is, K j,p is obtained by mapping and extending the oversampling area K j with p coarse grid blocks. Let ψ i j ∈ W aux be the auxiliary multiscale basis function, φ j i,glo be the solution of (28), and φ j i,ms be the solution of (26). Then we have where M 2 = 3 1 + 1 2 c −1 Proof By the above definition, we note that Subtracting the above equations, we obtain Then we have Subsequently, we will estimate the above inequality in three parts, first estimating φ j i,glo − χ p,p−1 separately, then further estimating the first part's results in the second part, and finally obtaining our final estimate of the second part in the third part. , we can get By the definition of cutoff function, note that 1 − χ p+1,p j ≤ 1, then the above inequality can be calculated as . (79) Note that for K ∈ T H , we pick the reflection P := I − π which satisfies . Therefore, we have We will estimate the second term π φ j i,glo − χ p,p−1 . Again by the definition of cutoff function, note that χ p+1,p j ≤ 1, we have Hence, combining (81) and (82), we can get . Then we note that by using the fact that By the definition of DG norm, we have Consequently, we will estimate the global multiscale basis functions φ j i,glo on domain Ω \ K j,p−1 and have Next, since Then we note that .
(85) Finally, we obtain . (87) Using the above inequality recursively, we have This completes the proof.
Finally, we prove the convergence.
Theorem 4.2 Let u be the solution of (6) and u of f ms be the solution of (30). Then we have By the definition of the global multiscale basis functions and the projection π, we have To simplify our inequalities, let b π (u gl ) , ψ l p compose a matrix B, b π φ j i,ms π ψ l p + a DG φ j i,ms , ψ l p compose a matrix A, and (d ij ) compose a matrix D. Then we have Then by the definition of the space V h and based on the method we used in this paper, for ψ = N j=1 Using the given ψ ∈ W aux , there is y ∈ V h such that Let v = y in (91), we can similarly note that Then we have . The rest proof is similarly with Theorem 4.1. This completes the proof.
We remark that Theorem 4.2 gives a better constraint than Theorem 4.1, because the convergence rate in Theorem 4.2 is independent of the constant M . This is supported by the numerical data shown below.

Numerical results
In this section, we present numerical simulation results to show the performance of the proposed method. We consider a complicated high-contrast elastic model whose Young's modulus E (x) is shown in Figure 5, besides, we let The Poisson's ratio v is 0. 25 In order to evaluate the accuracy of the multiscale solution, following errors are adopted: For the convenience of presenting the simulation results, we provide some notations in Table 1. In the standard test, the number of basis function and oversampled element are 6 and 5, respectively, the contrast of the model is 10 4 . The FEM solution and the relaxed CEM-GMsFEM solution with standard parameter settings are shown in Figure 6 and Figure 7. We can see there is obvious distortion in the reference solution due to the strong heterogeneity of the media, however, the coarse-grid solution still succeed in capturing those small details. Actually, it is hard to find any differences between the multiscale solution and the reference solution.   To study the effects of different number of basis functions on the accuracy of the GMsFEM solution, we display the errors with respect to the number of basis functions in Table 2 and Figure 8, we note that in these tests, we let N oc = 5, H = 1/60. It is clear that more number of basis functions will definitely yield more accurate GMsFEM solutions. For example, we can observe that if only 1 basis was used in the Relaxed CEM-GMsFEM, the H 1 error is about 26%, which decays to only 2% if there are 8 bases. We note the L 2 error the Relaxed version is 0.217% if N bf = 8, which is also much smaller than the case of N bf = 1.   It is also important to investigate the influence of oversampling size. To this end, we compare the error as a function of N oc in Table 3 and Figure 9 with fixed N bf = 4 and H = 1/15. It can be found there is sharp decay of the error once N oc reaches 3 for the Relaxed scheme and 5 for the Constraint scheme. Once N oc is larger than a threshold, accuracy improvements are no longer significant.  We also examine the robustness of the proposed method against the contrast of the model. We fix H and N bf and vary the ratio of the maximum value of E(x) over its minimum value and N oc. The simulation results are reported in Table 4. We find that larger N oc leads to more robust coarse-grid simulations, the Relaxed formulation is more robust than the Constraint version. Finally, we study the convergence behavior of the CEM-GMsFEM solution with respect to the coarsegrid size H. We set the number of oversampling layers be N oc ≈ 4 log (H) / log (1/7.5) and N bf = 4 to form the basis spaces. The convergence results are presented in Table 5. In Figure 10, we can see the H 1 error of CEM-GMsFEM solutions almost linearly converge against the coarse-mesh size H as predicted by the analysis.  In summary, the relaxed CEM-GMsFEM is more accurate and more robust than the Constraint CEM-GMsFEM. Besides, more basis functions and larger local problem size lead to more accurate coarse-grid solutions.

Conclusion
In this paper, we introduce the CEM-GMsFEM with the discontinuous Galerkin form, which is a local multiscale model reduction method in the discontinuous Galerkin framework and can be used to solve the elasticity problem of high contrast. Firstly, we use the main mode from the local spectrum problem to construct a multiscale basis function to locally represent the solution space in each coarse block. Then, the multiscale basis function in the coarse oversampling region is solved by constraining the energy minimization problem. After proving the convergence of the CEM-GMsFEM and the Relaxed CEM-GMsFEM in the DG form, we also give the numerical results of a simple elasticity equation to further illustrate the method's accuracy.