A Least Squares Method for Linear Elasticity using A Patch Reconstructed Space

TL;DR

This paper introduces a simple, robust, and efficient discontinuous least squares finite element method for linear elasticity, utilizing patch reconstruction with one unknown per element to achieve optimal convergence.

Contribution

The paper presents a novel discontinuous least squares finite element method for linear elasticity using patch reconstruction, simplifying implementation and improving efficiency.

Findings

Achieves optimal convergence order under energy norm.

Demonstrates robustness and simplicity in implementation.

Numerical results verify theoretical error estimates.

Abstract

We propose a discontinuous least squares finite element method for solving the linear elasticity. The approximation space is obtained by patch reconstruction with only one unknown per element. We apply the L 2 norm least squares principle to the stress-displacement formulation based on discontinuous approximation with normal continuity across the interior faces. The optimal convergence order under the energy norm is attained. Numerical results of linear elasticity are presented to verify the error estimates. In addition to enjoying the advantages of discontinuous Galerkin method, we illustrate the great simplicity in implementation, the robustness and the improved efficiency of our method.

Tables

Table 1. Table 1 . # ​ S ​ ( K ) # 𝑆 𝐾 \#S(K) for d = 2 𝑑 2 d=2 .

$m$	1	2	3	4	5
$S\left(K\right)$	4	8	13	19	26

Table 2. Table 2 . The convergence rates of J h ​ ( 𝝈 h , 𝒖 h ) 1 2 subscript 𝐽 ℎ superscript subscript 𝝈 ℎ subscript 𝒖 ℎ 1 2 J_{h}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})^{\frac{1}{2}} , ‖ 𝝈 − 𝝈 h ‖ L 2 ​ ( Ω ) subscript norm 𝝈 subscript 𝝈 ℎ superscript 𝐿 2 Ω \|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)} and ‖ 𝒖 − 𝒖 h ‖ L 2 ​ ( Ω ) subscript norm 𝒖 subscript 𝒖 ℎ superscript 𝐿 2 Ω \|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{L^{2}(\Omega)} with m = 1 𝑚 1 m=1 .

mesh level	1/10	1/20	1/40	1/80	1/160	1/320	1/640
$J_h{({𝝈}_h,{𝒖}_h)}^{\frac{1}{2}}$	1.265e+1	6.444e-0	3.289e-0	1.655e-0	8.349e-1	4.197e-1	2.105e-1
order	-	0.97	0.97	0.99	0.99	0.99	1.00
${\Vert 𝝈-{𝝈}_h\Vert }_{L^2\left(\Omega \right)}$	3.916e-0	2.943e-0	2.121e-0	1.496e-0	8.417e-1	3.971e-1	1.873e-1
order	-	0.41	0.47	0.50	0.83	1.08	1.08
${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	1.017e-1	8.073e-2	5.461e-2	3.371e-2	1.920e-2	7.781e-3	2.679e-3
order	-	0.33	0.56	0.69	0.81	1.30	1.53

Table 3. Table 3 . J h ​ ( 𝝈 h , 𝒖 h ) 1 2 subscript 𝐽 ℎ superscript subscript 𝝈 ℎ subscript 𝒖 ℎ 1 2 J_{h}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})^{\frac{1}{2}} , ‖ 𝝈 − 𝝈 h ‖ L 2 ​ ( Ω ) subscript norm 𝝈 subscript 𝝈 ℎ superscript 𝐿 2 Ω \|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)} and ‖ 𝒖 − 𝒖 h ‖ L 2 ​ ( Ω ) subscript norm 𝒖 subscript 𝒖 ℎ superscript 𝐿 2 Ω \|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{L^{2}(\Omega)} for different values of λ 𝜆 \lambda with m = 2 𝑚 2 m=2 .

	mesh level	1/10	1/20	1/40	1/80	1/160
$\lambda =5$	$J_h{({𝝈}_h,{𝒖}_h)}^{\frac{1}{2}}$	1.716e-0	4.183e-1	1.011e-1	2.426e-2	5.935e-3
	${\Vert 𝝈-{𝝈}_h\Vert }_{L^2\left(\Omega \right)}$	1.233e-0	2.505e-1	5.289e-2	1.171e-2	2.623e-3
	${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	1.030e-1	2.645e-2	5.869e-3	1.380e-3	3.371e-4
$\lambda =1000$	$J_h{({𝝈}_h,{𝒖}_h)}^{\frac{1}{2}}$	1.715e-0	4.170e-1	1.011e-1	2.426e-2	5.936e-3
	${\Vert 𝝈-{𝝈}_h\Vert }_{L^2\left(\Omega \right)}$	1.192e-0	2.935e-1	6.845e-2	1.389e-2	3.025e-3
	${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	1.143e-1	2.991e-2	6.317e-3	1.439e-3	3.487e-4
$\lambda =20000$	$J_h{({𝝈}_h,{𝒖}_h)}^{\frac{1}{2}}$	1.782e-0	4.203e-1	9.789e-2	2.358e-2	5.799e-3
	${\Vert 𝝈-{𝝈}_h\Vert }_{L^2\left(\Omega \right)}$	1.123e-0	2.781e-1	6.421e-2	1.356e-2	2.972e-3
	${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	1.083e-1	2.946e-2	6.192e-3	1.391e-3	3.366e-4

Table 4. Table 4 . J h ​ ( 𝝈 h , 𝒖 h ) 1 2 subscript 𝐽 ℎ superscript subscript 𝝈 ℎ subscript 𝒖 ℎ 1 2 J_{h}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})^{\frac{1}{2}} , ‖ 𝝈 − 𝝈 h ‖ L 2 ​ ( Ω ) subscript norm 𝝈 subscript 𝝈 ℎ superscript 𝐿 2 Ω \|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)} and ‖ 𝒖 − 𝒖 h ‖ L 2 ​ ( Ω ) subscript norm 𝒖 subscript 𝒖 ℎ superscript 𝐿 2 Ω \|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{L^{2}(\Omega)} for different values of λ 𝜆 \lambda with m = 3 𝑚 3 m=3 .

	mesh level	1/10	1/20	1/40	1/80	1/160
$\lambda =5$	$J_h{({𝝈}_h,{𝒖}_h)}^{\frac{1}{2}}$	6.215e-1	7.762e-2	1.022e-2	1.336e-3	1.729e-4
	${\Vert 𝝈-{𝝈}_h\Vert }_{L^2\left(\Omega \right)}$	4.866e-1	4.313e-2	4.789e-3	5.033e-4	5.193e-5
	${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	1.118e-2	6.008e-4	3.883e-5	2.216e-6	1.420e-7
$\lambda =1000$	$J_h{({𝝈}_h,{𝒖}_h)}^{\frac{1}{2}}$	6.321e-1	7.996e-2	1.048e-2	1.363e-3	1.750e-4
	${\Vert 𝝈-{𝝈}_h\Vert }_{L^2\left(\Omega \right)}$	4.592e-1	4.999e-2	5.412e-3	5.521e-4	6.013e-5
	${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	1.353e-2	6.991e-4	4.638e-5	2.612e-6	1.610e-7
$\lambda =20000$	$J_h{({𝝈}_h,{𝒖}_h)}^{\frac{1}{2}}$	6.198e-1	7.723e-2	1.019e-2	1.333e-3	1.718e-4
	${\Vert 𝝈-{𝝈}_h\Vert }_{L^2\left(\Omega \right)}$	4.761e-1	5.212e-2	5.881e-3	6.122e-4	6.321e-5
	${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	1.278e-2	6.431e-4	4.292e-5	2.598e-6	1.523e-7

Table 5. Table 5 . # ​ S ​ ( K ) # 𝑆 𝐾 \#S(K) for d = 3 𝑑 3 d=3 .

$m$	1	2	3	4
$S\left(K\right)$	9	17	35	57

Table 6. Table 6 . To achieve the same accuracy, the ratio of the DOFs involved in our method to the DOFs involved in continuous finite element method.

	$m=2$	$m=3$	$m=4$
$d=2$	$96.2\%$	$86.6\%$	$68.1\%$
$d=3$	$73.8\%$	$65.6\%$	$53.8\%$

Equations

$$\Vert v{\Vert }_{L^2(\partial K)}^2\le C\left(h_K^{-1}\Vert v{\Vert }_{L^2(K)}^2+h_K\Vert \nabla v{\Vert }_{L^2(K)}^2\right),\forall v\in H^1(K);$$

$$\Vert \nabla v{\Vert }_{L^2(K)}\le C{h}_K^{-1}\Vert v{\Vert }_{L^2(K)},v\in {\mathbb{P}}_m(K),$$

$$S_t(K)=\bigcup _{\begin{array}{rl}\tilde{K}\in {\mathcal{T}}_h,& \hat{K}\in S_{t-1}(K)\\ \tilde{K}\cap \hat{K}& =e\in {\mathcal{E}}_h\end{array}}\tilde{K},t=1,2,\cdots$$

$${\mathcal{I}}_K\triangleq \left\{{\mathbf{x}}_{\tilde{K}}|\forall \tilde{K}\in S(K)\right\}.$$

$${\mathcal{R}}_{K}g={\mathrm{argmin}}_{p\in {\mathbb{P}}_m(S(K))}$$

$$p({\mathbf{x}}_K)=g({\mathbf{x}}_K).$$

$$p{|}_{{\mathcal{I}}_K}=0\text{implies}p{|}_{S(K)}\equiv 0.$$

$$(\mathcal{R}g){|}_K=({\mathcal{R}}_{K}g){|}_K,\text{for }K\in {\mathcal{T}}_h.$$

$$w_K(\mathbf{x})=\{\begin{array}{ll}1,& \mathbf{x}={\mathbf{x}}_K,\\ 0,& \mathbf{x}\text{ outside }K.\end{array}$$

$$\mathcal{R}g=\sum _{K\in {\mathcal{T}}_h}g(\mathbf{x}){\lambda }_K(\mathbf{x}),\forall g\in C^0(\Omega ).$$

$$\Lambda (m,S(K))=\underset{p\in {\mathbb{P}}_m(S(K))}{\max }\frac{{\max }_{\mathbf{x}\in S(K)}|p(\mathbf{x})|}{{\max }_{\mathbf{x}\in {\mathcal{I}}_K}|p(\mathbf{x})|}.$$

$$\Vert {\mathcal{R}}_{m}g{\Vert }_{L^{\infty }(K)}\le \Lambda (m,S(K))\sqrt{\#S(K)}\underset{\mathbf{x}\in {\mathcal{I}}_K}{\max }|g(\mathbf{x})|,$$

$$\Vert g-{\mathcal{R}}_{K}g{\Vert }_{L^{\infty }(K)}\le {\Lambda }_m\underset{p\in {\mathbb{P}}_m(S(K))}{\inf }\Vert g-p{\Vert }_{L^{\infty }(S(K))},$$

$$\Vert g-\mathcal{R}g{\Vert }_{H^q(K)}$$

$$\Vert D^q(g-\mathcal{R}g){\Vert }_{L^2(\partial K)}$$

$$\mathcal{A}\mathbf{\sigma }-\mathbf{\epsilon }(u)$$

$$\nabla \cdot \mathbf{\sigma }+\mathbf{f}$$

$$\mathbf{u}$$

$$\mathrm{n}\cdot \mathbf{\sigma }$$

$${\epsilon }_{i,j}(\mathbf{u})=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$$

$$\mathcal{A}\mathbf{\tau }\triangleq \frac{1}{2\mu }\left(\mathbf{\tau }-\frac{\lambda }{d\lambda +2\mu }(\mathrm{tr}\mathbf{\tau })\mathbf{I}\right),\forall \mathbf{\tau }\in {\mathbb{R}}^{d\times d},$$

$$H_D^1(\Omega )=\left\{v\in H^1(\Omega )|v=0,\text{on}{\Gamma }_D\right\},H_N^1(\Omega )=\left\{v\in H^1(\Omega )|v=0,\text{on}{\Gamma }_N\right\}.$$

$$\Vert \varphi {\Vert }_{-1,D}=\underset{0\ne \psi \in H_D^1(\Omega )}{\sup }\frac{(\varphi ,\psi )}{\Vert \psi {\Vert }_{H^1(\Omega )}}.$$

$$L^2(E{)}^{\mathbb{S},d\times d}$$

$$H^s(E{)}^{\mathbb{S},d\times d}$$

$$H^s({\mathcal{T}}_h{)}^{\mathbb{S},d\times d}$$

$$\{\mathbf{v}\}=\frac{1}{2}\left(\mathbf{v}{|}_{K^{+}}+\mathbf{v}{|}_{K^{-}}\right),[[\mathbf{v}]]=\{\begin{array}{ll}\mathbf{v}{|}_{K^{+}}\otimes {\mathrm{n}}^{+}+\mathbf{v}{|}_{K^{-}}\otimes {\mathrm{n}}^{-},& \text{for vector }\mathbf{v},\\ \mathbf{v}{|}_{K^{+}}\cdot {\mathrm{n}}^{+}+\mathbf{v}{|}_{K^{-}}\cdot {\mathrm{n}}^{-},& \text{for tensor }\mathbf{v},\end{array}$$

$$\{\mathbf{v}\}=\mathbf{v},[[\mathbf{v}]]=\{\begin{array}{ll}\mathbf{v}{|}_K\otimes \mathrm{n},& \text{for vector }\mathbf{v},\\ \mathbf{v}{|}_K\cdot \mathrm{n},& \text{for tensor }\mathbf{v},\end{array}$$

$$J_h(\mathbf{\sigma },\mathbf{u})\triangleq$$

$$+\sum _{e\in {\mathcal{E}}_h^i}\left(\frac{1}{h_e}\Vert [[\mathbf{u}]]{\Vert }_{L^2(e)}^2+\frac{1}{h_e}\Vert [[\mathbf{\sigma }]]{\Vert }_{L^2(e)}^2\right)$$

$$+\sum _{e\in {\Gamma }_D}\frac{1}{h_e}\Vert \mathbf{u}-\mathbf{g}{\Vert }_{L^2(e)}^2+\sum _{e\in {\Gamma }_N}\frac{1}{h_e}\Vert \mathrm{n}\cdot \mathbf{\sigma }-\mathbf{h}{\Vert }_{L^2(e)}^2.$$

$${\mathrm{V}}_h=U_h^d,{\mathbf{\Sigma }}_h=\left\{\mathbf{\tau }\in U_h^{d\times d}|\mathbf{\tau }={\mathbf{\tau }}^T\text{in}\forall K\in {\mathcal{T}}_h\right\}.$$

$$J_h({\mathbf{\sigma }}_h,{\mathbf{u}}_h)=\underset{({\mathbf{\tau }}_h,{\mathbf{v}}_h)\in {\mathbf{\Sigma }}_h\times {\mathrm{V}}_h}{\inf }{J}_h({\mathbf{\tau }}_h,{\mathbf{v}}_h).$$

$$a_h({\mathbf{\sigma }}_h,{\mathbf{u}}_h;{\mathbf{\tau }}_h,{\mathbf{v}}_h)=l_h({\mathbf{\tau }}_h,{\mathbf{v}}_h),\forall ({\mathbf{\tau }}_h,{\mathbf{v}}_h)\in {\mathbf{\Sigma }}_h\times {\mathrm{V}}_h,$$

$$a_h({\mathbf{\sigma }}_h,{\mathbf{u}}_h;{\mathbf{\tau }}_h,{\mathbf{v}}_h)$$

Full text

keywords: Linear elasticity, Least squares method, Patch reconstruction, Discontinuous Galerkin method.

MSC2010: 65N30

1. Introduction

The stress-displacement formulation is the first-order system of the linear elastic problem providing a relation between the strain and the equilibrium equation. Compared to the pure displacement formulation, the stress-displacement formulation is more attractive when considering the nearly incompressible case. For many practical applications, the stress field is the quantity of particular interests. Therefore, there are many efforts devoted to the mixed finite element method based on the weak form of the stress-displacement formulation. We refer to [1, 5, 3, 15, 17] for conforming mixed elements and [16, 4, 24] for non-conforming elements and the references therein. The main challenge of the mixed finite element method is the construction of proper finite element spaces since the requirement of the stable combination of approximation spaces and the symmetric constraint of the stress tensor. It has been a long-standing open problem until Hu and Zhang’s recently significant progress [15, 17]. Their work requires a subtle structure on the geometry of the element to construct the finite element space particularly for the linear elastic problem.

To solve the linear elastic problem using common finite element spaces, the least squares finite element methods have been investigated in a sequence of papers [10, 11, 9, 7, 22]. The least squares finite element method is a sophisticated technique in numerical partial differential equations, and we refer to the survey paper [6] and the references therein. The least squares method based on a discrete minus one inner product is proposed in [7]. Cai and his coworkers developed the least squares finite element methods based on the $L^{2}$ norm residual for solving the stress-displacement system [10, 11, 9]. One of their advantages over the usual mixed finite element method is the selection of the approximation spaces is not subject to the stability condition.

In this paper, a new discontinuous least squares finite element method is proposed based on the stress-displacement formulation. The novel point is the new approximation space which is obtained by patch reconstruction with one unknown per element [20, 21, 19]. The new space could be regarded as a subspace of the common space used in discontinuous Galerkin finite element method. We follow the idea in [11] to define an $L^{2}$ norm least squares functional based on discontinuous approximation spaces and we derive the optimal convergence order under the energy norm. By a series of numerical examples, the error estimates are verified. As in [11], the $L_{2}$ error can only be proved sub-optimal though numerical results show an optimal convergence for odd orders in approximation to the displacement. The implementation of our method is very convenient that the coding for different polynomial degrees and meshes with elements in different geometry are most reusable. We present an example on polygonal mesh to show such advantages of our method. The main steps of the method are detailed in Appendix, which may be helpful to implement for any high order accuracy and any polygonal mesh in an easy manner. The numerical results with very large Lamé parameter are presented to exhibit the robustness in the incompressible limit. A remarkable advantage of this new space is its great efficiency. We make an efficiency comparison between our method and the method in [11] using continuous finite element space. It is clear that our method uses fewer degrees of freedom than continuous approximation to achieve the same accuracy, and for higher order approximation, the advantage in efficiency of our method is more significant.

The rest of this paper is organized as follows. In Section 2, we introduce the reconstruction operator and its corresponding approximation space, and we also present the basic properties of the approximation space. In Section 3, we propose the discontinuous least squares method for the stress-displacement formulation and we derive the error estimate in energy norm. In Section 4, we present a series of the numerical examples to verify the convergence of our method, and we also make a comparison to illustrate the efficiency of the proposed method.

2. Approximation Space

Let $\Omega\subset\mathbb{R}^{d}(d=2,3)$ be a bounded convex domain with smooth boundary $\partial\Omega$. We denote by $\mathcal{T}_{h}$ a collection of polygonal (polyhedral) elements which partition the domain $\Omega$. We define all interior faces of $\mathcal{T}_{h}$ as $\mathcal{E}_{h}^{i}$ and denote by $\mathcal{E}_{h}^{b}$ the set of faces lying on $\partial\Omega$. We then let $\mathcal{E}_{h}=\mathcal{E}_{h}^{i}\cup\mathcal{E}_{h}^{b}$ the set of all faces. Let $h_{K}$ be the diameter of the element $K$ and $h_{e}$ be the size of the face $e$ and we denote $h=h_{\max}=\max_{K\in\mathcal{T}_{h}}h_{K}$ as the mesh size. We assume that $\mathcal{T}_{h}$ is a shape-regular partition of $\Omega$ in the sense of that: there exist

• •

two positive numbers $N$ and $\sigma$ that are independent of the mesh size $h$;

• •

a compatible subdivision $\widetilde{\mathcal{T}}_{h}$ consisting of shape-regular simplexes;

such that

• •

any polygonal (polyhedral) element $K\in\mathcal{T}_{h}$ admits a decomposition $\widetilde{\mathcal{T}}_{h|K}$ into less than $N$ shape-regular simplexes;

• •

the shape-regularity of $\widetilde{K}\in\widetilde{\mathcal{T}}_{h}$ reads [12]: the ratio $h_{\widetilde{K}}/\rho_{\widetilde{K}}$ is bounded by $\sigma$ where $\rho_{\widetilde{K}}$ denotes the radius of the largest ball inscribed in $\widetilde{K}$.

The above regularity assumptions, which are common in finite difference scheme [2] and in DG framework [19], could bring many useful consequences:

• M1

there exists a positive constant $\sigma_{v}$ that is independent of $h$ such that $\sigma_{v}\rho_{K}\leq\rho_{e}$ for any element $K\in\mathcal{T}_{h}$ and any face $e\subset\partial K$;

• M2

[trace inequality] there exists a constant $C$ that is independent $h$ such that

[TABLE]

• M3

[inverse inequality] there exists a constant $C$ that is independent $h$ such that

[TABLE]

where $\mathbb{P}_{m}(\cdot)$ is the polynomial space of degree $\leq m$.

We then define a reconstruction operator with the given partition $\mathcal{T}_{h}$ as follows. First, in every element $K$ we assign a point $\boldsymbol{x}_{K}$ as its corresponding collocation point. The choice of $\boldsymbol{x}_{K}$ could be very flexible, particularly in this paper $\boldsymbol{x}_{K}$ is specified as the barycenter of the element $K$. Second, for each $K\in\mathcal{T}_{h}$ we would construct an element patch $S(K)$ which contains $K$ itself and some elements around $K$. To be specific, for element $K$, a threshold value $\#S(K)$ is given to control the cardinality of $S(K)$, and we construct $S(K)$ recursively. Let $S_{0}(K)=\left\{K\right\}$ and we define $S_{t}(K)$ as follow:

[TABLE]

We end the recursion if the cardinality of $S_{t}(K)$ is greater than $\#S(K)$. Then, we calculate all distances between the collocation points of every element in $S_{t}(K)$ and point $\boldsymbol{x}_{K}$. We select the $\#S(K)$ smallest values and collect the corresponding elements to form the patch $S(K)$. Obviously the cardinality of $S(K)$ is just $\#S(K)$. In Appendix A, we show the details of the algorithm of the construction of the element patch.

Further, we denote by $\mathcal{I}_{K}$ the set containing all collocation points correspond to the elements in $S(K)$:

[TABLE]

Then, for any function $g\in C^{0}(\Omega)$ and element $K\in\mathcal{T}_{h}$ we seek a polynomial of degree $m$ defined on $S(K)$ by solving the following least squares problem:

[TABLE]

The geometrical positions of the points in $\mathcal{I}_{K}$ totally decide the existence and uniqueness of the solution to (3). We follow [21] to make the following assumption:

Assumption 1**.**

For any element $K\in\mathcal{T}_{h}$ and $p\in\mathbb{P}_{m}(S(K))$, one has that

[TABLE]

The assumption in fact excludes the case that all the points in $\mathcal{I}_{K}$ lie on an algebraic curve and demands that the number $\#S(K)$ should be greater than $\text{dim}(\mathbb{P}_{m})$.

It must be notable that the solution to (3) has a linear dependence on the function $g$, which enables us to define a linear reconstruction operator $\mathcal{R}$ for $g$:

[TABLE]

With $\mathcal{R}$, the function $g\in C^{0}(\Omega)$ is mapped into a piecewise polynomial function of degree $m$ on $\mathcal{T}_{h}$. We denote by $U_{h}$ the image of the operator $\mathcal{R}$. Further, we define $w_{K}(\boldsymbol{x})\in C^{0}(\Omega)$ as

[TABLE]

It is clear that $U_{h}=\text{span}\left\{\lambda_{K}\ |\ \lambda_{K}=\mathcal{R}w_{K}\right\}$ and one could explicitly write the reconstruction operator $\mathcal{R}$ as

[TABLE]

In Appendix B, we present an example of linear reconstruction to illustrate the implementation of solving the least squares problem (3).

Then, we would investigate the approximation property of the operator $\mathcal{R}$. We define the constant $\Lambda(m,S(K))$ for all elements as

[TABLE]

With $\Lambda(m,S(K))$, we could state the following estimates.

Lemma 1**.**

For any function $g\in C^{0}(\Omega)$ and element $K\in\mathcal{T}_{h}$, the following inequalities hold true:

[TABLE]

and

[TABLE]

where $\Lambda_{m}\triangleq\max_{K\in\mathcal{T}_{h}}\left\{1+\Lambda(m,S(K))\sqrt{\#S(K)}\right\}$.

Proof.

The proof could be found in [21, Theorem 3.3]. ∎

Under some mild and practical conditions on element patch $S(K)$, $\Lambda_{m}$ could be bounded uniformly which plays a vital role in the convergence estimate. We refer to [21, 20] for these conditions and more detailed discussion about the uniform upper bound. One of the conditions we shall note is that the number $\#S(K)$ shall be far greater than $\text{dim}(\mathbb{P}_{m})$. In Section 4, we list the values of $\#S(K)$ with different $m$ for the numerical tests.

As a direct result of Lemma 1, we could state the following approximation properties of the operator $\mathcal{R}$.

Theorem 1**.**

Let $g\in H^{m+1}(\Omega)$, there exist constants $C$ that are independent of $h$ such that

[TABLE]

Proof.

The proof directly follows from [20, Lemma 4] and [20, Assumption A]. ∎

3. Discontinuous Least Squares Finite Element Method

The problem concerned in this paper is the first-order system formulation of the linear elasticity: seek the stress $\boldsymbol{\sigma}=(\sigma_{ij})_{d\times d}$ and the displacement $\boldsymbol{u}=(u_{1},\dots,u_{d})^{T}$ such that

[TABLE]

where $\boldsymbol{f}$ is a given body force and $\boldsymbol{g},\boldsymbol{h}$ are the boundary conditions. $\Gamma_{D}$ and $\Gamma_{N}$ are two disjoint parts of the boundary $\partial\Omega$ such that $\bar{\Gamma}_{D}\cup\bar{\Gamma}_{N}=\partial\Omega$. For simplicity, $\Gamma_{D}$ is assumed to be non-empty. We denote by $\boldsymbol{\varepsilon}(\boldsymbol{u})=(\varepsilon_{i,j}(\boldsymbol{u}))_{d\times d}$ the symmetric strain tensor:

[TABLE]

The constitutive law with Lamé parameters $\lambda,\mu>0$ is expressed by the linear operator $\mathcal{A}:\mathbb{R}^{d\times d}\rightarrow\mathbb{R}^{d\times d}$:

[TABLE]

where $\boldsymbol{I}$ is the identity operator and $\mathrm{tr}(\cdot)$ denotes the standard trace operator.

Hereafter, let us note that $C$ and $C$ with a subscript that are generic constants that may differ from line to line but are independent of the mesh size $h$ and the parameter $\lambda$, and we will use the standard notation and definition for the spaces $L^{2}(E),L^{2}(E)^{d}$, $L^{2}(E)^{d\times d}$, $H^{s}(E)$, $H^{s}(E)^{d},H^{s}(E)^{d\times d}$ with $s\geq 0$ and $E$ a bounded domain, and their associated inner products and norms. Let

[TABLE]

We denote by $H_{D}^{-1}(\Omega)$ the dual space of $H_{D}^{1}(\Omega)$ with the norm

[TABLE]

Moreover we would use, for the partition $\mathcal{T}_{h}$, the standard broken Sobolev spaces $H^{s}(\mathcal{T}_{h})$, $H^{s}(\mathcal{T}_{h})^{d}$, $H^{s}(\mathcal{T}_{h})^{d\times d}$ with $s\geq 0$ and its corresponding broken norms. For the tensor spaces $L^{2}(E)^{d\times d}$, $H^{s}(E)^{d\times d}$, $H^{s}(\mathcal{T}_{h})^{d\times d}$, we define their corresponding symmetric spaces as

[TABLE]

Then, we introduce the standard trace operators that are commonly used in DG framework. Let $\boldsymbol{v}$ be a vector- or tensor-valued function and $e$ be an interior face shared by elements $K^{+}$ and $K^{-}$ with the unit outward norm $\boldsymbol{\mathrm{n}}^{+}$ and $\boldsymbol{\mathrm{n}}^{-}$ corresponding to $\partial K^{+}$ and $\partial K^{-}$, respectively. The average operator $\{\cdot\}$ and the jump operator $[[\cdot]]$ are defined as follows:

[TABLE]

and in the case $e\in\mathcal{E}_{h}^{b}$, $\{\cdot\}$ and $[[\cdot]]$ are modified as

[TABLE]

where $\boldsymbol{\mathrm{n}}$ is the unit outward normal on $e$.

Now let us define the following least squares functional for the problem (8):

[TABLE]

We introduce two approximation spaces based on the reconstructed space $U_{h}$: $\boldsymbol{\mathrm{V}}_{h}$ for the displacement $\boldsymbol{u}$ and $\boldsymbol{\Sigma}_{h}$ for the stress $\boldsymbol{\sigma}$ as follows:

[TABLE]

Here we impose the symmetric condition of the stress field in the solution space with a strong sense.

In this paper, the discontinuous least squares finite element method reads: find $(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})\in\boldsymbol{\Sigma}_{h}\times\boldsymbol{\mathrm{V}}_{h}$ such that

[TABLE]

To solve the minimization problem (11), one may write its corresponding variational equation which reads: find $(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})\in\boldsymbol{\Sigma}_{h}\times\boldsymbol{\mathrm{V}}_{h}$ such that

[TABLE]

where the bilinear form $a_{h}(\cdot;\cdot)$ and the linear form $l_{h}(\cdot)$ are defined as

[TABLE]

and

[TABLE]

Below we would concentrate on the uniform continuity and ellipticity of the bilinear form $a_{h}(\cdot,\cdot)$. To do so, we first introduce two energy norms $\|\cdot\|_{\boldsymbol{\mathrm{\Sigma}}}$ and $\|\cdot\|_{\boldsymbol{\mathrm{V}}}$:

[TABLE]

Obviously, $\|\cdot\|_{\boldsymbol{\mathrm{\Sigma}}}$ actually defines a norm on the space $H^{1}(\mathcal{T}_{h})^{d\times d}$. The following lemma ensures $\|\cdot\|_{\boldsymbol{\mathrm{V}}}$ is indeed a norm on the space $H^{1}(\mathcal{T}_{h})^{d}$.

Lemma 2**.**

For any function $\boldsymbol{v}_{h}\in H^{1}(\mathcal{T}_{h})^{d}$, the following Korn’s inequality holds true

[TABLE]

Proof.

The proof could be found in [8]. ∎

Then we state the continuity result of the bilinear form with respect to the norms $\|\cdot\|_{\boldsymbol{\mathrm{\Sigma}}}$ and $\|\cdot\|_{\boldsymbol{\mathrm{V}}}$.

Lemma 3**.**

For the bilinear form $a_{h}(\cdot;\cdot)$, the following estimates holds:

[TABLE]

for any $(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h}),(\boldsymbol{\tau}_{h},\boldsymbol{v}_{h})\in H^{1}(\mathcal{T}_{h})^{d\times d}\times H^{1}(\mathcal{T}_{h})^{d}$.

Proof.

We only need to bound the term $\|\mathcal{A}\boldsymbol{\sigma}_{h}\|$:

[TABLE]

which directly gives us

[TABLE]

Applying the Cauchy-Schwarz inequality to (12) and using (15) could yield the estimate (14), which completes the proof. ∎

In order to prove the coercivity of the bilinear form $a_{h}(\cdot,\cdot)$, we may require the following lemmas.

Lemma 4**.**

For any $\boldsymbol{\tau}\in L^{2}(\Omega)^{d\times d}$, there exists a constant $C$ such that

[TABLE]

Proof.

We split the $\|\boldsymbol{\tau}\|_{L^{2}(\Omega)}$ into two parts:

[TABLE]

Thus, the estimate (16) demands a bound of $\|\mathrm{tr}\boldsymbol{\tau}\|_{L^{2}(\Omega)}^{2}$. Then we follow the idea in [10] to apply the Helmholtz decomposition and here we prove for the case $d=2$. Let $\boldsymbol{q}\in H_{D}^{1}(\Omega)$ be the solution of the problem

[TABLE]

whose weak formulation is that $\boldsymbol{q}$ is the only solution of

[TABLE]

Taking $\boldsymbol{\xi}=\boldsymbol{q}$, together with the Poincare inequality $\|q\|_{H^{1}(\Omega)}\leq C|\nabla\boldsymbol{q}|_{L^{2}(\Omega)}$, directly yields

[TABLE]

Further we apply [13, Corollary 2.1] to obtain that

[TABLE]

For vector-valued function $\boldsymbol{v}=(v_{1},v_{2})$, we let $\nabla\times\boldsymbol{v}=\partial_{x_{1}}v_{2}-\partial_{x_{2}}v_{1}$ and $\nabla^{\perp}$ be the formal adjoint of the curl:

[TABLE]

As $\boldsymbol{\tau}-\mathcal{A}^{-1}\nabla\boldsymbol{q}$ is divergence-free, [13, Theorem 3.1] implies a decomposition that there exists a unique solution $\boldsymbol{\phi}\in H_{N}^{1}(\Omega)^{2}$ of

[TABLE]

such that

[TABLE]

From the definition of $\mathcal{A}$ and combining with the regularity of $\boldsymbol{\phi}$, we observe that

[TABLE]

Applying the trace operator brings us that

[TABLE]

From the decomposition, $\nabla^{\perp}\boldsymbol{\phi}$ is divergence free which satisfies that

[TABLE]

Thus, for any $\boldsymbol{v}\in H_{D}^{1}(\Omega)^{2}$, one concludes that

[TABLE]

where $\boldsymbol{I}_{2\times 1}=(1,1)^{T}$. We again apply the estimate (17) to get that

[TABLE]

Collecting above estimates could directly lead to a bound of $\mathrm{tr}\boldsymbol{\tau}$ in $L^{2}$ norm, which completes the proof in the case $d=2$. Besides, the proof could be extended to the case $d=3$ without any difficulty. ∎

Lemma 5**.**

For any $\boldsymbol{\tau}_{h}\in\boldsymbol{\Sigma}_{h}$, there exists a constant $C$ such that

[TABLE]

Proof.

From the definition (9), we clearly have that

[TABLE]

and

[TABLE]

where the last inequality follows from Cauchy-Schwarz inequality and the trace inequality (1), which completes the proof. ∎

Now we are ready to state that the bilinear form $a_{h}(\cdot;\cdot)$ is coercive with respect to the energy norms.

Lemma 6**.**

For the bilinear form $a_{h}(\cdot;\cdot)$, there exists a constant $C$ such that

[TABLE]

for any $(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})\in\boldsymbol{\Sigma}_{h}\times\boldsymbol{\mathrm{V}}_{h}$.

Proof.

The key point is to prove that

[TABLE]

where $a_{h}^{\frac{1}{2}}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h};\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})$ denotes the square root of $a_{h}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h};\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})$.

By using Cauchy-Schwarz inequality, we could observe that

[TABLE]

We then apply the element-wise integration by parts to get

[TABLE]

By Cauchy-Schwarz inequality, we immediately obtain

[TABLE]

For $e\subset\partial K$, we employ the trace inequality (1) and the inverse inequality (2) to find that

[TABLE]

Then we conclude that

[TABLE]

Analogously, we could get that

[TABLE]

Further, by the triangle inequality and Lemma 13 we have

[TABLE]

For the last term, we directly observe that

[TABLE]

Combining above inequalities could yield a bound that

[TABLE]

We substitute this inequality into (20) and we could know that

[TABLE]

which, together with the fact $\|\nabla\cdot\boldsymbol{\sigma}_{h}\|_{-1,D}\leq Ca_{h}^{\frac{1}{2}}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h};\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})$ and Lemma 18, implies

[TABLE]

From the definition of $\|\cdot\|_{\boldsymbol{\mathrm{\Sigma}}}$, we conclude that

[TABLE]

Again we use Lemma 16 to estimate $\|\boldsymbol{u}_{h}\|_{\boldsymbol{\mathrm{V}}}$:

[TABLE]

Hence,

[TABLE]

which completes the proof. ∎

Since the bilinear form $a_{h}(\cdot;\cdot)$ is bounded and coercive, we have established the existence and uniqueness of the solution to the minimization problem (11). Ultimately, we state a priori error estimate of the method proposed in this section:

Theorem 2**.**

Let $(\boldsymbol{\sigma},\boldsymbol{u})\in H^{m+1}(\Omega)^{\mathbb{S},d\times d}\times H^{m+1}(\Omega)^{d}$ be the solution to the problem (8) and let $(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})$ be the solution to the problem (11), there exists a constant $C$ such that

[TABLE]

Proof.

It directly follows from (12) that for any $(\boldsymbol{\tau}_{h},\boldsymbol{v}_{h})\in\boldsymbol{\Sigma}_{h}\times\boldsymbol{\mathrm{V}}_{h}$, one has that

[TABLE]

For any $(\boldsymbol{\tau}_{h},\boldsymbol{v}_{h})\in\boldsymbol{\Sigma}_{h}\times\boldsymbol{\mathrm{V}}_{h}$, together with (19) and (14), we observe that

[TABLE]

By the triangle inequality, we get that

[TABLE]

We denote by $(\boldsymbol{\sigma}_{I},\boldsymbol{u}_{I})\in\boldsymbol{\Sigma}_{h}\times\boldsymbol{\mathrm{V}}_{h}$ be the interpolants of $(\boldsymbol{\sigma},\boldsymbol{u})$ and we only need to estimate the errors of $(\boldsymbol{\sigma}_{I},\boldsymbol{u}_{I})$ under norms $\|\cdot\|_{\boldsymbol{\mathrm{\Sigma}}}$ and $\|\cdot\|_{\boldsymbol{\mathrm{V}}}$, respectively. Using the trace inequality (1) and the approximation property (7), we arrive at

[TABLE]

Substituting the two estimates into (22) implies (21), which completes the proof. ∎

4. Numerical Results

In this section, we carry out a series of numerical results in two and three dimension to exhibit the accuracy and efficiency of the method proposed in Section 3.

4.1. Convergence order study

We first demonstrate the convergence behavior to examine the theoretical prediction and show the flexibility of the proposed method.

Example 1.

We consider a linear elasticity problem defined on the unit square $\Omega=[0,1]\times[0,1]$. We let $\Gamma_{N}$ be the boundary with $x=1$ and $\Gamma_{D}=\partial\Omega\backslash\Gamma_{N}$. The exact solution (see [14]) is taken as

[TABLE]

and the stress $\boldsymbol{\sigma}$, the source term $\boldsymbol{f}$ and the boundary conditions $\boldsymbol{g}$, $\boldsymbol{h}$ are taken accordingly. We fix $\mu=1$ and $\lambda=5$ to test the accuracy of the proposed method.

We solve this problem on a series of triangular meshes (see Fig. 1) with mesh size $h=1/10,1/20,\cdots,1/160$. We set the cardinality $\#S(K)$ uniformly and a group of reference values of $\#S(K)$ for the case $d=2$ are listed in Tab. 1. The values of functional $J_{h}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})$ and errors in $L^{2}$ norm in the approximation to the exact solution $(\boldsymbol{\sigma},\boldsymbol{u})$ are reported in Fig. 2. For fixed $m$, it is clear that the functional $J_{h}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})$, which is equivalent to the error $\left(\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{\boldsymbol{\mathrm{\Sigma}}}^{2}+\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{\boldsymbol{\mathrm{V}}}^{2}\right)^{\frac{1}{2}}$, converges to zero at the rate $O(h^{m})$ as the mesh size approaches to zero, and the error $\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}$ decreases to zero at the same speed. For odd $m$, the error $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{L^{2}(\Omega)}$ converges to zero optimally and for even $m$, the convergence order reduces to $m$. We note that for the case $m=1$ the convergence rate of the error $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{L^{2}(\Omega)}$ seems less than the expected value. The predicted convergence rates would recur when the mesh size $h$ is small enough (see Tab. 2). We refer to [10] for the possible reason and we mainly consider the case $m\geq 2$ in the rest of this section. In addition, all numerically detected convergence orders are in agreement with the theoretical analysis.

Moreover, we exhibit the robustness of the proposed by solving the problem with $\lambda\rightarrow\infty$. We still take $\boldsymbol{u}(x,y)$ as the exact solution but we select $\lambda=1000,20000$. The polynomial degree $m$ is chosen as $2,3$. The least squares functional and the errors under $L^{2}$ norm are gathered in Tab. 3 and Tab. 4 for decreasing mesh size $h$ and different values of the Lamé parameter $\lambda$. Clearly, the numerical solutions produced by our method converge uniformly as $\lambda$ increases, which confirms the convergence analysis.

Example 2.

In this test, we solve the same problem as in Example 1 but on a sequence of polygonal meshes. The meshes are generated by PolyMesher [23] and contain very general polygonal elements; see Fig. 3. The functional $J_{h}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})^{\frac{1}{2}}$ and the $L^{2}$ norm errors in approximation to $(\boldsymbol{\sigma},\boldsymbol{u})$ are displayed in Fig. 4. Again we observe the values of $J_{h}(\boldsymbol{\sigma}_{h},\boldsymbol{u}_{h})^{\frac{1}{2}}$ and $\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{L^{2}(\Omega)}$ tend to zero at the rate $O(h^{m})$. The convergence order of the error $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{L^{2}(\Omega)}$ is $O(h^{m+1})$ and $O(h^{m})$ for odd and even $m$, respectively. The numerical results validate our theoretical estimates and highlight the flexibility of the proposed method.

Example 3.

In this test, we compute an example in three dimension. We solve the linear elasticity problem on the unit cube $\Omega=[0,1]^{3}$ and we set $\Gamma_{N}$ is the boundary with $y=1$ and $\Gamma_{D}=\partial\backslash\Gamma_{N}$ . Let the exact displacement $u$(see [17]) be

[TABLE]

Then, the load function $\boldsymbol{f}$ and the boundary functions $\boldsymbol{g},\boldsymbol{h}$ are defined accordingly. The uniform cardinality $\#S(K)$ is given in Tab. 5 for the case $d=3$. We solve this example with a series of tetrahedral meshes with mesh size $h=1/3$, $1/6$, $1/12$, $1/24$. The Lamé parameters $\lambda$, $\mu$ in this test are $1$. The least squares functional and the errors in various norms with different $h$ and different $m$ are reported in Fig. 5, which distinctly coincide with the theoretical predicts.

4.2. Efficiency comparison

Now let us make a comparison between our method and the classical continuous least squares method proposed in [11]. According to [18], the number of the degrees of freedom of a specific discrete system could serve as a proper indicator for the scheme’s efficiency. For two dimension, we solve the problem taken from Example 1 by the two methods on a series of triangular meshes, respectively, and for three dimension we employ both methods to solve the problem in Example 3. Here for the continuous least squares method, we adopt the standard Lagrange finite element space of degree $2\leq m\leq 4$ for each component of the symmetric stress. To show the efficiency of our method, we compare the values of $J_{h}(\cdot,\cdot)$ and $\widetilde{J}_{h}(\cdot,\cdot)$, where $\widetilde{J}_{h}(\cdot,\cdot)$ is the least squares functional defined in [11]. The main difference between the two least squares functional is $\widetilde{J}_{h}(\cdot,\cdot)$ contains no jump term defined on the interior faces.

In Fig. 6 and Fig. 7, we plot the values of the least squares functions defined for two methods against the number of the degrees of freedom with $2\leq m\leq 4$. All convergence rates are perfectly consistent with the analysis. Clearly, our method is more efficient than the continuous least squares method. To achieve the same accuracy, fewer degrees of freedom are involved in our method for all $2\leq m\leq 4$, and the advantage of the efficiency of our method becomes more prominent with the increasing of $m$. More specifically, in Tab. 6 we list the ratio between the number of DOFs in our method and the number of DOFs in continuous finite element method when the two methods achieve the same accuracy. The saving of number of DOFs is more remarkable when adopting the high-order approximation.

5. Conclusion

We proposed a new discontinuous least squares method for the linear elasticity problem. The approximation space is reconstructed by solving the local least squares problem. We proved the optimal convergence order in energy norm. We conducted a sequence of numerical results that supported our theoretical results and exhibited the great flexibility, robustness and efficiency of the proposed method.

Acknowledgements

This research is supported by the National Natural Science Foundation of China (Grant No. 91630310, 11421110001, and 11421101) and the Science Challenge Project, No. TZ2016002.

Appendix A Algorithm of constructing element patch

In Algorithm 1 we show the details of the algorithm of constructing the element patch for every element in partition, which is very simple to implement.

Appendix B Solving local least squares problem

In this section, we present an example for solving the local least squares problem (3). We choose the linear reconstruction as an illustration. For element $K_{0}$ (the red element in Fig. 8), we collect $K_{0}$ and its face-neighbouring elements to form the element patch $S(K_{0})=\left\{K_{0},K_{1},K_{2},K_{3}\right\}$, see Fig. 8. We let $\mathcal{I}_{K_{0}}=\left\{\boldsymbol{x}_{0},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\boldsymbol{x}_{3}\right\}$, where $\boldsymbol{x}_{i}$ is the barycenter of $\boldsymbol{x}_{i}$.

Then for $g\in C^{0}(\Omega)$, the least squares problem on $S(K_{0})$ reads:

[TABLE]

Since the constraint of (23), $p(\boldsymbol{x})$ has the form

[TABLE]

where $\boldsymbol{x}=(x,y)$ and $\boldsymbol{x}_{i}=(x_{i},y_{i})$. Then, the problem (23) is equivalent to

[TABLE]

It is easy to get the unique solution to (24):

[TABLE]

where

[TABLE]

Hence,

[TABLE]

where $I_{3\times 1}=(1,1,1)^{T}$. It is noticeable that the matrix $M$ is independent on the function $g$ and contains all information of the function $\lambda_{K_{0}}$, $\lambda_{K_{1}}$, $\lambda_{K_{2}}$, $\lambda_{K_{3}}$ on the element $K_{0}$ according to the expression (4). Then we store the matrix $M$ on all elements to represent the approximation space $U_{h}$. The procedure of this implementation could be adapted to the case of greater $m$ without difficulties.