Open Access Published by De Gruyter March 26, 2022

On the Spectrum of an Operator Associated with Least-Squares Finite Elements for Linear Elasticity

Linda Alzaben , Fleurianne Bertrand and Daniele Boffi

From the journal Computational Methods in Applied Mathematics

https://doi.org/10.1515/cmam-2022-0044

Abstract

In this paper we provide some more details on the numerical analysis and we present some enlightening numerical results related to the spectrum of a finite element least-squares approximation of the linear elasticity formulation introduced recently. We show that, although the formulation is robust in the incompressible limit for the source problem, its spectrum is strongly dependent on the Lamé parameters and on the underlying mesh.

Keywords: Eigenvalue Problem; Linear Elasticity; Least-Squares Finite Elements

MSC 2010: 65N25; 74B05

1 Introduction

In this paper we continue the discussion started in [4, 5] about the spectrum of operators arising from least-squares finite element approximation of partial differential equations. In [4] several least-squares formulations associated with the Laplace problem and in [5] two formulations associated with the linear elasticity problem were considered. Here we continue the analysis of the two-field formulation presented in [5]. The two-field formulation was introduced in [8] for the approximation of the source problem and has the merit of providing a robust discretization also when the system approaches the incompressible limit.

Several methods are available for the computation of the eigenvalues and eigenfunctions in linear elasticity and the aim of this paper is not to compete with them. However, a good knowledge of the properties of the spectrum is useful for several applications. For instance, if a transient problem is approximated by using the two-field approach for the space semi-discretization, then the behavior of the solution depends essentially on the discrete eigenvalues of the operator corresponding to the least-squares model. Hence, we believe that a study of the approximation of the eigenvalues and the eigenfunctions may be interesting in view of a better understanding of the scheme under investigation.

The aim of this paper is two-fold. First, we complete the analysis presented in [5] where only a sketch of the main ideas were indicated. Actually, since the least-squares formulation is not symmetric, the convergence analysis should take into account the dual problem for which we provide a careful description in this paper. Second, we present a bunch of numerical experiments that highlight some properties of the spectrum. A peculiarity of our operators is that the continuous problem has positive and real eigenvalues, while its approximation may have eigenvalues everywhere in the complex plane. Our numerical tests show that for small values of the Lamé constant λ (i.e., when the considered elastic solid is far from being incompressible) the discrete spectrum is generally well-behaving, that is, it is distributed in the right half of the complex plane and has a small imaginary part. On the other hand, when the solid tends to the incompressible limit, as λ increases, the distribution of the discrete spectrum is more spread in the entire complex plane, including eigenvalues with negative real part. We present several examples of this behavior and discuss how it may depend on the chosen mesh sequence. We conclude that, although the two-field formulation has been proved to be robust for the approximation of the source problem corresponding to linear elasticity, the spectrum of the discrete solution operator is more stable when a small value of λ is considered.

It will be interesting in the future to discuss in more detail the asymptotic exactness of the least-squares operator associated to linear elasticity proved in [9], and to see how it is related to eigenvalue problems and to the situation presented in this paper.

In Section 2 we describe the first order system corresponding to linear elasticity and recall its two-field least-squares representation. Section 3 deals with its finite element discretization and Section 4 develops the convergence analysis. Finally, we report our numerical results in Section 5.

2 The Continuous Problem

Consider a polytopal domain Ω in ℝ d ( d = 2 , 3 ) . We partition the boundary of the domain into two open subsets Γ D and Γ N such that ∂ ⁡ Ω = Γ ¯ D ∪ Γ ¯ N and Γ D ∩ Γ N = ∅ with 𝐧 denoting the outward unit vector normal to Γ N . The linear elasticity problem that we are going to consider consists of finding a stress tensor 𝝈 ¯ = ( σ i , j ) d × d and a displacement field 𝐮 = ( u 1 , … , u d ) ⊤ for a given body force 𝐟 = ( f 1 , … , f d ) ⊤ that satisfy the system

(2.1) { 𝝈 ¯ - 𝒞 ⁢ ϵ ¯ ⁢ ( 𝐮 ) = 𝟎 ¯ in ⁢ Ω , 𝐝𝐢𝐯 𝝈 ¯ = - 𝐟 in ⁢ Ω , 𝐮 = 𝟎 on ⁢ Γ D , 𝐧 ⋅ 𝝈 ¯ = 𝟎 on ⁢ Γ N ,

where ϵ ¯ ⁢ ( 𝐯 ) is the symmetric gradient (also known as the strain tensor) given by

ϵ ¯ ⁢ ( 𝐯 ) = 1 2 ⁢ ( ∇ ¯ ⁢ 𝐯 + ( ∇ ¯ ⁢ 𝐯 ) ⊤ ) ,

and 𝒞 is the elasticity tensor (a symmetric operator) for isotropic, homogeneous material defined in terms of the Lamé constants μ and λ as

𝒞 ⁢ ϵ ¯ ⁢ ( 𝐯 ) = λ ⁢ tr ⁡ ( ϵ ¯ ⁢ ( 𝐯 ) ) ⁢ 𝐈 ¯ + 2 ⁢ μ ⁢ ϵ ¯ ⁢ ( 𝐯 ) ,

with 𝐈 ¯ being the identity tensor.

It is well known that there is a relation between stress and strain tensors such that

(2.2) 𝝉 ¯ = 𝒞 ⁢ ϵ ¯ ⁢ ( 𝐯 ) or ϵ ¯ ⁢ ( 𝐯 ) = 𝒜 ⁢ 𝝉 ¯ ,

where 𝒜 is the compliance tensor given by

(2.3) 𝒜 ⁢ 𝝉 ¯ = 1 2 ⁢ μ ⁢ ( 𝝉 ¯ - λ d ⁢ λ + 2 ⁢ μ ⁢ tr ⁡ ( 𝝉 ¯ ) ⁢ 𝐈 ¯ ) .

If λ is finite, then from (2.2) we have 𝒜 = 𝒞 - 1 . However, as λ approaches ∞ the elasticity tensor blows up and the material becomes nearly incompressible or incompressible. For this reason it may be more convenient to use the relation ϵ ¯ ⁢ ( 𝐯 ) = 𝒜 ⁢ 𝝉 ¯ if we are interested in formulations that are stable uniformly in λ. We can then rewrite system (2.1) as: find a symmetric d × d stress tensor 𝝈 ¯ and a displacement vector-field 𝐮 such that

(2.4) { 𝒜 ⁢ 𝝈 ¯ - ϵ ¯ ⁢ ( 𝐮 ) = 𝟎 ¯ in ⁢ Ω , 𝐝𝐢𝐯 𝝈 ¯ = - 𝐟 in ⁢ Ω , 𝐮 = 𝟎 on ⁢ Γ D , 𝐧 ⋅ 𝝈 ¯ = 𝟎 on ⁢ Γ N .

We are interested in the spectrum of the solution operator associated with (2.4), hence, we replace the source term 𝐟 by ω ⁢ 𝐮 , where ω is the eigenvalue associated with 𝐮 . Since our problem is symmetric, we are looking for real eigenvalues ω ∈ ℝ such that for non-vanishing 𝐮 and some 𝝈 ¯ the following set of equations is satisfied:

(2.5) { 𝒜 ⁢ 𝝈 ¯ - ϵ ¯ ⁢ ( 𝐮 ) = 𝟎 ¯ in ⁢ Ω , 𝐝𝐢𝐯 𝝈 ¯ = - ω ⁢ 𝐮 in ⁢ Ω , 𝐮 = 𝟎 on ⁢ Γ D , 𝐧 ⋅ 𝝈 ¯ = 𝟎 on ⁢ Γ N .

The eigenvalue problem (2.5) is compact due to the regularity properties of the solution of system (2.4). Therefore, the eigenvalues are real and positive and form an increasing sequence

0 < ω 1 ≤ ω 2 ≤ ⋯ ≤ ω i ≤ ⋯ with ⁢ lim i → ∞ ⁡ ω i = ∞ ,

repeated according to their multiplicity so that each eigenvalue ω i corresponds to a one-dimensional eigenspace E i .

We rewrite (2.4) following the least-squares principle approach introduced in [8]. As stated by the authors and as a consequence of the above considerations, this approach has the important advantage of automatically stabilizing the stress-displacement system in the incompressible limit when comparing it to other approaches, thus being more robust. The so called two-field formulation (in the unknowns 𝐮 and 𝝈 ¯ ) was achieved by applying the L 2 norm and minimizing the functional

(2.6) ℱ ( 𝝉 ¯ , 𝐯 ; 𝐟 ) = ∥ 𝒜 𝝉 ¯ - ϵ ¯ ( 𝐯 ) ∥ 0 2 + ∥ 𝐝𝐢𝐯 𝝉 ¯ + 𝐟 ∥ 0 2 ,

in 𝑿 ¯ N × H 0 , D 1 ⁢ ( Ω ) d , where

𝑿 ¯ = { 𝐇 ⁢ ( 𝐝𝐢𝐯 ; Ω ) d if ⁢ Γ N ≠ ∅ , { 𝝉 ¯ ∈ 𝐇 ⁢ ( 𝐝𝐢𝐯 ; Ω ) d : ∫ Ω tr ⁡ ( 𝝉 ¯ ) ⁢ d 𝐱 = 0 } if ⁢ Γ N = ∅ ,

and 𝑿 ¯ N being its subspace

𝑿 ¯ N = { 𝝉 ¯ ∈ 𝑿 ¯ : 𝐧 ⋅ 𝝉 ¯ = 𝟎 ⁢ on ⁢ Γ N } .

Applying the same strategy to (2.5) would lead to a non-linear problem. Following the original idea developed in [4] for the Laplace operator, in [5] the spectrum of the operators associated with the least-squares 𝑠𝑜𝑢𝑟𝑐𝑒 formulation was considered.

2.1 The Variational Formulation

The minimization of the functional corresponding to the source problem in (2.6), gives rise to the variational problem: find ( 𝝈 ¯ , 𝐮 ) ∈ 𝐗 ¯ N × H 0 , D 1 ⁢ ( Ω ) d such that

(2.7) { ( 𝒜 ⁢ 𝝈 ¯ , 𝒜 ⁢ 𝝉 ¯ ) + ( 𝐝𝐢𝐯 𝝈 ¯ , 𝐝𝐢𝐯 𝝉 ¯ ) - ( 𝒜 ⁢ 𝝉 ¯ , ϵ ¯ ⁢ ( 𝐮 ) ) = - ( 𝐟 , 𝐝𝐢𝐯 𝝉 ¯ ) for all ⁢ 𝝉 ¯ ∈ 𝑿 ¯ N , - ( 𝒜 ⁢ 𝝈 ¯ , ϵ ¯ ⁢ ( 𝐯 ) ) + ( ϵ ¯ ⁢ ( 𝐮 ) , ϵ ¯ ⁢ ( 𝐯 ) ) = 𝟎 for all ⁢ 𝐯 ∈ H 0 , D 1 ⁢ ( Ω ) d ,

so that the eigenvalue variational problem associated with the two-field formulation reads: find eigensolutions ( ω , 𝐮 ) ∈ ℝ × H 0 , D 1 ⁢ ( Ω ) d with 𝐮 ≠ 0 such that for some 𝝈 ¯ ∈ 𝑿 ¯ N we have

(2.8) { ( 𝒜 ⁢ 𝝈 ¯ , 𝒜 ⁢ 𝝉 ¯ ) + ( 𝐝𝐢𝐯 𝝈 ¯ , 𝐝𝐢𝐯 𝝉 ¯ ) - ( 𝒜 ⁢ 𝝉 ¯ , ϵ ¯ ⁢ ( 𝐮 ) ) = - ω ⁢ ( 𝐮 , 𝐝𝐢𝐯 𝝉 ¯ ) for all ⁢ 𝝉 ¯ ∈ 𝑿 ¯ N , - ( 𝒜 ⁢ 𝝈 ¯ , ϵ ¯ ⁢ ( 𝐯 ) ) + ( ϵ ¯ ⁢ ( 𝐮 ) , ϵ ¯ ⁢ ( 𝐯 ) ) = 𝟎 for all ⁢ 𝐯 ∈ H 0 , D 1 ⁢ ( Ω ) d .

As discussed in [5, 4], problem (2.8) has the non-symmetric structure

(2.9) ( A B ⊤ B C ) ⁢ ( x y ) = ω ⁢ ( 0 D 0 0 ) ⁢ ( x y ) ,

where the operators are associated to the bilinear forms as follows:

(2.10) { A : ( 𝒜 ⁢ 𝝈 ¯ , 𝒜 ⁢ 𝝉 ¯ ) + ( 𝐝𝐢𝐯 𝝈 ¯ , 𝐝𝐢𝐯 𝝉 ¯ ) , B : - ( 𝒜 ⁢ 𝝈 ¯ , ϵ ¯ ⁢ ( 𝐯 ) ) , C : ( ϵ ¯ ⁢ ( 𝐮 ) , ϵ ¯ ⁢ ( 𝐯 ) ) , D : - ( 𝐮 , 𝐝𝐢𝐯 𝝉 ¯ )

such that x and y are associated to 𝝈 ¯ and 𝐮 , respectively, in an abstract setting. Compared to the Laplacian case in [4], we can see that the operators B ⊤ and - D are different from each other, being related to the bilinear forms - ( 𝒜 ⁢ 𝝉 ¯ , ϵ ¯ ⁢ ( 𝐮 ) ) and ( 𝐮 , 𝐝𝐢𝐯 𝝉 ¯ ) . In particular, it is not possible to show that (2.9) corresponds to a symmetric problem. Thus, we might expect complex eigenvalues among the approximation of the solutions of (2.8) even if the eigenvalues of the continuous problem are real. This is different from the structure of the FOSLS of the Poisson equation in [4] which corresponds to a symmetric problem when applying the integration by part to the bilinear form associated to - D (see [4, Section 2.1]).

3 Galerkin Discretization

The discrete variational formulation associated with (2.8) is established by introducing finite-dimensional subspaces Σ h ⊂ 𝑿 ¯ N , U h ⊂ H 0 , D 1 ⁢ ( Ω ) d and by considering discrete variables 𝝈 ¯ h ⊂ Σ h and 𝐮 h ⊂ U h , respectively. Then the Galerkin approximation of the eigenvalue problem is: find ( ω h , 𝐮 h ) ∈ ℂ × U h with a non-vanishing 𝐮 h and some 𝝈 ¯ h ∈ Σ h such that

(3.1) { ( 𝒜 ⁢ 𝝈 ¯ h , 𝒜 ⁢ 𝝉 ¯ ) + ( 𝐝𝐢𝐯 𝝈 ¯ h , 𝐝𝐢𝐯 𝝉 ¯ ) - ( 𝒜 ⁢ 𝝉 ¯ , ϵ ¯ ⁢ ( 𝐮 h ) ) = - ω h ⁢ ( 𝐮 h , 𝐝𝐢𝐯 𝝉 ¯ ) for all ⁢ 𝝉 ¯ ∈ Σ h , - ( 𝒜 ⁢ 𝝈 ¯ h , ϵ ¯ ⁢ ( 𝐯 ) ) + ( ϵ ¯ ⁢ ( 𝐮 h ) , ϵ ¯ ⁢ ( 𝐯 ) ) = 𝟎 for all ⁢ 𝐯 ∈ U h .

This discrete eigenvalue problem has the same structure as in (2.9) with the natural definition of matrices that are associated with the bilinear forms defined in (2.10). Looking at the problem from an algebraic point of view, the solution of this problem satisfies some properties that are discussed in [5, 4]. In general, our discrete eigenvalue problem in (3.1) has the form of a generalized eigenvalue problem

(3.2) ℳ ⁢ x = ω ⁢ 𝒩 ⁢ x ,

where in our framework the matrix ℳ is symmetric and invertible, while 𝒩 is non-symmetric and singular. The aim of this paper is not to discuss how to solve efficiently our problem. We just indicate that a possible way to resolve the singularity of the system is to switch the roles of the two matrices ℳ , 𝒩 and to consider the solutions ( γ , x ) of the problem

𝒩 ⁢ x = γ ⁢ ℳ ⁢ x .

If γ = 0 , we will say that the corresponding eigenvalue of (3.2) is ω = ∞ . The remaining finite eigenmodes are ( ω , x ) where ω = 1 γ .

We report the following proposition from [5] that classifies the eigenvalues of our problem.

Proposition 1.

Consider the matrices associated with the bilinear forms defined in (2.10). Then the generalized eigenvalue problem

(3.3) ( A B ⊤ B C ) ⁢ ( 𝝈 ^ h 𝐮 ^ h ) = ω h ⁢ ( 0 D 0 0 ) ⁢ ( 𝝈 ^ h 𝐮 ^ h )

has three families of eigenvalues, specifically:

ω h = ∞ with multiplicity equal to dim ⁡ ( Σ h ) ,
ω h = ∞ with multiplicity equal to dim ker ( D ) ) if D is not full rank,
ω h = complex eigenvalues counted with their multiplicity equal to rank ⁡ ( D ) .

We are actually interested in the eigenpairs corresponding to a non-zero displacement 𝐮 h ≠ 𝟎 . This can be easily identified by considering the following Schur complement approach. We can rewrite the algebraic linear system of the block matrices in (3.3) as follows:

(3.4) { A ⁢ 𝝈 ^ h + B ⊤ ⁢ 𝐮 ^ h = ω h ⁢ D ⁢ 𝐮 ^ h , B ⁢ 𝝈 ^ h + C ⁢ 𝐮 ^ h = 𝟎 .

Looking at the pair of interest, and using the fact that A is invertible, we can rewrite the first equation in (3.4) as follows:

𝝈 ^ h = A - 1 ⁢ ( ω h ⁢ D - B ⊤ ) ⁢ 𝐮 ^ h .

Now by substituting the value of 𝝈 ^ h in the second equation of system (3.4) and by combining common terms together, we arrive at the Schur complement linked only to the displacement given by

(3.5) ( B ⁢ A - 1 ⁢ B ⊤ - C ) ⁢ 𝐮 ^ h = ω h ⁢ B ⁢ A - 1 ⁢ D ⁢ 𝐮 ^ h .

Remark 1.

In general, one can deduce two Schur complements associated with system (3.4), one linked with the 𝝈 ^ h variable and another to the variable 𝐮 ^ h . That is, we can initially seek for the pair ( ω h , 𝝈 ^ h ) and then obtain the eigenpair ( ω h , 𝐮 ^ h ) (see [2, Section 3.1]). If one naively picks the Schur complement associated with 𝝈 ^ h , artificial displacement (that is, 𝐮 ^ h = 0 ) which is not related to our problem occurs, thus, wrong solution would then be present. Hence, one needs to be careful while selecting the proper Schur complement associated to the variable of interest in order to ensure that we are considering only genuine pairs for the displacement.

4 Convergence Analysis

The convergence analysis of eigenvalue problems has the standard abstract setting presented in [7, 3]. Usually the analysis is divided into two steps: first considering the convergence of eigenmodes (with the absence of superior modes) that is a consequence of the convergence of the discrete solution operator to the continuous one, then discussing the rate of convergence.

To this end, we introduce the solution operator T : L 2 ⁢ ( Ω ) → L 2 ⁢ ( Ω ) associated with formulation (2.8). Given 𝐟 ∈ L 2 ⁢ ( Ω ) d , we define T ⁢ 𝐟 ∈ H 0 1 ⁢ ( Ω ) d being the second component 𝐮 of the solution of (2.7), which solves the following problem for some 𝝈 ¯ ∈ 𝑿 ¯ N :

(4.1) { ( 𝒜 ⁢ 𝝈 ¯ , 𝒜 ⁢ 𝝉 ¯ ) + ( 𝐝𝐢𝐯 𝝈 ¯ , 𝐝𝐢𝐯 𝝉 ¯ ) - ( 𝒜 ⁢ 𝝉 ¯ , ϵ ¯ ⁢ ( T ⁢ 𝐟 ) ) = - ( 𝐟 , 𝐝𝐢𝐯 𝝉 ¯ ) for all ⁢ 𝝉 ¯ ∈ 𝑿 ¯ N , - ( 𝒜 ⁢ 𝝈 ¯ , ϵ ¯ ⁢ ( 𝐯 ) ) + ( ϵ ¯ ⁢ ( T ⁢ 𝐟 ) , ϵ ¯ ⁢ ( 𝐯 ) ) = 𝟎 for all ⁢ 𝐯 ∈ H 0 , D 1 ⁢ ( Ω ) d .

Since the range of the operator T is included in the space H 0 1 ⁢ ( Ω ) d which is compact in L 2 , it follows that T is compact.

The discrete solution operator T h : L 2 ⁢ ( Ω ) → L 2 ⁢ ( Ω ) with 𝐟 ∈ L 2 ⁢ ( Ω ) d is defined as the second component T h ⁢ 𝐟 ∈ U h of the solution of the Galerkin approximation of (2.7) which solves the following problem for some 𝝈 ¯ h ∈ Σ h :

(4.2) { ( 𝒜 ⁢ 𝝈 ¯ h , 𝒜 ⁢ 𝝉 ¯ ) + ( 𝐝𝐢𝐯 𝝈 ¯ h , 𝐝𝐢𝐯 𝝉 ¯ ) - ( 𝒜 ⁢ 𝝉 ¯ , ϵ ¯ ⁢ ( T h ⁢ 𝐟 ) ) = - ( 𝐟 , 𝐝𝐢𝐯 𝝉 ¯ ) for all ⁢ 𝝉 ¯ ∈ Σ h , - ( 𝒜 ⁢ 𝝈 ¯ h , ϵ ¯ ⁢ ( 𝐯 ) ) + ( ϵ ¯ ⁢ ( T h ⁢ 𝐟 ) , ϵ ¯ ⁢ ( 𝐯 ) ) = 𝟎 for all ⁢ 𝐯 ∈ U h .

In [5] it was stated that the uniform convergence of T h to T holds true, so that the discrete eigenvalues converge to the continuous ones without spurious modes. In the next theorem we will give a rigorous proof of the uniform convergence, by giving more details than what was previously published.

We start by recalling the expression of the generic dual problem associated with the variational formulation of the source problem as follows: given 𝐠 ∈ L 2 ⁢ ( Ω ) d find 𝝌 ¯ ∈ 𝑿 ¯ N and 𝐩 ∈ H 0 1 ⁢ ( Ω ) d such that

(4.3) { ( 𝒜 ⁢ 𝝌 ¯ , 𝒜 ⁢ 𝝃 ¯ ) + ( 𝐝𝐢𝐯 𝝌 ¯ , 𝐝𝐢𝐯 𝝃 ¯ ) - ( 𝒜 ⁢ 𝝃 ¯ , ϵ ¯ ⁢ ( 𝐩 ) ) = 𝟎 for all ⁢ 𝝃 ¯ ∈ 𝑿 ¯ N , - ( 𝒜 ⁢ 𝝌 ¯ , ϵ ¯ ⁢ ( 𝐯 ) ) + ( ϵ ¯ ⁢ ( 𝐩 ) , ϵ ¯ ⁢ ( 𝐯 ) ) = ( 𝐠 , 𝐯 ) for all ⁢ 𝐯 ∈ H 0 , D 1 ⁢ ( Ω ) d .

Theorem 2.

Let s > 1 2 be such that u ∈ H 1 + s ⁢ ( Ω ) d , where u is the second component of the solution to (2.7), and such that the following stability of the dual problem holds true:

(4.4) ∥ 𝐩 ∥ s + 1 + ∥ 𝝌 ¯ ∥ s + ∥ 𝐝𝐢𝐯 𝝌 ¯ ∥ s + 1 ≤ C ⁢ ∥ 𝐠 ∥ 0 .

Let u h ∈ U h be the numerical solution corresponding to u . Assume that the finite element spaces Σ h and U h satisfy the following approximation properties:

(4.5) inf 𝝉 ¯ ∈ Σ h ⁡ ∥ 𝝌 ¯ - 𝝉 ¯ ∥ 𝐝𝐢𝐯 ≤ C ⁢ h s ⁢ ( ∥ 𝝌 ¯ ∥ s + ∥ 𝐝𝐢𝐯 𝝌 ¯ ∥ 1 + s ) ,

(4.6) inf 𝐯 ∈ U h ⁡ ∥ 𝐩 - 𝐯 ∥ 1 ≤ C ⁢ h s ⁢ ∥ 𝐩 ∥ 1 + s .

Then we have

∥ 𝐮 - 𝐮 h ∥ 0 ≤ Ch s ⁢ ∥ 𝐟 ∥ 0 .

Proof.

Consider the dual problem (4.3) and recall the stability bound (4.4). By choosing 𝐠 = 𝐮 - 𝐮 h and taking the test functions in (4.3) to be 𝝃 ¯ = 𝝈 ¯ - 𝝈 ¯ h , 𝐯 = 𝐮 - 𝐮 h , we have

∥ 𝐮 - 𝐮 h ∥ 0 2 = ( 𝒜 ⁢ 𝝌 ¯ , 𝒜 ⁢ ( 𝝈 ¯ - 𝝈 ¯ h ) ) + ( 𝐝𝐢𝐯 𝝌 ¯ , 𝐝𝐢𝐯 ( 𝝈 ¯ - 𝝈 ¯ h ) )

- ( 𝒜 ⁢ ( 𝝈 ¯ - 𝝈 ¯ h ) , ϵ ¯ ⁢ ( 𝐩 ) ) - ( 𝒜 ⁢ 𝝌 ¯ , ϵ ¯ ⁢ ( 𝐮 - 𝐮 h ) ) + ( ϵ ¯ ⁢ ( 𝐩 ) , ϵ ¯ ⁢ ( 𝐮 - 𝐮 h ) ) .

Now using the error equations of (2.7), for all 𝝉 ¯ h ∈ Σ h and 𝐯 h ∈ U h we arrive at

∥ 𝐮 - 𝐮 h ∥ 0 2 = ( 𝒜 ⁢ ( 𝝌 ¯ - 𝝉 ¯ h ) , 𝒜 ⁢ ( 𝝈 ¯ - 𝝈 ¯ h ) ) + ( 𝐝𝐢𝐯 ( 𝝌 ¯ - 𝝉 ¯ h ) , 𝐝𝐢𝐯 ( 𝝈 ¯ - 𝝈 ¯ h ) )

- ( 𝒜 ⁢ ( 𝝈 ¯ - 𝝈 ¯ h ) , ϵ ¯ ⁢ ( 𝐩 - 𝐯 h ) ) - ( 𝒜 ⁢ ( 𝝌 ¯ - 𝝉 ¯ h ) , ϵ ¯ ⁢ ( 𝐮 - 𝐮 h ) ) + ( ϵ ¯ ⁢ ( 𝐩 - 𝐯 h ) , ϵ ¯ ⁢ ( 𝐮 - 𝐮 h ) ) .

After using the Cauchy–Schwarz inequality, the definitions of ∥ ⋅ ∥ 𝐝𝐢𝐯 , ∥ ⋅ ∥ 1 , and collecting common terms, we arrive at

∥ 𝐮 - 𝐮 h ∥ 0 2 ≤ C ⁢ ( ∥ 𝝌 ¯ - 𝝉 ¯ h ∥ 𝐝𝐢𝐯 + ∥ 𝐩 - 𝐯 h ∥ 1 ) ⁢ ( ∥ 𝝈 ¯ - 𝝈 ¯ h ∥ 𝐝𝐢𝐯 + ∥ 𝐮 - 𝐮 h ∥ 1 )

≤ C ⁢ ( ∥ 𝝌 ¯ - 𝝉 ¯ h ∥ 𝐝𝐢𝐯 + ∥ 𝐩 - 𝐯 h ∥ 1 ) ⁢ ∥ 𝐟 ∥ 0 .

Applying the approximation properties of Σ h and U h to the above and using the regularity of the dual problem in (4.4), we get

∥ 𝐮 - 𝐮 h ∥ 0 2 ≤ C ⁢ h s ⁢ ∥ 𝐮 - 𝐮 h ∥ 0 ⁢ ∥ 𝐟 ∥ 0 . ∎

Remark 2.

The regularity assumed in (4.4) requires some explanation. Actually, it has been recently observed in [10] that a dual problem arising from least-squares formulations does not necessarily share the same regularity properties as the original problem. By inspecting the variational form (4.3), it can be observed that formally the solution of the dual problem satisfies the following strong formulation:

(4.7)

𝒜 ⁢ 𝝌 ¯ = ϵ ¯ ⁢ ( 𝐩 ) - 𝒞 ⁢ ∇ ⁡ 𝐰 ,

- 𝐝𝐢𝐯 𝝌 ¯ = 𝐰 ,

where 𝐰 satisfies

- 𝐝𝐢𝐯 ( 𝒞 ⁢ ϵ ¯ ⁢ ( 𝐰 ) ) = 𝐠 .

We can then deduce the existence of a suitable s so that (4.4) is satisfied. However, we cannot conclude that if 𝐮 is more regular, then the same regularity can be automatically transferred to the solution of the dual problem. We will go back to this discussion when estimating the rate of convergence of the approximation of the eigenvalues.

Remark 3.

The approximation properties (4.5) and (4.6) stated in Theorem 2 are satisfied by the most natural finite element spaces that can be used for the approximation of our problem. For instance Raviart–Thomas finite element enjoy (4.5) and standard Lagrangian finite element achieve (4.6). An interesting question is whether such estimates depend on the Lamé coefficients. Actually, the constants C are independent of the coefficients, while it can be expected that the norms of the dual solution might be affected by them.

The uniform convergence implies the convergence of the discrete eigenvalues which we recall in the next proposition (see [7] and [3]).

Proposition 3.

Assume that the operator T h converges in norm to T as h goes to zero, that is,

(4.8) ∥ T - T h ∥ ℒ ⁢ ( X ) → 0 .

Let κ ~ i = κ ~ i + 1 = … = κ ~ i + m - 1 be an eigenvalue of multiplicity m of the operator T. Then, for h small enough (so that the total number of discrete finite eigenvalues is greater than or equal to i + m - 1 ), the m discrete eigenvalues κ ~ j , h ( j = i , … , i + m - 1 ) of with T h converge to κ ~ j . Moreover, let E = ⊕ j = i i + m - 1 E j be the continuous eigenspace spanned by { u i , … , u i + m - 1 } , and let E h = ⊕ j = i i + m - 1 E j , h be the direct sum of the corresponding discrete eigenspaces spanned by { u i , h , … , u i + m - 1 , h } . Then the corresponding eigenfunctions converge, that is,

(4.9) δ ⁢ ( E , E h ) → 0 ,

where δ denotes the gap between Hilbert subspaces.

Moving to the rate of convergence, we start by stating the classical result for the convergence of the eigenfunctions (see, for instance, [3, Theorem 7.1]).

Theorem 4.

Let X be L 2 ⁢ ( Ω ) d or H 1 ⁢ ( Ω ) d . With the same notation and assumptions of Proposition 3 we have

δ ⁢ ( E , E h ) ≤ C ⁢ ∥ ( T - T h ) | E ∥ ℒ ⁢ ( X ) ,

where the gap δ is evaluated in the norm of X.

The estimation of the rate of convergence for the eigenvalues requires a more careful analysis since the discrete eigenvalue problem is not symmetric. To this aim, we embed our problem in the framework of general variationally posed eigenvalue problems as follows. A general variationally posed eigenvalue problem reads: given a Hilbert space 𝒳 and suitably defined bilinear forms 𝒜 ~ : 𝒳 × 𝒳 → ℝ and ℬ ~ : 𝒳 × 𝒳 → ℝ , find eigenvalues ξ ∈ ℂ and non-vanishing eigenfunctions 𝒰 ∈ 𝒳 such that

𝒜 ~ ⁢ ( 𝒰 , 𝒱 ) = ξ ⁢ ℬ ~ ⁢ ( 𝒰 , 𝒱 ) for all ⁢ 𝒱 ∈ 𝒳 .

We then introduce the solution operator 𝒯 : 𝒳 → 𝒳 associated with the general formulation so that if ℱ ∈ 𝒳 , then 𝒯 ⁢ ℱ ∈ 𝒳 is the solution of the problem

𝒜 ~ ⁢ ( 𝒯 ⁢ ℱ , 𝒱 ) = ℬ ~ ⁢ ( ℱ , 𝒱 ) for all ⁢ 𝒱 ∈ 𝒳 .

For the moment we proceed formally and assume that the operator 𝒯 is well defined. It is out of the scope of this paper to recall all the detail of the general theory which we are going to apply to our specific formulation. We identify 𝒳 with its dual so that the adjoint operator 𝒯 * : 𝒳 → 𝒳 is associated with the dual problem, that is, given 𝒢 ∈ 𝒳 we have that 𝒯 * ⁢ 𝒢 ∈ 𝒳 is equal to the solution of the problem

𝒜 ~ ⁢ ( 𝒱 , 𝒯 * ⁢ 𝒢 ) = ℬ ~ ⁢ ( 𝒱 , 𝒢 ) for all ⁢ 𝒱 ∈ 𝒳 .

Given a finite element space 𝒳 h ⊂ 𝒳 , we can consider the discrete eigenvalue problem: find eigenvalues ξ h ∈ ℂ and non-vanishing eigenfunctions 𝒰 h ∈ 𝒳 h such that

𝒜 ~ ⁢ ( 𝒰 h , 𝒱 ) = ξ h ⁢ ℬ ~ ⁢ ( 𝒰 h , 𝒱 ) for all ⁢ 𝒱 ∈ 𝒳 h .

Analogously, we can consider the discrete solution operator 𝒯 h : 𝒳 → 𝒳 that satisfies 𝒯 h ⁢ ℱ ∈ 𝒳 h and

𝒜 ~ ⁢ ( 𝒯 h ⁢ ℱ , 𝒱 ) = ℬ ~ ⁢ ( ℱ , 𝒱 ) for all ⁢ 𝒱 ∈ 𝒳 h ,

for ℱ ∈ 𝒳 and its adjoint 𝒯 h * : 𝒳 → 𝒳 which satisfies 𝒯 h * ⁢ 𝒢 ∈ 𝒳 h and

(4.10) 𝒜 ~ ⁢ ( 𝒱 , 𝒯 h * ⁢ 𝒢 ) = ℬ ~ ⁢ ( 𝒱 , 𝒢 ) for all ⁢ 𝒱 ∈ 𝒳 h

for 𝒢 ∈ 𝒳 .

We assume that the eigenmodes of 𝒯 h are approximating correctly the ones of 𝒯 and we discuss the rate of convergence. We consider now an eigenvalue ξ i = ξ i + 1 = … = ξ i + m - 1 of 𝒯 of multiplicity m and the corresponding discrete eigenvalues ξ j , h ( j = i , … , i + m - 1 ) of 𝒯 h approximating it. We denote by ℰ j the eigenspace of 𝒯 corresponding to ξ j and by ℰ j * the one of 𝒯 * . The following theorem summarizes the rate of convergence of the discrete eigenvalues (see, for instance, [3, Theorem 7.3]).

Theorem 5.

Let { U i , … , U i + m - 1 } form a basis for the eigenspace E = ⊕ j = i i + m - 1 E j and let { U i * , … , U i + m - 1 * } be the dual basis of the corresponding eigenspace E * = ⊕ j = i i + m - 1 E j * associated with T * . Then for j = i , … , i + m - 1 ,

(4.11) | ξ j - ξ j , h | ≤ C ⁢ { ∑ k , ℓ = 1 m | ( ( 𝒯 - 𝒯 h ) ⁢ 𝒰 k , 𝒰 ℓ * ) | + ∥ ( 𝒯 - 𝒯 h ) | ℰ ∥ ℒ ⁢ ( 𝒳 ) ⁢ ∥ ( 𝒯 * - 𝒯 h * ) | ℰ * ∥ ℒ ⁢ ( 𝒳 ) } .

Remark 4.

In the general case of non-symmetric problems, the previous theorem should take into account also the ascent multiplicities of the eigenvalues. In our case, however, the continuous problem is symmetric so that we know that the ascent multiplicities of our eigenvalues are always equal to one.

We now conclude this section by showing how the abstract results apply to our specific problem. As noted in Section 2, the eigenvalues ω i of the continuous problem (2.5) are real, positive and form an increasing sequence with each being repeated according to its multiplicity. One the contrary, the discrete eigenvalues ω i , h of (3.1) are complex and can be ordered according to their modules

0 ≤ | ω 1 , h | ≤ | ω 2 , h | ≤ | ω 3 , h | ≤ ⋯

and again repeated according to their multiplicities. The convergence of the eigenvalues and absence of spurious modes can, for instance, be stated as follows (see [7, Theorem 9.1]).

Theorem 6.

Let us assume that the convergence in norm (4.8) is satisfied. For all compact sets K in the complex plane that do not contain any eigenvalue of the continuous problem, there exists h 0 such that for all h < h 0 no eigenvalue of the discrete problem belongs to K.

As a consequence of the previous theorem, a more concrete representation of the eigenvalue convergence in the complex plane can be given as follows.

Property 1.

Let B be a ball with radius R > 0 and let S = { ω i } i = 1 n be the set of n (depending on R) eigenvalues, counted with their multiplicity, being inside B (i.e. | ω i | < R ). Then, for all R > 0 and all ϵ > 0 , there exists h 0 such that for h < h 0 exactly n discrete eigenvalues, counted with their multiplicity, are inside B, that is, | ω i , h | < R . Moreover, the n discrete eigenvalues can be sorted such that

(4.12) | ω i - ω i , h | < ϵ , i = 1 , 2 ⁢ … , n .

In [5] the following theorem was stated concerning the rate of convergence of the eigenvalues. We are now going to prove it rigorously by using the abstract setting that we have introduced.

Theorem 7.

Let ω i = ω i + 1 = … = ω i + m - 1 be an eigenvalue of (2.8) of multiplicity m and let E = ⊕ j = i i + m - 1 E j the corresponding eigenspace. Let E * be the corresponding eigenspace of the adjoint problem. If the discrete spaces Σ h and U h satisfy the approximation properties in Theorem 2, then for h small enough (so that dim ⁡ U h ≥ i + m - 1 ) the m discrete eigenvalues ω i , h , ω i + 1 , h , … , ω i + m - 1 , h of (3.1) converge to ω i . Denote by E h the direct sum of the discrete eigenspaces and consider the following quantities:

ρ ⁢ ( h ) = sup 𝐮 ∈ E ∥ 𝐮 ∥ = 1 ⁡ inf 𝝉 ¯ ∈ Σ h 𝐯 ∈ U h ⁡ ( ∥ 𝝈 ¯ - 𝝉 ¯ ∥ 𝐝𝐢𝐯 + ∥ 𝐮 - 𝐯 ∥ 1 ) ,

ρ * ⁢ ( h ) = sup 𝐩 ∈ E * ∥ 𝐩 ∥ = 1 ⁡ inf 𝝃 ¯ ∈ Σ h 𝐪 ∈ U h ⁡ ( ∥ 𝝌 ¯ - 𝝃 ¯ ∥ 𝐝𝐢𝐯 + ∥ 𝐩 - 𝐪 ∥ 1 ) ,

where 𝛔 ¯ is the other component of the solution of (2.8) corresponding to u and analogously 𝛏 ¯ is the other component corresponding to p for the dual problem. Then the following error estimates hold true:

δ ⁢ ( E , E h ) ≤ C ⁢ ρ ⁢ ( h ) , | ω j - ω j , h | ≤ C ⁢ ρ ⁢ ( h ) ⁢ ρ * ⁢ ( h ) j = i , … , i + m - 1 ,

where the gap δ is evaluated in the H 1 ⁢ ( Ω ) norm.

Proof.

The estimate for the eigenspaces is a direct consequence of Theorem 4 and of the a priori estimates for the approximation of the source problem (2.7). Indeed, since T is the solution operator associated with 𝐮 (with T ⁢ 𝐟 = κ ~ ⁢ 𝐟 = 𝐮 ), by the definition of the operator norm and the standard energy estimates, and by considering the other component 𝝈 ¯ of the solution of (2.7), we have

δ ⁢ ( E , E h ) ≤ C ⁢ ∥ ( T - T h ) | E ∥ ℒ ⁢ ( H 1 )

= C ⁢ sup 𝐟 ∈ E ∥ 𝐟 ∥ 1 = 1 ⁡ ∥ ( T - T h ) ⁢ 𝐟 ∥ 1

≤ C ⁢ sup 𝐮 ∈ E ∥ 𝐮 ∥ 1 = 1 ⁡ ∥ 𝐮 - 𝐮 h ∥ 1

≤ C ⁢ sup 𝐮 ∈ E ∥ 𝐮 ∥ 1 = 1 ⁡ ( ∥ 𝐮 - 𝐮 h ∥ 1 + ∥ 𝝈 ¯ - 𝝈 ¯ h ∥ 𝐝𝐢𝐯 )

≤ C ⁢ sup 𝐮 ∈ E ∥ 𝐮 ∥ 1 = 1 ⁡ inf 𝝉 ¯ ∈ Σ h 𝐯 ∈ U h ⁡ ( ∥ 𝐮 - 𝐯 ∥ 1 + ∥ 𝝈 ¯ - 𝝉 ¯ ∥ 𝐝𝐢𝐯 ) = C ⁢ ρ ⁢ ( h ) .

In order to show the estimate for the eigenvalues we are going to use the results presented in Theorem 5. The correspondence between the abstract setting and our problem is given by the following identifications:

𝒳 = 𝐗 ¯ N × H 0 , D 1 ⁢ ( Ω ) d , 𝒰 = ( 𝝈 ¯ , 𝐮 ) , 𝒱 = ( 𝝉 ¯ , 𝐯 ) ,

𝒜 ~ ⁢ ( 𝒰 , 𝒱 ) = ( 𝒜 ⁢ 𝝈 ¯ , 𝒜 ⁢ 𝝉 ¯ ) + ( 𝐝𝐢𝐯 𝝈 ¯ , 𝐝𝐢𝐯 𝝉 ¯ ) - ( 𝒜 ⁢ 𝝉 ¯ , ϵ ¯ ⁢ ( 𝐮 ) ) - ( 𝒜 ⁢ 𝝈 ¯ , ϵ ¯ ⁢ ( 𝐯 ) ) + ( ϵ ¯ ⁢ ( 𝐮 ) , ϵ ¯ ⁢ ( 𝐯 ) ) ,

ℬ ~ ⁢ ( 𝒰 , 𝒱 ) = - ( 𝐮 , 𝐝𝐢𝐯 𝝉 ¯ ) .

With these identifications we are now in a position to give a representation of the solution operator 𝒯 and of its adjoint 𝒯 * . It is interesting to remark that the bilinear form 𝒜 ~ is symmetric, so that the non-symmetry of the problem comes from the fact that the bilinear form on the right hand side ℬ ~ is not symmetric. For this reason, according to (4.10), the adjoint operator will be constructed by keeping the same left hand side of our equation and by considering the transpose of the right hand side.

More precisely, given ℱ = 𝒢 = ( 𝐆 ¯ , 𝐟 ) ∈ 𝒳 , 𝒯 ⁢ ℱ will be given by ( 𝝈 ¯ , 𝐮 ) ∈ 𝒳 solution of (2.7), while 𝒯 * ⁢ 𝒢 is given by the solution ( 𝝌 ¯ , 𝐩 ) ∈ 𝒳 of the following dual problem:

(4.13) { ( 𝒜 ⁢ 𝝌 ¯ , 𝒜 ⁢ 𝝃 ¯ ) + ( 𝐝𝐢𝐯 𝝌 ¯ , 𝐝𝐢𝐯 𝝃 ¯ ) - ( 𝒜 ⁢ 𝝃 ¯ , ϵ ¯ ⁢ ( 𝐩 ) ) = 𝟎 for all ⁢ 𝝃 ¯ ∈ 𝑿 ¯ N , - ( 𝒜 ⁢ 𝝌 ¯ , ϵ ¯ ⁢ ( 𝐯 ) ) + ( ϵ ¯ ⁢ ( 𝐩 ) , ϵ ¯ ⁢ ( 𝐯 ) ) = - ( 𝐝𝐢𝐯 𝐆 ¯ , 𝐯 ) for all ⁢ 𝐯 ∈ H 0 , D 1 ⁢ ( Ω ) d .

The definition of the discrete operators 𝒯 h and 𝒯 h * is completely analogous and makes use of the discrete space 𝒳 h = ( Σ h , U h ) .

Theorem 5 bounds the error in the eigenvalues with the sum of two terms. The bound for the first term uses the same arguments as in the proof of Theorem 2, together with the definitions of ρ ⁢ ( h ) and ρ * ⁢ ( h ) . Let us focus on the second term.

We have to estimate the product of the following two quantities:

∥ ( 𝒯 - 𝒯 h ) | ℰ ∥ ℒ ⁢ ( 𝒳 ) and ∥ ( 𝒯 * - 𝒯 h * ) | ℰ * ∥ ℒ ⁢ ( 𝒳 )

Actually, the result will follow by observing that the first term is bounded by C ⁢ ρ ⁢ ( h ) and the second one by C ⁢ ρ * ⁢ ( h ) .

If ( 𝝈 ¯ , 𝐮 ) is an eigenfunction in ℰ , then the following a priori estimate is valid:

∥ 𝝈 ¯ - 𝝈 𝒉 ¯ ∥ 𝐝𝐢𝐯 + ∥ 𝐮 - 𝐮 h ∥ 1 ≤ C ⁢ ρ ⁢ ( h ) ⁢ ∥ 𝐮 ∥ 0 .

Since 𝐮 is an eigenfunction, its L 2 ⁢ ( Ω ) -norm is bounded by the H 1 ⁢ ( Ω ) -norm and we have that ∥ 𝐮 ∥ 1 ≤ ∥ 𝒰 ∥ 𝒳 , so that we can conclude the desired result

∥ ( 𝒯 - 𝒯 h ) | ℰ ∥ ℒ ⁢ ( 𝒳 ) ≤ C ⁢ ρ ⁢ ( h ) .

The estimate for the adjoint operator is performed analogously observing that, if ( 𝝌 ¯ , 𝐩 ) is an eigenfunction in ℰ * , then the following a priori bound holds true:

∥ 𝝌 ¯ - 𝝌 𝒉 ¯ ∥ 𝐝𝐢𝐯 + ∥ 𝐩 - 𝐩 h ∥ 1 ≤ C ρ * ( h ) ∥ 𝐝𝐢𝐯 𝝌 ¯ ∥ 0 .

By the standard relations between the eigenspaces ℰ and ℰ * and the fact that 𝝌 ¯ is an eigenfunction of the dual problem we can finally conclude that

∥ ( 𝒯 * - 𝒯 h * ) | ℰ * ∥ ℒ ⁢ ( 𝒳 ) ≤ C ⁢ ρ * ⁢ ( h ) . ∎

Remark 5.

A natural question is whether the results of the previous theorem guarantee the double order of convergence for the approximation of the eigenvalues. The double order of convergence would follow if we can show that ρ * ⁢ ( h ) is asymptotically equivalent to ρ ⁢ ( h ) . This is a tricky question because it involves the regularity of the dual problem. In general, as already mentioned in Remark 2, we cannot expect from the solution of the dual problem the same regularity as for the original one. On the other hand, we do not need the regularity of a generic solution, but only of the solution corresponding to an eigenfunction. At the moment, only partial results in this direction are available; in particular, we can prove that the dual solution of the FOSLS approximation of the Poisson equation has the same regularity as the original one when an eigenfunction is considered [4]. The same result can be proved for some DPG formulations [6]. On the other hand, the numerical results of the next section confirm that in our examples the double order of convergence is achieved.

5 Numerical Results

In this section we present several numerical results related to the formulation that we have described. Some preliminary results were already present in [5]. However in that case, for simplicity, a special constitutive law was considered so that the problem was equivalent to a Stokes system. Now we solve the true linear elasticity problem and we are interested in showing how the reported numerical results are robust with respect to the Lamé constants, in particular when tending to the incompressible limit. The behavior of the computed spectrum will give us useful information on some properties of the solution operator.

We are not focusing on the efficiency of the numerical solver, but we are mainly interested in representing the computed spectrum in the complex plane as a function of the mesh and of the Lamé constants.

To this end, we consider a 2 ⁢ D linear elasticity problem with homogeneous Dirichlet boundary conditions, where the compliance tensor is given by

(5.1) 𝒜 ⁢ 𝝉 ¯ = 1 2 ⁢ μ ⁢ ( 𝝉 ¯ - λ 2 ⁢ λ + 2 ⁢ μ ⁢ tr ⁡ ( 𝝉 ¯ ) ⁢ 𝐈 ¯ ) .

In all our tests the Lamé constants are μ = 1 and λ varying from 1, 100, 10 4 to 10 8 . Moreover, different mesh structures are studied and investigated.

We start by presenting the rate of convergence for the first and second eigenvalues on different meshes. Then the spread of eigenvalues in the complex plane is explored for each mesh as it is refined.

In order to avoid artifacts introduced by the fact that we are rewriting the elasticity equation as a first order system, we considered only genuine eigenvalues of (2.8) in the sense that they correspond to eigenfunctions for which the displacement is not zero. These eigenvalues are related to the Schur complement described in (3.5).

A first order scheme based on Raviart–Thomas elements ( R ⁢ T 0 ) for the approximated space Σ h is considered with continuous Galerkin of degree one for the space U h , with h being the mesh size or step size.

These finite element spaces satisfy the conditions of Theorem 2. The general matrices in system (3.3) were constructed with FEniCS [1] and then matrices corresponding to the Schur complement in (3.4) were extracted and used to solve the eigenvalue problem.

5.1 Rate of Convergence on the Square Domain

In order to confirm the convergence rate stated in Section 4, we first consider a square domain Ω = ] 0 , 1 [ 2 and three different mesh sequences.

Figure 1

Meshes on the unit square with N = 4 .

(a)

Right

(b)

Crossed

(c)

Nonuniform

Table 1

Estimates of the first and second eigenvalues on a square.

	λ = 1	λ = 100	λ = 10 4	λ = 10 8
ω 1	37.266072200953786	52.31315105053875	52.3443693	52.344691
ω 2	37.2660721997643	91.4778227239564	92.11827609964527	92.12439336305897

Figure 2

Rate of convergence for the first and second eigenvalues on squared domain.

(a)

Right mesh

(b)

Crossed mesh

(c)

Nonuniform mesh

(d)

Right mesh

(e)

Crossed mesh

(f)

Nonuniform mesh

We use the following meshes: a non-symmetric uniform mesh (Right), a symmetric uniform mesh (Crossed), and a non-symmetric and non-structured mesh (Nonuniform). The integer N refers to the number of subdivisions on each side of the square.An example of such meshes is given in Figure 1. We computed with FEniCS the reference solutions providing an estimate of the first and second eigenvalues for the different cases we want to approximate. We used a high order scheme and an adaptive algorithm with the classical SOLVE-ESTIMATE-MARK-REFINE strategy. The values of the reference eigenvalues are reported in Table 1. Figure 2 shows the rate of convergence for both the first and the second approximated eigenvalues. Clearly both eigenvalues match the theoretical result in Theorem 7 with the rate being of second order on all meshes provided.

5.2 Rate of Convergence on the L-Shaped Domain

The second domain that is considered is the L-shaped domain obtained by removing a quarter of the square ] 0 , 1 [ 2 . It is known that on such domain singularities due to the re-entrant corner may occur. We also consider three sequences of triangulations: (Left), (Crossed or Uniform), and (Nonuniform) meshes. These are illustrated in Figure 3 with N = 4 . In order to compare the approximated eigenvalues with a reference value, the same adaptive code, as in the previous case, was used to find the estimate for the first and second eigenvalues for different λ. These estimated values are reported in Table 2. It is well known that the rate of convergence is reduced when the eigenfunction is singular. This is expected in particular in the case of the first eigenvalue.

Figure 3

Meshes for L-shaped domain with N = 4 .

(a)

Left

(b)

Uniform

(c)

Nonuniform

Table 2

Estimates of the first and second eigenvalues on an L-shaped domain.

	λ = 1	λ = 100	λ = 10 4	λ = 10 8
ω 1	54.36578831544661	127.990463	128.52363885700083	128.5293816640767
ω 2	69.08352886532845	147.44431322194393	148.06735898422232	148.07344440728022

Figure 4

Rate of convergence for the first and second eigenvalues L-shaped domain.

(a)

Left mesh

(b)

Uniform mesh

(c)

Nonuniform mesh

(d)

Left mesh

(e)

Uniform mesh

(f)

Nonuniform mesh

When looking at the computed results, a particular phenomenon is present when approximating the first eigenvalue close to the incompressible limit. This is clearly seen in Figures 4 (a), (b), and (c) for all meshes. Actually, the expected rate of convergence is achieved but there is a strange pre-asymptotic behavior when the value of λ increases. By looking at the computed eigenfunctions it looks like this phenomenon is related to a sort of locking caused by the location of the degrees of freedom in the proximity to the re-entered corner. When the system approaches the incompressible limit, the displacement presents a recirculation close to the re-entrant corner. A coarse mesh is not capable to represent such vortex. For λ = 1 this situation does not occur. In order to better describe this phenomenon, we report in Figures 5 (a) and (c) the eigenfunction corresponding to the first eigenvalue on the mesh after three refinements for λ = 100 . We see that the correct vortex is not well approximated, possibly because of the mesh being coarse so that there are not enough degrees of freedom close to the re-entrant corner. As we refine, specifically starting from the fifth iteration and onwards, the expected circulation is noticed as Figures 5 (b) and (d) illustrate. This might be the reason why we need a fine mesh to observe the expected rate for the first eigenvalue. On the contrary, the re-entrant corner has no negative effects on the second eigenvalue even for coarse meshes. This can be seen, for example, in Figure 6 which illustrates the second eigenfunction on a Nonuniform mesh with N = 16 and λ = 100 .

Figure 5

First eigenfunction for λ = 100 on L-shaped domain.

(a)

Left mesh with N = 16

(b)

Left mesh with N = 64

(c)

Nonuniform mesh with N = 16

(d)

Nonuniform mesh with N = 64

$Figure 6 Second eigenfunction for λ = 100 {\lambda=100} on Nonuniform L-shaped domain.$

Figure 6

Second eigenfunction for λ = 100 on Nonuniform L-shaped domain.

5.3 The Distribution of the Eigenvalues in the Complex Plane

In this subsection we discuss the distribution of the spectrum of the discrete operator in the complex plane. This is the main motivation of our paper and provide us with some properties of the solution operator that are not visible when considering the source problem. Needless to say, the knowledge of the spectrum of the operator has important consequences for the design of numerical schemes, for instance when a transient problem is approximated.

We consider the distribution of the eigenvalues in the case of the two domains previously considered (square and L-shaped) and the same range of Lamé constants.

We note that we are taking specific outer zoom with the values on the axes being large in order to show the spread of a large portion of the spectrum. Moreover, different scales are considered depending on each case in order to better highlight the behavior of the specific approximation.

Figure 7

Eigenvalues of a square on a refined Right mesh.

$(a) λ = 1 {\lambda=1}$

(a)

λ = 1

$(b) λ = 100 {\lambda=100}$

(b)

λ = 100

$(c) λ = 10 4 {\lambda=10^{4}}$

(c)

λ = 10 4

$(d) λ = 10 8 {\lambda=10^{8}}$

(d)

λ = 10 8

Figure 8

Eigenvalues of a square on a refined Crossed mesh.

$(a) λ = 1 {\lambda=1}$

(a)

λ = 1

$(b) λ = 100 {\lambda=100}$

(b)

λ = 100

$(c) λ = 10 4 {\lambda=10^{4}}$

(c)

λ = 10 4

$(d) λ = 10 8 {\lambda=10^{8}}$

(d)

λ = 10 8

We recall that the exact eigenvalues are real positive, so that we would expect that the discrete eigenvalues, for h small enough, distribute along the positive real axis of the complex plane. This is actually a naive interpretation of the convergence results presented in the previous sections. Indeed, by inspecting carefully the result presented in Property 1, it can be seen that computed eigenvalues can also be spread in different regions of the complex plane, as far as they lie outside a ball centered at the origin. The larger the ball, the smaller in general should h be in order to detect this behavior.

Figure 7 shows the spectrum of our operator on the Right mesh for different values of λ. Each plot reports with different colors the results computed on four successfully refined meshes.

As shown in Figure 7 (a), we observe that for a small value of λ, the spectrum is well behaving close to the positive real axis with some non-real eigenvalues appearing as we refine. In this case, all values are located in the right half of the complex plane, that is, they have positive real part.

When λ increases, we can see from Figure 7 that eigenvalues with negative real part show up. Also in this case some non-real eigenvalues are present, even if their imaginary part is not as large as we will see in the next examples.

Figure 9

Eigenvalues of a square on a refined Nonuniform mesh.

$(a) λ = 1 {\lambda=1}$

(a)

λ = 1

$(b) λ = 100 {\lambda=100}$

(b)

λ = 100

$(c) λ = 10 4 {\lambda=10^{4}}$

(c)

λ = 10 4

$(d) λ = 10 8 {\lambda=10^{8}}$

(d)

λ = 10 8

The results for the Crossed mesh are reported in Figure 8. In this case, all eigenvalues have positive real part and are located to the right side of the complex plane. The more we refine, the more the scheme produces complex eigenvalues moving away from the region of interest and diverging. For small values of λ (Figure 8 (a)), most eigenvalues are positive and real with some complex being present for finer meshes. As λ becomes larger in Figure 8, the spectrum is spread more and more with complex eigenvalues having positive real part.

Finally, the spectrum computed with the Nonuniform mesh has a different structure. The eigenvalues appear to be more spread all over the complex plane as Figure 9 shows. As in the previous observations, the values of the eigenvalues, for λ = 1 , are concentrated to the right of the complex plane. However, in this case, the eigenvalues are scattered much more when comparing them to the Right and Crossed meshes for λ = 1 . The spectrum gets surprisingly more symmetric with respect to the origin when moving towards the incompressible limit.

Figure 10

Eigenvalues of a Left L-shaped domain.

$(a) λ = 1 {\lambda=1}$

(a)

λ = 1

$(b) λ = 100 {\lambda=100}$

(b)

λ = 100

$(c) λ = 10 8 {\lambda=10^{8}}$

(c)

λ = 10 8

In what follows we present some results for the L-shaped domain. We only present the cases where λ = 10 r for ( r = 0 , 2 , 8 ) since there are no significant differences between r = 4 and 8.

We start by exploring the Left mesh structure in Figure 10. Looking at the spectrum in this case, we see a similar behavior as for the case of the Right mesh on the square.

Figure 11

Eigenvalues of a Uniform L-shaped domain

$(a) λ = 1 {\lambda=1}$

(a)

λ = 1

$(b) λ = 100 {\lambda=100}$

(b)

λ = 100

$(c) λ = 10 8 {\lambda=10^{8}}$

(c)

λ = 10 8

Figure 12

Eigenvalues of Nonuniform L-shaped domain.

$(a) λ = 1 {\lambda=1}$

(a)

λ = 1

$(b) λ = 100 {\lambda=100}$

(b)

λ = 100

$(c) λ = 10 8 {\lambda=10^{8}}$

(c)

λ = 10 8

For Uniform meshes, as Figure 11 shows, occurrence of complex eigenvalues for λ = 1 are more present than in the previous case. Moreover, complex eigenvalues appear for higher values of λ and grow when refining. The significant eigenvalues are those close to the origin which are positive and real or close to being real in the limit.

The Nonuniform mesh shows again a similar picturesque behavior as in the case of the square domain. Figure 12 shows the distribution of eigenvalues for this case. Even for small λ there are eigenvalues with large imaginary part in this case, although with positive real part.

Figure 13

Eigenvalues within a circle of radius R for Nonuniform mesh on a square and L-shaped domains.

(a)

Square domain

(b)

L-shaped domain

An interesting observation, that matches the convergence of eigenvalues in Property 1, is clearly observed on Nonuniform meshes. To show this, we choose a Nonuniform mesh with λ = 10 8 for both domains. As seen in Figure 13, when refining the mesh, eigenvalues are distributed allover the complex plane, within a circle of radius R centered at the origin. The more we refine, the larger the circle becomes, indicating that the number of positive and finite eigenvalues of interest are present inside that circle. The circle can be taken larger and larger as the mesh size h approaches zero.

Figure 14

Eigenvalues for meshes zooming to region close to the origin.

(a)

Zoom out

(b)

Intermediate zoom

(c)

Zoom in

(d)

Zoom out

(e)

Intermediate zoom

(f)

Zoom in

We conclude our numerical experiments with more pictures for different values of the zoom parameters, so that the convergence behavior corresponding to Property 1 can be better appreciated. Figure 14 reports eigenvalues of both domains with λ = 10 8 . The outer zoom lens in Figure 14 (a) and (d) gives a clear picture on the spread of eigenvalues in the complex plane. Looking at an intermediate zoom lens as Figure 14 (b) and (e) show, the structure of eigenvalues for both domains is still not aligned with the region of interest as the spread is apparent. On the other hand, a zoomed in picture as Figure 14 (c) gives eigenvalues, on a square, which are real and positive with respect to the previously reported pictures. Similar level of zoom show the same nice behavior in the case of the L-shaped domain (see Figure 14 (f)). Thus, the scheme produces real and positive eigenvalues in the region of interest with no spurious modes as expected.

Funding source: Deutsche Forschungsgemeinschaft

Award Identifier / Grant number: BE 6511/1-1

Funding statement: Fleurianne Bertrand gratefully acknowledges support by the Deutsche Forschungsgemeinschaft in the Priority Program SPP 1748 Reliable simulation techniques in solid mechanics, Development of non-standard discretization methods, mechanical and mathematical analysis under the project number BE 6511/1-1. Daniele Boffi is member of INdAM Research group GNCS and his research is partially supported by PRIN/MIUR and by IMATI/CNR.

References

[1] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes and G. N. Wells, The FEniCS project version 1.5, Arch. Numer. Software 3 (2015), no. 100. Search in Google Scholar

[2] L. Alzaben, F. Bertrand and D. Boffi, Computation of eigenvalues in linear elasticity with least-squares finite elements: Dealing with the mixed system, 14th World Congress on Computational Mechanics (WCCM ECCOMAS Congress 2020), https://www.scipedia.com/public/Alzaben_et_al_2021a. 10.23967/wccm-eccomas.2020.095Search in Google Scholar

[3] I. Babuška and J. Osborn, Eigenvalue problems, Handbook of Numerical Analysis, Vol. II, Handb. Numer. Anal. II, North-Holland, Amsterdam (1991), 641–787. 10.1016/S1570-8659(05)80042-0Search in Google Scholar

[4] F. Bertrand and F. Boffi, First order least-squares formulations for eigenvalue problems, preprint (2020), https://arxiv.org/abs/2002.08145; IMA J. Appl. Math., to appear. 10.1093/imanum/drab005Search in Google Scholar

[5] F. Bertrand and D. Boffi, Least-squares formulations for eigenvalue problems associated with linear elasticity, Comput. Math. Appl. 95 (2021), 19–27. 10.1016/j.camwa.2020.12.013Search in Google Scholar

[6] F. Bertrand, F. Boffi and H. Schneider, DPG approximation of eigenvalue problems, preprint (2020), https://arxiv.org/abs/2012.06623. Search in Google Scholar

[7] D. Boffi, Finite element approximation of eigenvalue problems, Acta Numer. 19 (2010), 1–120. 10.1017/S0962492910000012Search in Google Scholar

[8] Z. Cai and G. Starke, Least-squares methods for linear elasticity, SIAM J. Numer. Anal. 42 (2004), no. 2, 826–842. 10.1137/S0036142902418357Search in Google Scholar

[9] C. Carstensen and J. Storn, Asymptotic exactness of the least-squares finite element residual, SIAM J. Numer. Anal. 56 (2018), no. 4, 2008–2028. 10.1137/17M1125972Search in Google Scholar

[10] B. Keith, A priori error analysis of high-order LL* (FOSLL*) finite element methods, Comput. Math. Appl. 103 (2021), 12–18. 10.1016/j.camwa.2021.10.015Search in Google Scholar

Received: 2022-02-17

Accepted: 2022-02-18

Published Online: 2022-03-26

Published in Print: 2022-07-01

This work is licensed under the Creative Commons Attribution 4.0 International License.

On the Spectrum of an Operator Associated with Least-Squares Finite Elements for Linear Elasticity

Abstract

1 Introduction

2 The Continuous Problem

2.1 The Variational Formulation

3 Galerkin Discretization

Proposition 1.

Remark 1.

4 Convergence Analysis

Theorem 2.

Proof.

Remark 2.

Remark 3.

Proposition 3.

Theorem 4.

Theorem 5.

Remark 4.

Theorem 6.

Property 1.

Theorem 7.

Proof.

Remark 5.

5 Numerical Results

5.1 Rate of Convergence on the Square Domain

5.2 Rate of Convergence on the L-Shaped Domain

5.3 The Distribution of the Eigenvalues in the Complex Plane

References

Journal and Issue

Articles in the same Issue