Journal of Computational and Applied Mathematics

We construct finite element spaces of symmetric stress tensors that are exactly divergence-free. Moreover, their basis functions can be chosen so that they have small supports. These properties are highly desired in a number of important applications. Approximation properties of finite element spaces of divergence-free tensor functions are derived from properties of C 1 finite elements. © 2023TheAuthor(s).PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCC BYlicense(http://creativecommons.org/licenses/by/4.0/).


Introduction
In many problems of mathematical physics we have to deal with divergence-free functions, e.g., when solving Maxwell equations, Navier-Stokes equations, Einstein's equations, bending moments of a clamped plate, heat fluxes by a dual variational formulation of a steady-state heat conduction problem, and various continuity equations.In [1][2][3][4][5], we showed how to construct finite element spaces whose vector functions are exactly divergence-free in 2D and 3D.In this paper, we shall construct finite element spaces whose tensor functions are exactly divergence-free in two-dimensional domains.Such spaces are needed to calculate e.g. the stress tensor of a linear elasticity problem by means of a dual variational formulation.
Let Ω ⊂ R 2 be a bounded domain with Lipschitz boundary ∂Ω.We shall use the standard Sobolev space notation.
The symbol H k (Ω) stands for the Sobolev spaces of functions whose generalized derivatives up to the order k are square integrable over Ω.The norm in (H k (Ω)) p (p ≥ 1 is integer) is denoted by ∥ • ∥ k and the scalar product in the Lebesgue space (L 2 (Ω)) p is denoted by (•, •) 0 .The space of symmetric 2 × 2 tensors (L 2 (Ω)) 2×2 sym = {τ ∈ (L 2 (Ω)) 2×2  | τ = τ ⊤ } will be equipped with the scalar product Let us introduce the operator ε : ) , Denote by C ∞ 0 (Ω) the space of infinitely differentiable functions with compact support in Ω. 2 .Then we say that the divergence of the tensor function τ exists in the sense of distributions in Ω and define Div τ = d.
Obviously, for any smooth tensor τ we find that i.e., each component of Div τ is, in fact, the usual divergence of the corresponding row of τ .Due to the symmetry of τ we can also define the divergence of tensors by mean of columns.

Dual formulation of the linear elasticity problem
When solving real-life technical problems, the knowledge of the stress tensor is more important than the knowledge of displacements.We can, of course, first calculate the strain tensor by differentiating the computed displacements, and then construct the stress tensor by Hooke's law.This is the so-called primal approach.Alternatively, we can apply a dual formulation that allows us to calculate the stress tensor directly.
Let Γ D and Γ N be disjoint and relatively open sets of ∂Ω such that 2 be given body forces and g = (g 1 , g 2 ) ⊤ ∈ (L 2 (Γ N )) 2 given surface forces.In addition, if Γ D = ∅, we assume that the following equilibrium condition for forces f and g and their moments is satisfied: where P = {v ∈ (H 1 (Ω)) 2  | ε(v) = 0} is a three-dimensional space with basis (1, 0) ⊤ , (0, 1) ⊤ , and (x 2 , −x 1 ) ⊤ , see [6, p. 95].Define the set of statically admissible stresses as where and From (1) and (2) we find that Div τ + f = 0 in Ω and τ n = g on Γ N for a sufficiently smooth τ ∈ T (f , g).Now let us recall very briefly a primal formulation of the classical linear elasticity problem, see [6, p. 64].The strain tensor ε(v) is coupled with the stress tensor τ via the generalized Hooke's law where the linear elasticity coefficients and that there exists a constant c > 0 such that where ∥ • ∥ is the standard spectral norm.Further, assume that u ∈ (H 1 (Ω)) 2 is a suitable displacement function which satisfies the Dirichlet boundary conditions on Γ D , namely u = u on Γ D in the sense of traces.Then the primal problem of linear elasticity consists of finding u ∈ (H 1 (Ω)) 2 which minimizes the functional of potential energy I : (H 1 (Ω)) 2  → R 1 defined by In what follows, we concentrate on a dual variational formulation of the linear elasticity problem [6, p. 106].According to [6, p. 65], the generalized Hooke's law can be inverted (see also [7]).It gives a relation between strain and stress tensor for a nonhomogeneous and anisotropic material of the elastic body: A ijkl τ kl ) , and that there exists a constant c > 0 such that Definition.The dual problem of linear elasticity consists of finding a 2 × 2 symmetric tensor σ which minimizes the functional (of the complementary energy) J : Obviously, the bilinear form a(•, •) is symmetric and uniformly elliptic.Further, we introduce an equivalent formulation of the dual problem.Let some particular solution τ ∈ T (f , g) be given, i.e.Div τ + f = 0 in Ω and τ n = g on Γ N in the sense of distributions, where n is the outward unit normal to ∂Ω.Furthermore, we shall assume that all functions appearing below are sufficiently smooth so that the corresponding operations are correctly defined.Using the substitution τ = τ 0 + τ , we can reformulate the dual problem as follows: Find a 2 × 2 symmetric tensor σ 0 which minimizes the functional J 0 : (L 2 (Ω)) 2×2 sym → R 1 defined by over the space of divergence-free tensor functions If σ 0 ∈ T is sufficiently smooth, then Div σ 0 = 0 in Ω and σ 0 n = g on Γ N .Since τ is a given function, J(τ ) is a fixed constant.Then one immediately sees that and therefore, Remark 1.There are many possibilities how to construct a particular solution τ ∈ T (f , g) appearing in ( 4) and (6).For instance, when Γ N = ∅ we define and set τ 12 = τ 21 = 0, where ( Then from ( 1) and ( 2) we find that Div τ + f = 0 in Ω.For Γ N ̸ = ∅ see Appendix and [8, p. 447].

Construction of divergence-free tensor functions
First we introduce a method concerning how divergence-free (i.e.solenoidal) vector functions from the space can be described.The definition of Q is based on the following Green's formula: for a sufficiently smooth vector function q.In this case we have div q = 0 in Ω and q ⊤ n = 0 on Γ N if q ∈ Q due to density theorems for the Sobolev spaces (3), see [9,10].Note that the above Green's formula ( 8) is valid also for functions q whose divergence exists in the sense of distributions, see [11].Further, we introduce the operator where Theorem 1.Let Γ D and Γ N be connected.Then Q = curl W .
Proof.Let q ∈ Q .Then by Green's formula (8) we find that div q = 0 in Ω in the sense of distributions if all test functions v are equal to zero on whole boundary ∂Ω.Taking v ≡ 1 in ( 8) we find that ∫ ∂Ω q ⊤ n ds = 0.
According to [11, p. 22], there exists a stream function w ∈ H 1 (Ω) unique apart from an additive constant (which will be chosen later) such that q = curl w. ( Let Γ N ̸ = ∅ and let s = (n 2 , −n 1 ) ⊤ , i.e., s is the unit tangential vector to ∂Ω.Then by ( 7) and Green's formula (8) we for all ϕ ∈ C ∞ (Ω) such that ϕ = 0 on Γ D .This implies that the tangential derivative of w on Γ N is zero in the sense of distributions.Since Γ N is connected, we can choose the stream function w in ( 9) to be zero on Γ N .If Γ N = ∅, then the additive constant can be chosen arbitrarily.Hence, q ∈ curl W . Conversely, let w ∈ W and let ϕ ∈ C ∞ (Ω) with ϕ = 0 on Γ D .Similarly to (10) we get However, the last integral vanishes, since w = 0 on Γ N .Hence, curl w ∈ Q by (7).

□
Note that there are analogous results based on a theorem of De Rahm, see e.g.[12].Now recall the definition of the Hessian operator hes: ) .
Further, we introduce the operator inv : ) and let where ∂ n stands for the normal derivative.Notice that invz is the inverse tensor to hes z up to det(hes z).
Theorem 2. Let Γ D and Γ N be connected.Then where T is from (5).
Proof.Let τ = (τ ij ) ∈ T be arbitrary.Denote its columns by q 1 and q 2 .Clearly, q i ∈ Q for j = 1, 2. Since Γ D and Γ N are connected, we obtain by Theorem 1 that there exist two stream functions w 1 , w 2 ∈ W such that Due to the symmetry τ 12 = τ 21 we obtain −∂ 1 w 1 = ∂ 2 w 2 .Setting w = (w 1 , w 2 ), we see that w is also divergence-free.Hence, w ∈ Q , since w ⊤ n = 0 on Γ N .Using Theorem 1 once again, we find that there exists the third stream function z ∈ W such that w = curl z.
Since the derivatives ∂ 1 z and ∂ 2 z also belong to W , we get that z ∈ Z .Moreover, we see that τ = inv z, for instance,

Finite element approximation of the dual problem
Using Theorem 1, we can define finite element spaces of divergence-free vector functions by Q h = curl W h , where W h is an arbitrary finite element space of W .Similarly, Theorem 2 gives us a possibility to construct finite element spaces of divergence-free tensor functions as follows.Let Z h ⊂ Z be a finite element space generated by C 1 -elements over some triangulation [13,14].Then we set Remark 3.For construction of finite element spaces Z h we refer to [14].If basis functions {z i } in Z h , i = 1, . . ., dim Z h , have small supports, then {inv z i } ⊂ T h have also small supports.We can take for instance piecewise quadratic composite C 1 -triangular elements invented by G. Heindl [15] (see also [4,16]).Applying the operator inv to any function z h ∈ Z h , we get a piecewise constant tensor inv(z h ) ∈ T h which is exactly divergence-free in the sense of distributions.The basis functions inv z i can be calculated analytically by double differentiation of C 1 piecewise quadratic basis functions z i .Then we can easily handle, e.g., possible inequality constraints that occur in plasticity and related limit or shakedown analysis, in problems of Signorini's type or in stress hardening [6].We could also apply the composite piecewise cubic Hsieh-Clough-Tocher C 1 triangular elements (see [17, p. 341]).Then the corresponding divergence-free basis functions inv z i are piecewise linear and discontinuous, in general.
A conforming finite element approximation of the dual linear elasticity problem consists of finding σ 0 h ∈ T h which minimizes the same functional (4) over T h and the function is called then an approximate solution.By ( 6), (10), and the well-known Cea's lemma there exists a constant where σ 0 = inv z 0 , z 0 ∈ Z , and τ h = inv z h due to Theorem 2. Now we can apply standard interpolation results for C 1 -elements to prove the convergence (or some rate of convergence) of σ h to σ as h → 0, see e.g.[14,17].

Conclusions
Divergence-free vector or tensor functions have many important application.For instance, a simultaneous use of the primal and dual formulations enables us to apply the hypercircle method which produces exact evaluation of the discretization error when the primal and dual finite element method are employed, see [6, p. 259].We can also obtain two-sided bounds of energy and a posteriori error estimates, see [18, p. 64].Some other applications are mentioned in Introduction.
In this paper, we show how to construct finite element stress tensor basis functions that are exactly divergence-free and have small supports which is necessary to get sparse resulting systems of linear algebraic equations.Our construction is based on a special second order differential operator inv which is applied to C 1 finite elements.The convergence and the rate of convergence of the corresponding finite element approximations follow directly from approximation properties of C 1 elements used.

□Remark 2 .
and thus, inv z ∈ T .For a given divergence-free tensor τ ∈ T the corresponding potential function z such that τ = inv z is called the Airy function, see[6, p. 161].It is unique apart from a linear function.To see this we take z 1 , z 2 ∈ Z such that inv z 1 = inv z 2 = τ .Then inv(z 1 − z 2 ) = 0 implying that z 1 − z 2 is a linear polynomial.