A posteriori estimates for conforming Kirchhoff plate elements
Tom Gustafsson Rolf Stenberg Juha Videman
math.NA 65N30


Gaussian Processes for Individualized Continuous Treatment Rule Estimation
    Pavel Shvechikov 

 and Evgeniy Riabenko The reported study was funded by RFBR, according to the research project No. 16-07-01192 А.

 Department of Computer Science, Higher School of Economics
    December 30, 2023
=====================================================================================================================================================================================================




 We derive a residual a posteriori estimator for the Kirchhoff plate bending problem. We consider the problem with a combination of clamped, simply supported and free boundary conditions subject to both distributed and concentrated (point and line) loads. Extensive numerical computations are presented to verify the functionality of the estimators.


 
    Kirchhoff plate, C^1 elements, a posteriori estimates


65N30




§ INTRODUCTION

The purpose of this paper is to perform an a posteriori error analysis of conforming finite element methods for the classical Kirchhoff plate bending model. So far this has not been done in full generality as it comes to the boundary conditions. Most papers deal only with clamped or simply supported boundaries, see <cit.> for conforming C^1 elements,  <cit.> for the mixed Ciarlet–Raviart method (<cit.>), and <cit.> for discontinuous Galerkin (dG) methods. The few papers that do address more general boundary conditions, in particular free, are <cit.> in which the nonconforming Morley element is analysed, <cit.> where a new mixed method is introduced and analysed, and <cit.> where a continuous/discontinuous Galerkin method is considered. One should also note that the Ciarlet–Raviart method cannot even be defined for general boundary conditions. Free boundary conditions could be treated using dG methods following an analysis similar to the one presented here. 

In this study, we will derive a posteriori estimates using conforming methods and allowing for a combination of clamped, simply supported and free boundaries. In addition, we will investigate the effect of concentrated point and line loads, which are not only admissible in our H^2-conforming  setting but of great engineering interest, on our a posteriori bounds in numerical experiments. We note that finite element approximation of elliptic problems with loads acting on lower-dimensional manifolds has been considered by optimal control theory, see <cit.> and all the references therein.

The outline of the paper is the following. In Section 2, we recall the Kirchhoff-Love plate model by presenting its variational formulation and the corresponding boundary value problem. We perform this in detail for the following reasons. First, as noted above, general boundary conditions are rarely considered in the numerical analysis literature. 
Second, the free boundary conditions consist of a vanishing normal moment and a vanishing Kirchhoff shear force. 
These arise from the variational formulation via successive integrations by parts. It turns out that the same 
 steps are needed in the a posteriori analysis in order to obtain a sharp estimate, i.e. both reliable and efficient. 
In the following two sections, we present the classical conforming finite element methods and derive new a posteriori error estimates. In the last section, we present the results of our numerical experiments computed with the triangular Argyris element. We consider the point, line and square load cases with simply supported boundary conditions in a square domain as well as solve the problem in an L-shaped domain with uniform loading  using different combinations of boundary conditions.



§ THE KIRCHHOFF PLATE MODEL


The dual kinematic and force variables in the model are the curvature and the
moment tensors. Given the deflection u of the midsurface of the plate, the
curvature is defined through 

    Κ(u) = -(∇ u),

with the infinitesimal strain operator defined by

    (v)
        = 1/2*∇v+∇v^T,

where (∇v)_ij = ∂ v_i/∂ x_j.
The dual force variable, the moment tensor ,  is related to
Κ through the  constitutive relation

    (u) = d^3/12Κ(u),

where d denotes the plate thickness and where we have assumed an isotropic linearly elastic material, i.e.

    A
        = E/1 + ν*A
                   + ν/1 - ν(tr A) I,
            ∀A∈ℝ^2 × 2.

Here E and ν are the Young's modulus and Poisson ratio, respectively. 
The shear force is denoted by
=(u). The moment equilibrium equation reads as 

    (u) = (u),

where  is the vector-valued divergence operator applied to tensors.
The transverse shear equilibrium equation is

    -(u) = l,

with  l denoting the transverse loading. Using the constitutive
relationship (<ref>), a straightforward elimination yields the well-known Kirchhoff–Love plate equation:

    (u) := D Δ^2 u = l,

where the so-called bending stiffness D is defined as

    D = E d^3/12(1-ν^2).




Let Ω⊂ℝ^2 be a  polygonal domain that describes
the midsurface of the plate.  The plate is considered to be clamped on ⊂∂Ω, simply supported on ⊂∂Ω and
free on ⊂∂Ω as depicted in Fig. <ref>. The
loading is assumed to consist of a distributed load f∈ L^2(Ω), a load
g∈ L^2(S) along the line S⊂Ω, and a point load F at an interior 
point x_0 ∈Ω.

Next, we will  turn to the boundary conditions, which are best understood from
the variational formulation. (Historically, this was also how they were first
discovered by Kirchhoff, cf. <cit.>.) The elastic energy of the plate as a function of the
deflection v is 1/2 a(v,v), with the bilinear form a defined by

    a(w,v) = ∫_Ω(w) : Κ(v)  dx = ∫_Ωd^3/12(∇ w) : (∇ v) dx,

and the potential energy due to the loading is

    l(v)= ∫_Ω f v  dx + ∫_S g v  d s + F v(x_0).

Defining the space of kinematically admissible deflections

    V={  v ∈ H^2(Ω) :  v|_∪=0,  ∇ v ·n|_=0      },
 
minimization of the total energy 

    u=_v∈ V{1/2 a(v,v)-l(v)   }

leads to the following problem formulation.



[Variational formulation]
Find u ∈ V such that

    a(u,v) = l(v)   ∀ v ∈ V.





To derive the  corresponding boundary value problem, we recall the following
integration by parts formula, valid in any domain ⊂Ω   

    ∫_   (w) : Κ(v)  dx 
       = ∫_(w) ·∇ v  dx - ∫_∂(w)n·∇ v ds 
       = ∫_(w) v  dx + ∫_∂(w) ·n v ds - ∫_∂(w)n·∇ v ds ,

At the boundary ∂ R, the correct physical quantities are the components
in the normal n and tangential s directions. 
Therefore, we write 

    ∇ v=∂ v/∂ nn + ∂ v/∂ ss

and  define the normal shear force and the normal and twisting moments as

    (w)=(w)·n,   (w) =  n·(w) n, 
       (w) =(w)=
     s·(w) n .

 With this notation, we can write
 
    ∫_∂(w) ·n v ds - ∫_∂(w)n·∇ v ds
          
       =  ∫_∂(w)v  ds - ∫_∂ ((w)∂ v/∂ n+    (w)∂ v/∂ s)ds,

and thus rewrite the integration by parts formula (<ref>) as

    ∫_   (w) : Κ(v)  dx 
       = ∫_(w) v  dx + 
               ∫_∂(w)v  ds
               
           - ∫_∂ ((w)∂ v/∂ n+    (w)∂ v/∂ s)ds.


The key observation for deriving the correct boundary conditions is that, at any boundary point, a
value of v specifies also ∂ v /∂ s.
Defining the  Kirchhoff shear force (cf. <cit.>)

    (w)= (w) +∂(w)/∂ s

an integration by parts on a smooth part  of ∂ R yields
 
    ∫_(w)v  ds
     - ∫_ (w)∂ v/∂ sds
     = ∫_(w)   ds -|_a^b(w) v,

where a and b are the endpoints of . 



We are now in the position to state the boundary value problem for the Kirchhoff plate model.
Assuming   a smooth solution u in (<ref>), we have
  
    a(u,v)=   ∫ _Ω𝒜(u)  v  dx+  ∫_∂Ω(u)v  ds
     
        -∫_∂Ω ((u)∂ v/∂ n+    (u)∂ v/∂ s)ds.

With the combination of clamped, simply supported and free boundary  conditions at
∂Ω = ∪∪,  we   have for any v∈ V,

    ∫_∂Ω(u)v  ds - ∫_∂Ω ((u)∂ v/∂ n+    (u)∂ v/∂ s)ds
    
       = ∫_(u)v  ds
     - ∫_ (u)∂ v/∂ sds
       - ∫_∪ (u)∂ v/∂ nds.

In the final step, we integrate by parts at the free part of the boundary. To this end, let =∪_i=1^m+1^i, with ^i smooth.
Integrating by parts over ^i yields

    ∫_^i (u)v  ds
     - ∫_^i (u)∂ v/∂ sds
     = ∫_^i (u) v  ds -|_c_i-1^ c_i(u) v

where c_0 and c_m+1 are the end points of  and c_i, i=1,…, m, its successive interior corners.  Combining equations (<ref>)–(<ref>), and noting that v(c_0)=v(c_m+1)=0, gives finally 

    a(u,v)=∫ _Ω   𝒜(u)  v  dx   
     - ∫_∪ (u)∂ v/∂ nds
    
       
     + ∑_i=1^m+1∫_^i (u)   v  ds - ∑_i=1^m{( (u)|_c_i+ -(u)|_c_i-} v(c_i),

 where M_ns(u)|_c_i ± = lim_ϵ→ 0 +  M_ns(u)|_c_i + ϵ(c_i± 1 - c_i)
 
 


Choosing  v∈ V in such a way that three of the four terms in (<ref>) vanish and the test function in the fourth  term remains arbitrary and repeating this for each term,  we arrive at the following boundary value problem:


  * In the domain we have the distributional differential equation

    𝒜(u) = l   Ω,

where l is the distribution defined by (<ref>).



  * On the clamped part we have the conditions

    u=0   ∂ u/∂ n=0   .



  * On the simply supported part it holds

    u=0  (u)=0    .
 


  * On the free part  it holds 

    (u)=0 (u)=0    ^i ,  i=1,…, m.


  * At the interior corners on the free part, we have the   matching condition on the twisting moments

    (u)|_c_i+ =(u)|_c_i-   c_i,   i=1,…,m.

 

 



 




§ THE FINITE ELEMENT METHOD AND THE A POSTERIORI ERROR ANALYSIS
 


The finite element method is defined on a mesh  consisting of shape regular triangles.
We assume that the point load is applied on a node of the mesh. Further, we
assume that the triangulation is such that the applied line load is on element
edges. We denote the edges in the mesh by   and divide them into the
following parts: the edges in the interior ^i, the edges on the curve of
the line load ^S ⊂^i, and the edges on the free and simply supported boundary,
^f and ^s, respectively. The conforming finite element space is
denoted by V_h. Different choices for V_h are presented in Section <ref>. Note that we often write a ≲ b (or a ≳ b ) when a ≤ Cb (or  a ≥ Cb) for some positive constant C independent of the finite element mesh.



[The finite element method]
    
    Find u_h ∈ V_h such that
    
    a(u_h,v_h) = l(v_h)   ∀ v_h ∈ V_h.



Let K and K' be two adjoining triangles with normals n and n', respectively, and  with the common edge  E=K∩ K'. On E we define the following jumps 

    (v)  |_E =   (v)  - M_n'n'(v)

and

    (v)  |_E =   (v)   +V_n'(v).







In the analysis, we will need the Girault–Scott <cit.> interpolation operator   Π_h : V →
V_h for which  the following estimate  holds
 
    ∑_K ∈ h_K^-4w-Π_h w_0,K^2 +∑_E ∈ h_E^-1∇(w-Π_h w)_0,E^2
         +∑_E ∈ h_E^-3w-Π_h w_0,E^2 ≲w_2^2     and    Π_h w_2 ≲w_2.

Note that the Girault–Scott  interpolant uses point values at the
vertices of the mesh.
We use this property in the proof of Theorem <ref>  to
derive a proper upper bound for the error in terms of the  edge  residuals.


Next, we formulate an a posteriori estimate for
Problem <ref>. 
The local error indicators are the following:


  * The residual on each element

    h_K^2 (u_h) -f _0,K ,    K∈.



  * The  residual of the normal moment jump along interior edges  

    h_E^1/2(u_h) _0,E ,    E∈^i .


  * The residual of the jump in the effective shear force along interior   edges  

    h_E^3/2(u_h)  - g _0,E,      E∈^S,
 

    h_E^3/2(u_h)  _0,E,     E∈^i∖^S.
 

  * The   normal moment  on edges at the free and simply supported boundaries

    h_E^1/2(u_h)  _0,E,    E∈^f∪^s.


  * The  effective shear force along edges at the free boundary

    h_E^3/2(u_h)   _0,E,     E∈^f.
 

The global error estimator is then defined through

    η^2=   ∑_K∈  h_K^4 (u_h) -f _0,K^2 + ∑_ E∈^S h_E^3 (u_h)   -g_0,E^2
        
          
       + ∑_E∈^i∖^S h_E^3 (u_h)  _0,E^2    
         +∑_ E∈^i h_E (u_h) _0,E^2
          
          
       + ∑_ E∈^f h_E^3 (u_h)   _0,E^2    
         +∑_ E∈^f∪^s h_E (u_h)  _0,E^2.
 

    The following estimate holds
    
    u-u_h_2 ≲η.

 


    Let w = u-u_h and w := Π_h w ∈ V_h be its interpolant.
    In view of the well-known coercivity of the bilinear form a and Galerkin orthogonality, we have
    
    u-u_h_2^2 
                           ≲ a(u-u_h,w) 
                        = a(u-u_h,w-w) 
                        
       = l(w-w)-a(u_h,w-w) .

   Since  x_0 is a mesh node and  the interpolant uses nodal values, we have
    
    F (w(x_0)-w(x_0))=0,
 
    and hence
    
    l(w-w)=(f, w-w)      +⟨ g, w-w⟩_S .

From integration by parts over the element edges, using the fact that the
interpolant uses values at the nodes, it then follows that
   
    u-u_h_2^2 
       ≲  (f, w-w)      +⟨ g, w-w⟩_S-a(u_h,w-w) 
                        
       = (f, w-w)      +⟨ g, w-w⟩_S 
                        
            
                        -  ∑_K ∈{((u_h),w-w)_K +  
                        ⟨(u_h) ,w-w⟩_∂ K
                -
                     ⟨(u_h),∂∂ s (w-w)⟩_∂ K 
                     -   ⟨(u_h),∂∂ n (w-w)⟩_∂ K}
       = (f, w-w)      +⟨ g, w-w⟩_S
                          
           
                          -  ∑_K ∈{((u_h),w-w)_K +  
                        ⟨(u_h) ,w-w⟩_∂ K
            - ⟨(u_h),∂∂ n (w-w)⟩_∂ K}
                        .

    Regrouping and recalling definitions (<ref>) and (<ref>),  yields
    
    u-u_h_2^2
                         
    ≲   ∑_K ∈(f-(u_h),w-w)_K  
                          
         
                       -∑_E∈^S⟨(u_h)  - g ,w-w⟩_E -
                       ∑_E∈^i ∖^S⟨(u_h)  ,w-w⟩_E
          - ∑_E∈^i⟨(u_h), ∂∂ n_E (w-w)⟩_E
          -    ∑_E∈^f⟨(u_h) ,w-w⟩_E
     - ∑_ E∈^f∪^s⟨(u_h), ∂∂ n_E (w-w)⟩_E 
                        .

    The asserted a posteriori estimate now follows by applying the Cauchy–Schwarz inequality and the interpolation estimate (<ref>).
  
 
   
Instead of the jump terms in the estimator  η, we could consider the normal and twisting moment jumps
   
    h_E^1/2(u_h)_0,E,    h_E^1/2(u_h) _0,E,

   and the normal shear force jumps
   
    h_E^3/2(u_h)  _0,E,     
     h_E^3/2(u_h)  - g _0,E.
 
In this case we cannot, however, prove the efficiency, i.e. the lower bounds.

Next, we will consider the question of efficiency.
Let f_h∈ V_h be the interpolant of f and define

    _K(f) = h_K^2‖ f-f_h‖_0,K.

Similarly, for   a polynomial approximation  g_h of g on E ⊂ S we define

 


    _E(g) = h_E^3/2‖ g-g_h‖_0,E.

 In the following theorem,  ω_E stands for the union of elements sharing an edge E.
In its proof, we will adopt some of the techniques used in <cit.>.

    For all v_h ∈ V_h it holds
    
    h_K^2 (v_h) - f_0,K   ≲u-v_h_2,K + _K(f),  K∈, 
     
           h_E^1/2(v_h) _0,E   ≲u-v_h_2,ω_E  
                                                                                                              + ∑_K ⊂ω_E_K(f),    E∈^i,
      
     
           h_E^3/2(v_h) _0,E   ≲u-v_h_2,ω_E 
     + ∑_K ⊂ω_E_K(f),    E∈^i∖^S,
     
    
     h_E^3/2(v_h)  -g_0,E   ≲u-v_h_2,ω_E 
     + ∑_K ⊂ω_E_K(f)+_E(g),  E∈^S,
           
    
      h_E^1/2(v_h)   _0,E   ≲u-v_h_2,ω_E  
                                                                                                              + ∑_K ⊂ω_E_K(f),    E∈^f∪^s,
            
    
            h_E^3/2(v_h)   _0,E   ≲u-v_h_2,ω_E 
        + ∑_K ⊂ω_E_K(f),    E∈^f.

    

    Denote by b_K ∈ P_6(K) the sixth order bubble that, together with its
    first-order derivatives, vanishes on ∂ K, i.e. let b_K =
    (λ_1,Kλ_2,Kλ_3,K)^2, where λ_j,K are
    the barycentric coordinates for K. Then we define
    
    γ_K = b_K h_K^4((v_h)-f_h)  in K  and  γ_K=0  in Ω∖ K,

    for  v_h ∈ V_h. The problem statement gives
    
    a_K(u,γ_K) = (f,γ_K)_K,

    where a_K(u,γ_K) = ∫_K (u) : Κ(γ_K)  dx. We have
    
    h_K^4 (v_h)-f_h_0,K^2
                   ≲ h_K^4 √(b_K)((v_h)-f_h) _0,K^2 
       = ((v_h)-f_h,γ_K)_K 
       = ((v_h),γ_K)_K - (f,γ_K)_K + (f-f_h,γ_K)_K 
    
                   = a_K(v_h-u,γ_K)+(f-f_h,γ_K)_K.

The local bound (<ref>) now follows from applying the continuity of a, the Cauchy–Schwarz inequality and inverse estimates.
    

    Next, consider inequality (<ref>). Suppose  E=K_1 ∩ K_2 for  the triangles K_1 and K_2; thus  ω_E=K_1∪ K_2. Let λ_E ∈ P_1(ω_E) be the linear
    polynomial satisfying
    
    λ_E|_E=0,  and  ∂λ_E/∂ n_E =1,

  and let p_1 be the polynomial that satisfies p_1 |_ E= 
  M_nn(v_h) |_E and ∂ p_1/∂ n_E|_E=
  0. Moreover, let p_2 ∈ P_8(ω_E) be the eight-order bubble that
  takes value one at the midpoint of the edge E and, together with its first-order
  derivatives, vanishes on ∂ω_E.
    Define  w = λ_E p_1 p_2. Since    
    
    ∂ w/∂ n_E|_E = ∂λ_E/∂ n_E M_nn(v_h)  p_2= M_nn(v_h)  p_2,

      scaling yields the equivalence
    
    M_nn(v_h) _0,E^2    ≈∂ w/∂ n_E_0,E^2
                ≈√( p_2)  M_nn(v_h) _0,E^2
               
       
       = ⟨ M_nn(v_h),∂ w∂ n_E⟩_E.
 
    Furthermore, since 
     
    ∂ w/∂ s|_E = 0,    w|_E ∪∂ω_E=0   and  ∇ w|_∂ω_E =0,

  the integration by parts formula  (<ref>) yields
     
    ⟨ M_nn(v_h),∂ w∂ n_E⟩_E 
            
                                                      =-∫_ω_E(v_h):Κ(w) dx + ((v_h),w)_ω_E.

        Extending w by zero to Ω∖ω_E,  we obtain from the problem statement (<ref>)  
        
    ∫_ω_E(u):Κ(w) dx -(f,w)_ω_E=0.

        Hence, using 
        
         the Cauchy–Schwarz inequality, we get from (<ref>)
        
    ⟨ M_nn   (v_h), ∂ w∂ n_E⟩_E 
            
       =  ∫_ω_E(u-v_h):Κ(w) dx+  ((v_h)-f,w)_ω_E
       ≲u-v_h_2,ω_E |w|_2,ω_E + (v_h)-f_0,ω_Ew_0,ω_E.

    By scaling, one easily shows that
            
    |w|_2,ω_E≲    h_E^-1/2∂ w/∂ n_E_0,E  and  w_0,ω_E≲  h_E^3/2∂ w/∂ n_E_0,E.

    The estimate (<ref>) 
     then follows from (<ref>),  (<ref>), (<ref>), (<ref>), and the already proved bound (<ref>).
       

    Since (<ref>) follows from (<ref>) with
    g=0, we prove the latter.  Due to the regularity condition imposed on the
    mesh there exists for each edge E a symmetric pair of smaller triangles
    (K^'_1,K^'_2) that satisfy ω^'_E=K_1^'∪
    K_2^'⊂ω_E, see Fig. <ref>.      Let
    w^' = p_2^' ((v_h) -g_h) where
    p_2^' is the eight-order bubble that takes value one at the midpoint
    of E and, together its first-order derivatives, vanishes on ∂ω_E^'.
    By the norm equivalence, we first have
   
    (v_h)   -g_h_0,E^2   ≈   w^'_0,E^2
            ≲⟨(v_h)  -g_h, w^'⟩_E .

    Next, we write 
       
    ⟨(v_h)  -g_h, w^'⟩_E
            =
             ⟨(v_h)  -g, w^'⟩_E +  ⟨ g -g_h, w^'⟩_E.

     Due to  symmetry, ∂ w^'/∂ n|_E = 0, and hence (<ref>) and (<ref>) give
    
    ⟨(v_h)  -g, w^'⟩_E = ⟨(v_h), w^'⟩_E-⟨   g, w^'⟩_E
                    
       =∫_ω_E^'(v_h):Κ(w^' ) dx-  ((v_h),w^' )_ω_E^' 
                    -⟨   g, w^'⟩_E.

               Extending w^' by zero to Ω∖ω^'_E, the variational form (<ref>) implies that
       
    ∫_ω_E^'(u):Κ(w^') dx -(f,w^')_ω_E^'-⟨ g,w^'⟩_E=0    .

        Hence, 
         
    ⟨(v_h)  -g, w^'⟩_E=  ∫_ω_E^'(v_h-u):Κ(w^') dx +(f-(v_h) ,w^')_ω_E^'

        and the Cauchy–Schwarz inequality, scaling estimates and (<ref>) give
        
    ⟨(v_h)  -g, w^'⟩_E    
       ≲ h_E^-3/2( u-v_h_2,ω_E^' +h_K^2 (v_h)-f_0,ω_E^')  w^'_0,E
       ≲ h_E^-3/2( u-v_h_2,ω_E + ∑_K ∈ω_E_K(f) )
              w^'_0,E  .

 The asserted estimate then follows from (<ref>), (<ref>) and (<ref>).   
     
     
     The estimates  (<ref>),
 (<ref>), are proved similarly to the bounds  (<ref>) and  (<ref>), respectively.

The above estimates provide  the following global bound.
  
It holds
    
    η≲u-u_h_2 + (f)+ (g),

where 

    (f) = √(∑_K ∈_K(f)^2 )    (g) = √(∑_E∈^S_E(g)^2).







§ THE CHOICE OF V_H



Let us briefly discuss some possible choices of conforming finite elements for
the plate bending problem. Each choice consists of a polynomial space
𝒫 and of a set of N degrees of freedom defined through a functional ℒ : C^∞→ℝ. We denote by x^k, k∈{1,2,3}, the vertices of the
triangle and by e^k, k∈{1,2,3}, the midpoints of the edges, i.e.

    e^1 = 1/2(x^1+x^2),   e^2 = 1/2(x^2+x^3),   e^3 =
        1/2(x^1+x^3).


The simplest H^2-conforming triangular finite element that is locally H^4(K)
in each K is the Bell triangle.


    
    𝒫   = { p ∈ P_5(K) : ∂ p∂ n∈ P_3(E) ∀ E ⊂ K }
    ℒ(w)    = 
            w(x^k),    for 1≤ k ≤ 3,
    ∂ w/∂ x_i(x^k),    for 1 ≤ k ≤ 3 and 1 ≤ i ≤ 2,
    ∂^2 w/∂ x_i ∂ x_j(x^k),    for 1 ≤ k ≤ 3 and 1 ≤ i, j ≤ 2.



Even though the polynomial space associated with the Bell triangle is not the whole P_5(K) it is still larger
than P_4(K). This can in some cases complicate the implementation.
Moreover, the asymptotic interpolation estimates for P_5(K) are not obtained.
This can be compensated by adding three degrees of freedom at the midpoints of the
edges of the triangle and increasing accordingly the size of the polynomial space.


    
    𝒫   = P_5(K), 
    ℒ(w)    = 
            w(x^k),    for 1≤ k ≤ 3,
    ∂ w/∂ x_i(x^k),    for 1 ≤ k ≤ 3
        and 1 ≤ i ≤ 2,
    ∂^2 w/∂ x_i ∂ x_j(x^k),    for 1 ≤ k ≤ 3
        and 1 ≤ i, j ≤ 2,
    ∂ w/∂ n(e^k),    for 1 ≤
                    k ≤ 3.



The Argyris triangle can be further generalized to higher-order polynomial
spaces, cf. P. Šolín <cit.>. Triangular macroelements such as
the Hsieh–Clough–Tocher triangle are not locally H^4(K) and therefore
additional jump terms are present inside the elements.  Various conforming
quadrilateral elements have been proposed in the literature for the plate
bending problem cf. Ciarlet <cit.>. The proofs of the  lower bound  that we presented
do not directly apply to quadrilateral elements,  but the  techniques can be adapted to them as well.





































§ NUMERICAL RESULTS




In our examples, we will use the fifth degree Argyris triangle. On a uniform mesh for a solution u∈ H^r (Ω), with r≥ 2,  we thus have the error estimate <cit.> 

    ‖ u-u_h‖_2 ≲ h^s | u |_r,

with s=min{r-2,4}.
Since the mesh length is related to the number of degrees of freedom N by  h∼ N^-1/2 on a uniform mesh,  we can also write

    ‖ u-u_h‖_2 ≲ N^-s/2| u |_r.

If the solution is smooth, say r≥ 6, we thus have the estimates 

    ‖ u-u_h‖_2 ≲ h^4  ‖ u-u_h‖_2 ≲ N^-2.

In fact, the  rate N^-2 is  optimal also on a general mesh since, except for a polynomial solution, it holds <cit.>

    ‖ u-u_h‖_2 ≳ N^-2.



In the adaptive  computations we use the following strategy for marking the
elements that will be refined <cit.>.


    Given a partition , error indicators η_K, K ∈ and a
    threshold θ∈ (0,1), mark K for refinement if η_K ≥θmax_K^'∈η_K^'.


The parameter θ has an effect on the portion of elements that are
marked, i.e. for θ = 0 all elements are marked and for θ = 1 only
the element with the largest error indicator value is marked. We simply take
θ = 0.5 which has proven to be a feasible choice in most cases.

The set of marked elements are refined using 
Triangle <cit.>, version 1.6, by requiring additional
vertices at the edge midpoints of the marked elements and by allowing the mesh
generator to improve mesh quality through extra vertices. The default minimum
interior angle constraint of 20 degrees is used.


The regularity of the solution depends on the regularity of the load and the corner singularities, cf. <cit.>. Below we consider two sets of problems, one where the regularity is mainly restricted by the load, and another one where the load is uniform and the corner singularities dominate.



 §.§ Square plate, Navier solution


A classical series solution to the Kirchhoff plate bending problem, the Navier
solution <cit.>, in the special case of a unit square with simply supported boundaries
and the loading

    f(x) = 
            f_0,    if x∈ [12-c, 12+c] × [12-d, 12+d],
    
            0,    otherwise,

reads

    u(x,y) = 16 f_0/D π^6∑_m=1^∞∑_n=1^∞sinmπ/2sinnπ/2sinm π csinnπ d/mn(m^2+n^2)^2sinm π xsinn π y.

In the limit c ⟶ 0 and 2 c f_0 ⟶ g_0 we get the
line load solution

    u(x,y) = 8 g_0/D π^5∑_m=1^∞∑_n=1^∞sinmπ/2sinnπ/2sinnπ d/n(m^2+n^2)^2sinm π xsinn π y,

and in the limit c, d ⟶ 0 and 4 c d f_0 ⟶ F_0 we obtain
the point load solution

    u(x,y) = 4 F_0/D π^4∑_m=1^∞∑_n=1^∞sinmπ/2sinnπ/2/(m^2+n^2)^2sinm π xsinn π y.

From the series we can infer that the solution is in  H^3-ϵ(Ω), H^7/2-ϵ(Ω) and H^9/2-ϵ(Ω), for any ϵ>0, for the point load, line load and the square load, respectively.  
in the three cases. On a uniform mesh, one should thus observe the convergence rates N^-0.5,    N^-0.75, and N^-1.25.

An unfortunate property of the series solutions is that the partial
sums converge very slowly. This makes computing the difference between
the finite element solution and the series solution in H^2(Ω) and L^2(Ω)-norms
a challenging task since the finite element solution quickly ends up
being more accurate than any reasonable partial sum. In fact, the "exact" series solution is practically useless, for example,
 for computing the shear force which is an important design parameter.
 
 
The H^2(Ω)-norm is equivalent to the energy norm,

    v= √(a(v,v)),

with which the error  is straightforward to compute.
In view of the  Galerkin
orthogonality and symmetry, one obtains
    
    u-u_h^2 = a(u-u_h,u) 
                              = l(u-u_h) ,

    i.e. the error is given by
     
    u-u_h= √(
                                l(u-u_h))  .
 
 This is especially useful for the point load for which 
  
    u-u_h= √(
                              F_0(u(12,12)-u_h(12,12)))  .
 
Evaluating the series solution at the point of maximum deflection gives <cit.> 

    u(12,12)    = 4 F_0/D π^4∑_m=1^∞∑_n=1^∞(sinmπ/2sinnπ/2)^2/(m^2+n^2)^2
       = 4 F_0/D π^4∑_m=1^∞(sinmπ/2)^2 ∑_n=1^∞(sinnπ/2)^2/(m^2+n^2)^2
       = F_0/2 D π^3∑_m=1^∞(sinmπ/2)^2(sinhmπ-mπ)/m^3(1+coshmπ).


We first consider a point load with F_0=1, d=1, E=1 and ν=0.3,
and compare the true error with the estimator η. In this case, we have the
approximate maximum displacement u(12,12) ≈ 0.1266812,  computed by evaluating and summing the first 10 million terms of the series
(<ref>).  Starting with an initial mesh shown in
Fig. <ref>, we repeatedly mark and refine the mesh to obtain
a sequence of meshes, see Fig. <ref> where the values of the
elementwise error estimators are depicted for four consecutive meshes. Note that the estimator and the adapted
marking strategy initially refine heavily in the neighborhood of the point load
as one might expect based on the regularity of the solution in the vicinity of
the point load.

In addition to the adaptive strategy, we solve the problem using a uniform mesh
family where we repeatedly split each triangle into four subtriangles starting
from the initial mesh of Fig. <ref>. The
energy norm error and η versus the number of degrees of freedom N are
plotted in Fig. <ref>.  The results show that the adaptive
meshing strategy improves  significantly  the rate of convergence in the energy
norm. In Fig. <ref>,  we have also plotted, for reference, the slopes corresponding to the expected convergence
rate O(N^-0.5) for uniform refinement and the optimal
convergence rate for P_5 elements, O(N^-2).

In Fig. <ref> is is further revealed that the energy norm
error and the estimator η follow similar trends. This is exactly what one
would expect given that the estimator is an upper and lower bound for the true
error modulo an unknown constant. This is better seen by drawing the normalized
ratio η over u-u_h, see Fig. <ref>. 
Since the estimator correctly follows the true error and an accurate
computation of norms like u-u_h_2 is expensive, the rest of the
experiments document only the values of η and N for the purpose of
giving idea of  the convergence rates.

We continue with the line load case taking g_0=1 and d=1/3, and using the same
material parameter values as before.  The initial and final meshes are shown in
Fig. <ref>. The estimator can be seen to primarly focus on
the end points of the line load. The values of η and N are visualized in Fig. <ref>
together with the expected and the optimal rates of convergence. Again the
adaptive strategy improves the convergence of the total error in comparison to the uniform
refinement strategy. The local error estimators and the adaptive process are
presented in Fig. <ref>.

We finish this subsection by solving the square load case with f_0=1,
c=d=1/3 and the same material parameters as before. The initial and the final
meshes are shown in Fig. <ref>. The convergence
rates are visualized in Fig. <ref> and the local error
estimators in Fig. <ref>. An improvement in the
convergence rate is again visible in the results.






    ndofs nelems energynorm eta
    70 8 0.0334469831818 1.03051270004
    206 32 0.0169100336631 0.493682375884
    694 128 0.00843662159162 0.247183724801
    2534 512 0.00421783761299 0.123606218404



    ndofs nelems energynorm eta
    70 8 0.0334469831818 1.03051270004
    206 32 0.02365310132 0.82181176783
    278 48 0.0119938888183 0.416323779618
    350 64 0.00610693490316 0.212543389914
    422 80 0.0032646626221 0.114554845207
    494 96 0.00199970066259 0.0715111566944
    566 112 0.00139625006813 0.0557785320735





























































    ndofs nelems eta
    219 38 0.112263251802
    774 152 0.0392580419907
    2910 608 0.0138990740805
    11286 2432 0.00491278489704



    ndofs nelems eta
    219 38 0.0874824228262
    318 60 0.0172422350444
    732 152 0.00563138145735
    1158 242 0.00175456884437
    1712 362 0.000633039818195
    2270 486 0.00041008970146
    3394 728 0.000183071024098















































    ndofs nelems eta
    208 34 0.0380160176302
    716 136 0.0080971107166
    2650 544 0.00102523606838
    10190 2176 0.000137845398212



    ndofs nelems eta
    208 34 0.0380160176302
    407 72 0.00957105089795
    452 82 0.00538652794625
    653 122 0.00227568555901
    957 188 0.00119173613498
    1453 292 0.0004456239461
    1478 296 0.000394849345563
    1907 389 0.000227748630462
    2867 593 9.98047955912e-05
    3074 639 7.70829141819e-05

















































 


 §.§ L-shaped domain


Next we solve the Kirchhoff plate problem in L-shaped domain with uniform loading f=1
and the following three sets of boundary conditions:

    
  * Simply supported on all boundaries.
    
  * Clamped on all boundaries.
    
  * Free on the edges sharing the re-entrant corner and simply supported
        along the rest of the boundary.

Due to the presence of a re-entrant corner, the solutions belong to
H^2.33(Ω), H^2.54(Ω) and H^2.64(Ω) in the cases 1,
2 and 3, respectively(see <cit.>). As before, we use fifth-order Argyris
elements to demonstrate the effectiveness of the adaptive solution strategy.
The initial and the final meshes are shown in Fig. <ref>. The resulting
total error estimators and unknown counts are visualized in Fig. <ref>.




    ndofs nelems eta
    593 111 2.02577144144
    2179 444 1.46810976277
    8348 1776 1.12381491681
    32674 7104 0.871033873118



    ndofs nelems eta
    593 111 2.02577144144
    638 121 1.86546429104
    727 140 1.59759130521
    824 160 1.23885967393
    921 180 0.969188164204
    1018 200 0.762907587751
    1115 220 0.603440209985
    1212 240 0.479294415861



    ndofs nelems eta
    593 111 0.924193320451
    2179 444 0.595889771428
    8348 1776 0.40426945804
    32674 7104 0.275514833081



    ndofs nelems eta
    593 111 0.924193320451
    638 121 0.725209767782
    727 140 0.579349194803
    824 160 0.392423628718
    921 180 0.269295860608
    1018 200 0.187398321262
    1115 220 0.133243366253
    1212 240 0.0982385833818



    ndofs nelems eta
    593 111 0.672667391309
    2179 444 0.388381670884
    8348 1776 0.229683630386
    32674 7104 0.137832309494



    ndofs nelems eta
    593 111 0.672667391309
    690 131 0.391628712819
    778 149 0.233378180614
    866 167 0.143258369597
    954 185 0.0917709133812
    1042 203 0.0633457066741
    1156 226 0.0454579296491
    1304 255 0.032930857058






















































































































 









§ ACKNOWLEDGEMENTS

The authors thank the two anonymous referees and Prof. A. Ern  for comments that improved the final version of the paper. 




siamplain
