Abstract
The symmetric C0 interior penalty method is one of the most popular discontinuous Galerkin methods for the biharmonic equation. This paper introduces an automatic local selection of the involved stability parameter in terms of the geometry of the underlying triangulation for arbitrary polynomial degrees. The proposed choice ensures a stable discretization with guaranteed discrete ellipticity constant. Numerical evidence for uniform and adaptive mesh refinement and various polynomial degrees supports the reliability and efficiency of the local parameter selection and recommends this in practice. The approach is documented in 2D for triangles, but the methodology behind can be generalized to higher dimensions, to non-uniform polynomial degrees, and to rectangular discretizations. An appendix presents the realization of our proposed parameter selection in various established finite element software packages.
1 Introduction
A classical conforming finite element method for the plate problem originates from the work of Argyris [4] with a quintic polynomial ansatz space P5(T) and 21 degrees of freedom on each triangle T. The practical application appears less prominent due to the higher computational efforts [28], although recent works [23, 33] establish a highly efficient adaptive multilevel solver for the hierarchical Argyris FEM. Due to their easier implementation, classical nonconforming FEMs [19, 25, 26, 30] and discontinuous Galerkin schemes appear advantageous compared to conforming discretizations. At least the growing number of publications on these schemes indicates a higher popularity.
Discontinuous Galerkin finite element methods [6, 16, 19, 27, 44] such as the C0 interior penalty method are a popular choice for fourth-order problems in order to avoid C1-conforming finite elements. Their main drawback is the dependence of the stability on a sufficiently large penalty parameter. The choice of this parameter is often heuristical and based on the individual experience of the user. This paper proposes an explicit geometry-dependent and local selection of the penalty parameter that guarantees the stability of the resulting scheme.
Given a source term f ∈ L2(Ω) in a bounded polygonal Lipschitz domain Ω ⊂ ℝ2, let u ∈ V :=
The Riesz representation theorem applies in the Hilbert space (V, a) and proves well-posedness of the formulation (1.1). Elliptic regularity theory verifies that f ∈ L2(Ω) implies u ∈ H2+α(Ω) ∩
For a regular triangulation 𝓣 of Ω into closed triangles and the polynomial degree k ⩾ 2, the symmetric C0 interior penalty method (C0IP) seeks
The three contributions to the bilinear form Ah : Vh × Vh → ℝ defined, for uIP, vIP ∈ Vh, by
originate from the piecewise integration by parts in the derivation of (1.1) (see [16]). While the boundedness of Ah follows from standard arguments, its coercivity is subject to the assumption of a sufficiently large penalty parameter σIP,E > 0 in the bilinear form cIP : Vh × Vh → ℝ defined, for uIP, vIP ∈ Vh, by
with the normal jumps [∇ vIP ⋅ νE]E (see [16, 24]). This paper introduces a local parameter
The penalty parameter σIP,E in (1.5) contains a prefactor a > 1. Theorem 3.1 below establishes that every choice of a > 1 leads to guaranteed stability with stability constant at least ϰ = 1 − 1/
Numerical experiments confirm the guaranteed stability and exhibit rate-optimal convergence of an adaptive C0IP for various polynomial degrees in Section 4 below. The computation of the discrete inf-sup constants (as certain eigenvalues of the discrete operator) reveals little overestimation only and recommends the proposed local parameter selection in practise. A detailed investigation of the influence of the parameter a reveals that large penalty parameters lead to a substantial increase of the condition numbers of the system matrix. Hence, we recommend a small choice of a, e.g., a = 2, in order to avoid large condition numbers. Alternatively, a strong penalization requires the application of suitable preconditioners for C0IP [11, 13, 17, 18, 20, 21].
Example 1.1
On a uniform triangulation into right isosceles triangles of the same area, the penalty parameter (1.5) for the quadratic C0IP method reads
A choice of a close to 1 leads to smaller penalty parameters for the majority of interior edges compared to values from the literature, e.g., σIP,E = 5 in [14].
The remaining parts of the paper are organized as follows. Section 2 specifies the notation such that Section 3 can present the main result of the paper in Theorem 3.1. Numerical experiments investigate the stability of the scheme with the suggested automatic penalty selection in Section 4. The coercivity constants are computed numerically and compared with the theoretically established value. Moreover, the performance of an adaptive mesh-refinement algorithm is examined. Appendix A presents the numerical realization of this automatic choice of the penalty parameter to existing FEM software packages such as the unified form language [35, 36, 39], deal.II [5], and NGSolve [42]. The implementation of the C0 interior penalty method in these packages turns out to be very compact. A complete and accessible documentation of a self-contained MATLAB code for the lowest-order discretization is available in [22, App. B].
2 Notation
Standard notation of Lebesgue and Sobolev spaces, their norms, and L2 scalar products applies throughout the paper. Instead of the space V in the weak formulation (1.1), discontinuous Galerkin methods employ the piecewise Sobolev space
with respect to a shape-regular triangulation 𝓣 of the domain Ω into closed triangles T ∈ 𝓣 with interior ∫(T) ⊂ Ω. Define the space of piecewise polynomials
of total degree at most k ∈ ℕ0. For v ∈ H2(𝓣), the piecewise application of the distributional derivatives leads to the piecewise Hessian
The piecewise Sobolev functions v ∈ H2(𝓣) allow for the evaluation of averages 〈v〉E and jumps [v]E across an edge E ∈ 𝓕. Each interior edge E ∈ 𝓕(Ω) of length hE := |E| is the common edge of exactly two triangles T+, T− ∈ T, written E = ∂ T+ ∩ ∂ T−. Then
The unit normal vector νE is oriented such that νE ⋅ νT±|E = ± 1 for the outward unit normal vectors νT± of T±. For each boundary edge E ∈ 𝓕(∂Ω), let T+ ∈ 𝓣 denote the unique triangle with edge E ∈ 𝓕(T+) and set [v]E := 〈 V 〉E := (v |T+)|E. Analogous definitions apply for vector- or matrix-valued polynomials. Let 𝓣 denote a shape regular triangulation of the polygonal Lipschitz domain Ω into closed triangles and 𝓥 (resp. 𝓥(Ω) or 𝓥(∂ Ω)) the set of all (resp. interior or boundary) vertices [10, 15]. Let 𝓕 (resp. 𝓕(Ω) or 𝓕(∂ Ω)) be the set of all (resp. interior or boundary) edges. For each triangle T ∈ 𝓣 of area |T|, let 𝓥(T) denote the set of its three vertices and 𝓕(T) the set of its three edges. Abbreviate the edge patch by ωE := ∫(T+ ∪ T−) ⊆ Ω for an interior edge E = ∂ T+ ∩ ∂ T− ∈ 𝓕(Ω) and by ωE := ∫(T+) ⊆ Ω for a boundary edge E ∈ 𝓕(T+) ∩ 𝓕(∂Ω).
The remaining contributions to the bilinear form Ah in (1.3) consist of the piecewise energy scalar product apw : Vh × Vh → ℝ and the jump term 𝓙: Vh × Vh → ℝ defined, for uIP, vIP ∈ Vh, by
Note that a different convention for the orientation of the normal vector νE in [16] leads to other signs in the bilinear form Ah compared to (1.3).
The notation A ≲ B abbreviates A ⩽ C B for a positive, generic constant C, solely depending on the domain Ω, the polynomial degree, and the shape regularity of the triangulation 𝓣; A ≈ B abbreviates A ≲ B ≲ A.
For two indices j, k ∈ ℕ, let δjk ∈ {0, 1} denote the Kronecker symbol defined by δjk := 1 if and only if j = k. The enclosing single bars |•| apply context-sensitively and denote not only the modulus of real numbers, the Euclidean norm of vectors in ℝ2, but also the cardinality of finite sets, the area of two-dimensional Lebesgue sets, and the length of edges.
3 Stability
This section develops a novel stability analysis with a penalty parameter σIP,E in the bilinear form cIP from (1.4). Let the discrete space Vh be equipped with the mesh-dependent norm ∥•∥h defined, for vIP ∈ Vh, by
with the piecewise semi-norm
Theorem 3.1
Given any a > 1, define the penalty parameter σIP,E > 0 as in (1.5). For every choice of a > 1, the constant ϰ := 1 − 1/
Four remarks are in order before the proof of the Theorem 3.1 concludes this section.
Remark 3.1
(choice of energy norm). The σ-independence of the stability constant ϰ in (3.1) justifies that ∥•∥h is a suitable choice for the energy norm on Vh. Any equivalent norm lead to a stable discretization as well. However, the stability constant ϰ might depend on σIP,E ≈ 1. For instance, the stability (3.1) remains valid for all alternative weights ensuring σ͠IP,E ⩾ σIP,E for all edges E ∈ 𝓕. In the case of a uniform triangulation into right isosceles triangles and a ⩾ 4/3, Example 1.1 shows that 1 ⩽ σIP,E and discrete stability holds for the choice of the following norm not including σIP,E
Remark 3.2
(other types of boundary conditions). The stability analysis generalizes to other types of boundary conditions. In the case of simply supported plates, i.e., u = Δ u = 0 on ∂ Ω, the bilinear form Ah from (1.3) is replaced by [12, Rem. 2]:
with the modified contributions
The proof of Theorem 3.1 applies verbatim in this case. For Cahn–Hillard boundary conditions ∂ u / ∂ n = ∂ Δ u / ∂ n = 0, the bilinear form is equal to the bilinear form in (1.3), but the space Vh consists of Pk(𝓣) functions vanishing at a single point in Ω [12, Rem. 2]. Since the proof of Theorem 3.1 does not exploit the boundary conditions of vh ∈ Vh, the estimate (3.1) holds for all piecewise polynomials vh ∈ Pk(𝓣) for k ⩾ 2. Consequently, the result generalizes to Cahn–Hillard boundary conditions as well.
Remark 3.3
(variable polynomial degrees). The assumption of uniform polynomial degree can be weakened, as follows for a different polynomial degree kT on every triangle T ∈ 𝓣. For an edge E ∈ 𝓕, an analogous argumentation leads to the choice of the edge-dependent parameter
Remark 3.4
(generalization to rectangles). Theorem 3.1 can be generalized to rectangular meshes. To this end, let 𝓣 denote a regular triangulation into closed rectangles and Qk(𝓣) the space of continuous piecewise polynomials of partial degree at most k ∈ ℕ. The sharp discrete trace inequality from [34, Eq. (C.20)] for rectangles T ∈ 𝓣 with edge E ∈ 𝓕(T) reads, for all qk ∈ Qk(T),
Then, for edge E ∈ 𝓕 and a > 1, an analog of Theorem 3.1 leads to the following choice of the edge-dependent parameter
The proof of the Theorem 3.1 employs a discrete trace inequality.
Lemma 3.1
(see [45, Thm. 3]). Every polynomial qk ∈ Pk(T) of degree at most k ∈ ℕ0 satisfies
The multiplicative constant is sharp in the sense that it cannot be replaced by any smaller constant in the absence of further conditions on the function qk.
Proof of Theorem 3.1
The point of departure is the difference, for every ϰ > 0,
The first step of the proof bounds the jump term 𝓙 with the weighted Young inequality, for every ε > 0,
Inserting this estimate into (3.3) confirms the lower bound
The second step estimates the average
In accordance with the bilinear form apw from (2.3), the L2 norm employs the Frobenius norm |•| for matrix-valued functions. A Cauchy–Schwarz inequality in ℝ2, the submultiplicativity of the Frobenius norm, and the normalization |νE| = 1 prove
for almost every x ∈ E. For an interior edge E = ∂ T+ ∩ ∂ T− ∈ 𝓕(Ω) with the neighboring triangles T+ and T−, this and the definition of the average from (2.2) show
The combination with a triangle inequality, another weighted Young inequality with γ > 0, and (3.5) results in
The multiplicative constant max {(1 + γ)/|T+|, (1 + γ−1)/|T−|} is minimal in the case of the equality (1 + γ)/|T+| = (1 + γ−1)/|T−|. Hence the optimal value γ = | T+|/| T−| leads to
The multiplication with hE/(ε σIP,E) and the definition of σIP,E in (1.5) show
For a boundary edge E ∈ 𝓕(∂ Ω) with adjacent triangle T+ ∈ 𝓣, the definition of the average, the equality |νE| = 1, and (3.5) verify
Consequently, the multiplication with hE/(ε σIP,E) and the definition of the penalty parameter σIP,E result in the analogous estimate (3.6) for E ∈ 𝓕(∂Ω).
The sum of (3.6) over all edges E ∈ 𝓕 and the finite overlap of the edge patches (ω(E): E ∈ 𝓔) lead to
The third step concludes the proof with a combination of the lower bound (3.4), and the estimate (3.7) for
Every choice of 0 < ϰ ⩽ min{1 − 1/(a ε), 1 − ε} leads to a nonnegative lower bound in (3.8), and so proves the claim (3.1); ε = 1/
4 Numerical experiments
This section investigates an implementation of the C0IP method from (1.2) using an unpublished in-house software package AFEM2 in MATLAB version 9.9.0.1718557 (R2020b) Update 6. The investigation considers four computational benchmark examples with domains displayed in Fig. 1. If not mentioned otherwise, the numerical experiments employ a = 4.
Visualization of the initial triangulations 𝓣0 for the four benchmark problems.
4.1 Adaptive mesh refinement
The experiments include a comparison of uniform and adaptive mesh refinement driven by the a posteriori error estimator η2(𝓣) := ∑T∈𝓣 η2(T) with, for all T ∈ 𝓣,
from [12, Sect. 7.1]. The software employs the Dörfler marking strategy [29] and newest-vertex bisection according to [43] in the adaptive Algorithm 1 to compute a sequence of nested triangulations (𝓣ℓ)ℓ ∈ ℕ0. The marking strategy with minimal cardinality of the set 𝓜ℓ of marked triangles can be realized in linear computational complexity [40].
Another main aspect of the numerical investigation concerns the influence of the parameter a > 1 in (1.5) to the stability of the scheme. To this end, the principal eigenvalue λ1,ℓ to the following discrete eigenvalue problem is computed: Seek
The Rayleigh–Ritz principle shows that ϰ = 1 − 1/
Algorithm 1
Adaptive C0 interior penalty method
Input: regular triangulation 𝓣0 and bulk parameter 0 < ϑ ⩽ 1.
for ℓ = 0, 1, 2, … do
Solve (1.2) with respect to triangulation 𝓣ℓ for solution
Compute refinement indicator η2(𝓣ℓ, T) from (4.1) for all T ∈ 𝓣ℓ.
Mark a minimal subset 𝓜ℓ ⊆ 𝓣ℓ by the Dörfler criterion
Refine 𝓣ℓ to 𝓣ℓ+1 by newest-vertex bisection such that 𝓜ℓ ⊆ 𝓣ℓ ∖ 𝓣ℓ+1.
end for
Ouput: sequence of triangulations 𝓣ℓ with uℓ and η(𝓣ℓ) for ℓ ∈ ℕ0.
4.2 Singular solution on L-shaped domain
This benchmark considers the L-shaped domain Ω = (−1, 1)2 ∖ [0, 1)2 with interior angle ω = 3π/2 at the reentrant corner. This determines the noncharacteristic solution α := 0.5444837 to sin2(α ω) = α2 sin2(ω) in
The exact singular solution in polar coordinates from [31, p. 107] reads
The right-hand side f := Δ2 u is computed accordingly. The reduced regularity of the exact solution u leads to the empirical convergence rate α/2 for uniform mesh refinement displayed by dashed lines in Fig. 3. The adaptive refinement strategy results in local mesh refining at the reentrant corner in Fig. 2. For all polynomial degrees k = 2, …, 5, the adaptive algorithm recovers the optimal convergence rates with respect to the number of degrees of freedom (ndof) as displayed by the solid lines in Fig. 3.
Adaptively refined mesh (ϑ = 0.5) with polynomial degree k = 2 and 10406 degrees of freedom for the L-shaped domain in Section 4.2.
Convergence history plots for the L-shaped domain in Section 4.2.
4.3 Singular solution on 1/8 cusp domain
This benchmark problem investigates the 1/8 cusp domain Ω := (−1, 1)2 ∖conv{(0, 0), (1, −1), (1, 0)} with interior angle ω = 7π/4 as depicted in Fig. 1b. The right-hand side f := Δ2 u is given by the exact singular solution in polar coordinates
analogously to (4.5) with gα,ω from (4.4) for the parameter α = 0.50500969. The singularity causes an empirical convergence rate α/2 for uniform mesh refinement in Fig. 4, while the adaptive algorithm recovers the optimal rates for the considered polynomial degrees. In the case k = 5, the highly adapted meshes lead to nearly singular matrices close to the number of 106 degrees of freedom. Undisplayed triangulation plots confirm the increased adaptive refinement towards the reentrant corner.
Convergence history plots for the 1/8 cusp domain in Section 4.3.
4.4 Dumbbell-slit domain
The benchmark problem considers the uniform force f ≡ 1 on the dumbbell-slit domain
depicted in Fig. 1c. Although Ω is no Lipschitz domain, Figures 5 confirm the optimal convergence rates with adaptive mesh refinement. Undisplayed plots of the adaptive triangulations show an increased refinement towards the slit as well as the two reentrant corners at (1, −0.75) and (3, −0.75).
Convergence history plots for the dumbbell domain in Section 4.4.
4.5 Four-slit domain
Let
denote the square domain with one slit at each edge as depicted in Fig. 1d. This benchmark also considers the uniform load f ≡ 1. Similarly to the previous benchmarks, optimal convergence rates can be observed in the case the adaptive mesh refinement for the polynomial degrees k = 2, …, 5 in Fig. 6. Undisplayed triangulation plots indicate that the adaptive algorithm detects the singularities of the unknown exact solution at tips of the four slits.
Convergence history plots for the fourslit domain in Section 4.5.
4.6 Investigation of parameter a
Given the triangulation 𝓣ℓ on the level ℓ ∈ ℕ0 from the adaptive Algorithm 1, the principle eigenvalue λ1,ℓ of the discrete eigenvalue problem (4.2) with respect to 𝓣ℓ is the stability constant in (3.1). The computation applies MATLAB’s eigs function to the eigenvalue problem of the matrices B(𝓣ℓ) and N(𝓣ℓ) from (4.3).
Figure 7 displays λ1,ℓ in dependence of the number of degrees of freedom for different choices of a > 1 using adaptive and uniform mesh refinement on the L-shaped domain in Section 4.2; a large parameter a results in a larger stability constant tending towards 1. The guaranteed lower bound ϰ = 1 − 1/
Plot of stability constant λ1,ℓ with various selections of parameter a for the L-shaped domain in Subsection 4.2.
An interesting empirical observation is that adaptive mesh refinement and higher-order polynomials even lead to slightly more stable discretizations. An over-stabilization as in the theoretical analysis from [9] is not necessary in the computational benchmarks. This supports the advantage of the automated choice of the penalization with σIP,E for guaranteed stability.
The condition number cond1(B(𝓣ℓ)) with respect to the 1-norm of the system matrix B(𝓣ℓ) from (4.3) depends on a as displayed in Fig. 8. The presented lower bounds for cond1(B(𝓣ℓ)) are computed using MATLAB’s condest function. Naturally, the number of degrees of freedom has the biggest influence and the rate of the increase depends on the polynomial degree k. While the importance of a seems negligible, the difference between the condition number for a = 1 and a = 100 is still about two orders of magnitude. Figure 9b illustrates this linear relation between the parameter a and the condition number on a fixed uniform mesh of the L-shaped domain. This suggests that a small a is advantageous.
Plot of condition number cond1(B(𝓣ℓ)) of the system matrix B(𝓣ℓ) from (4.3) with various selections of a for the L-shaped domain in Subsection 4.2.
Plots of a-dependence of quantities on a uniform mesh with 6 144 triangles of the L-shaped domain in Subsection 4.2.
4.7 Conclusions
The numerical experiments validate the theoretically established stability result with the proposed penalty parameter from (1.5) for polynomial degrees k = 2, 3, 4, 5. Figures 3–6 exhibit rate-optimal convergence of the adaptive C0 interior penalty method.
A detailed investigation of the influence of the parameter a > 1 recommends the choice of a as small as possible in order to avoid large condition numbers of the linear system of equations. We suggest a = 2 with ϰ = 1 − 1/
Funding statement: This paper has been supported by SPARC project (ID 235) ‘The Mathematics and Computation of Plates’ and by the Austrian Science Fund (FWF) standalone project P33216 ‘Computational nonlinear PDEs’.
Acknowledgment
It is our pleasure to thank Zhaonan Dong (of Inria, Paris) for fruitful discussions on the discrete trace inequality.
References
[1] S. Agmon, Lectures on Elliptic Boundary Value Problems, AMS Chelsea Publishing, Providence, RI, 2010.10.1090/chel/369Search in Google Scholar
[2] M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells, Unified form language: a domain-specific language for weak formulations and partial differential equations, ACM Trans. Math. Software 40 (2014), No. 2, Art. 9, 37.10.1145/2566630Search in Google Scholar
[3] 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, Archive of Numerical Software 3 (2015), 9–23.Search in Google Scholar
[4] J. H. Argyris, I. Fried, and D. W. Scharpf, The TUBA family of plate elements for the matrix displacement method, Aeronautical J. 72 (1968), 701–709.10.1017/S000192400008489XSearch in Google Scholar
[5] D. Arndt, W. Bangerth, M. Feder, M. Fehling, R. Gassmöller, T. Heister, L. Heltai, M. Kronbichler, M. Maier, P. Munch, J.-P. Pelteret, S. Sticko, B. Turcksin, and D. Wells, The deal.II Library, Version 9.4, J. Numer. Math. 30 (2022), No. 3, 231–246.10.1515/jnma-2022-0054Search in Google Scholar
[6] G. A. Baker, Finite element methods for elliptic equations using nonconforming elements, Math. Comp. 31 (1977), No. 137, 45–59.10.1090/S0025-5718-1977-0431742-5Search in Google Scholar
[7] P. Bastian, M. Blatt, A. Dedner, N.-A. Dreier, C. Engwer, R. Fritz, C. Grüninger, D. Kempf, R. Klöfkorn, M. Ohlberger, and O. Sander, The Dune framework: basic concepts and recent developments, Comput. Math. Appl. 81 (2021), 75–112.10.1016/j.camwa.2020.06.007Search in Google Scholar
[8] H. Blum and R. Rannacher, On the boundary value problem of the biharmonic operator on domains with angular corners, Math. Methods Appl. Sci. 2 (1980), No. 4, 556–581.10.1002/mma.1670020416Search in Google Scholar
[9] A. Bonito and R. H. Nochetto, Quasi-optimal convergence rate of an adaptive discontinuous Galerkin method, SIAM J. Numer. Anal. 48 (2010), No. 2, 734–771.10.1137/08072838XSearch in Google Scholar
[10] D. Braess, Finite Elements, 3rd ed., Cambridge University Press, Cambridge, 2007.Search in Google Scholar
[11] S. C. Brenner, J. Cui, T. Gudi, and L.-Y. Sung, Multigrid algorithms for symmetric discontinuous Galerkin methods on graded meshes, Numer. Math. 119 (2011), No. 1, 21–47.10.1007/s00211-011-0379-ySearch in Google Scholar
[12] S. C. Brenner, C0 interior penalty methods, In: Frontiers in Numerical Analysis — Durham 2010, Lect. Notes Comput. Sci. Vol. 85, Springer, Heidelberg, 2012, pp. 79–147.10.1007/978-3-642-23914-4_2Search in Google Scholar
[13] S. C. Brenner, C. B. Davis, and L.-Y. Sung, Additive Schwarz preconditioners for the obstacle problem of clamped Kirchhoff plates, Electron. Trans. Numer. Anal. 49 (2018), 274–290.10.1553/etna_vol49s274Search in Google Scholar
[14] S. C. Brenner, T. Gudi, and L.-Y. Sung, An a posteriori error estimator for a quadratic C0-interior penalty method for the biharmonic problem, IMA J. Numer. Anal. 30 (2010), No. 3, 777–798.10.1093/imanum/drn057Search in Google Scholar
[15] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, 3rd ed., Texts in Applied Mathematics, Vol. 15, Springer, New York, 2008.10.1007/978-0-387-75934-0Search in Google Scholar
[16] S. C. Brenner and L.-Y. Sung, C0 interior penalty methods for fourth order elliptic boundary value problems on polygonal domains, J. Sci. Comput. 22-23 (2005), 83–118.10.1007/s10915-004-4135-7Search in Google Scholar
[17] S. C. Brenner and L.-Y. Sung, Multigrid algorithms for C0 interior penalty methods, SIAM J. Numer. Anal. 44 (2006), No. 1, 199–223.10.1137/040611835Search in Google Scholar
[18] S. C. Brenner, L.-Y. Sung, and K. Wang, Additive Schwarz preconditioners for C0 interior penalty methods for the obstacle problem of clamped Kirchhoff plates, Numer. Methods Partial Differential Equations 38 (2022), No. 1, 102–117.10.1002/num.22834Search in Google Scholar
[19] S. C. Brenner, L.-Y. Sung, H. Zhang, and Y. Zhang, A Morley finite element method for the displacement obstacle problem of clamped Kirchhoff plates, J. Comput. Appl. Math. 254 (2013), 31–42.10.1016/j.cam.2013.02.028Search in Google Scholar
[20] S. C. Brenner and K. Wang, Two-level additive Schwarz preconditioners for C0 interior penalty methods, Numer. Math. 102 (2005), No. 2, 231–255.10.1007/s00211-005-0641-2Search in Google Scholar
[21] S. C. Brenner and J. Zhao, Convergence of multigrid algorithms for interior penalty methods, Appl. Numer. Anal. Comput. Math. 2 (2005), No. 1, 3–18.10.1002/anac.200410019Search in Google Scholar
[22] P. Bringmann, C. Carstensen, and J. Streitberger, Local parameter selection in the C0 interior penalty method for the biharmonic equation, arXiv:2209.05221, 2023.10.1515/jnma-2023-0028Search in Google Scholar
[23] C. Carstensen and J. Hu, Hierarchical Argyris finite element method for adaptive and multigrid algorithms, Comput. Methods Appl. Math. 21 (2021), No. 3, 529–556.10.1515/cmam-2021-0083Search in Google Scholar
[24] C. Carstensen, G. Mallik, and N. Nataraj, A priori and a posteriori error control of discontinuous Galerkin finite element methods for the von Kármán equations, IMA J. Numer. Anal. 39 (2019), No. 1, 167–200.10.1093/imanum/dry003Search in Google Scholar
[25] C. Carstensen and N. Nataraj, Adaptive Morley FEM for the von Kármán equations with optimal convergence rates, SIAM J. Numer. Anal. 59 (2021), No. 2, 696–719.10.1137/20M1335613Search in Google Scholar
[26] C. Carstensen and N. Nataraj, A priori and a posteriori error analysis of the Crouzeix–Raviart and Morley FEM with original and modified right-hand sides, Comput. Methods Appl. Math. 21 (2021), No. 2, 289–315.10.1515/cmam-2021-0029Search in Google Scholar
[27] D. A. Di Pietro and A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, Vol. 69, Springer, 2011.10.1007/978-3-642-22980-0Search in Google Scholar
[28] V. Dominguez and F.-J. Sayas, Algorithm 884: A Simple Matlab implementation of the Argyris element, ACM Trans. Math. Softw. 35 (2008), No. 2, 1–11.10.1145/1377612.1377620Search in Google Scholar
[29] W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (1996), No. 3, 1106–1124.10.1137/0733054Search in Google Scholar
[30] D. Gallistl, Morley finite element method for the eigenvalues of the biharmonic operator, IMA J. Numer. Anal. 35 (2015), No. 4, 1779–1811.10.1093/imanum/dru054Search in Google Scholar
[31] D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, Classics in Mathematics, Springer-Verlag, Berlin, 2001.10.1007/978-3-642-61798-0Search in Google Scholar
[32] P. Grisvard, Singularities in Boundary Value Problems, Recherches en Mathématiques Appliquées, Vol. 22, Masson, Paris; Springer-Verlag, Berlin, 1992.Search in Google Scholar
[33] B. Gräßle, Optimal multilevel adaptive FEM for the Argyris element, Comput. Methods Appl. Mech. Engrg. 399 (2022), 115352.10.1016/j.cma.2022.115352Search in Google Scholar
[34] K. Hillewaert, Development of the discontinuous Galerkin method for high-resultion, large scale CFD and acoustics in industrial geometries, Ph.D. thesis, Université Catholique de Louvain, 2013.Search in Google Scholar
[35] R. C. Kirby, A general approach to transforming finite elements, SMAI J. Comput. Math. 4 (2018), 197–224.10.5802/smai-jcm.33Search in Google Scholar
[36] R. C. Kirby and L. Mitchell, Code generation for generally mapped finite elements, ACM Trans. Math. Softw. 45 (2019), No. 4, 1–23.10.1145/3361745Search in Google Scholar
[37] J. Nečas, Direct Methods in the Theory of Elliptic Equations, Springer Monographs in Mathematics, Springer, Heidelberg, 2012.10.1007/978-3-642-10455-8Search in Google Scholar
[38] K. B. Ølgaard, A. Logg, and G. N. Wells, Automated code generation for discontinuous Galerkin methods, SIAM J. Sci. Comput. 31 (2008/09), No. 2, 849–864.Search in Google Scholar
[39] K. B. Ølgaard, A. Logg, and G. N. Wells, Automated code generation for discontinuous Galerkin methods, SIAM J. Sci. Comput. 31 (2009), No. 2, 849–864.10.1137/070710032Search in Google Scholar
[40] C.-M. Pfeiler and D. Praetorius, Dörfler marking with minimal cardinality is a linear complexity problem, Math. Comp. 89 (2020), No. 326, 2735–2752.10.1090/mcom/3553Search in Google Scholar
[41] F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. T. McRae, G.-T. Bercea, G. R. Markall, and P. H. J. Kelly, Firedrake: automating the finite element method by composing abstractions, ACM Trans. Math. Software 43 (2016), No. 3, 1–27.10.1145/2998441Search in Google Scholar
[42] J. Schöberl, C++11 implementation of finite elements in NGSolve, Institute for Analysis and Scientific Computing, TU Wien, 30 (2014).Search in Google Scholar
[43] R. Stevenson, The completion of locally refined simplicial partitions created by bisection, Math. Comp. 77 (2008), No. 261, 227–241.10.1090/S0025-5718-07-01959-XSearch in Google Scholar
[44] E. Süli and I. Mozolevski, hp-version interior penalty DGFEMs for the biharmonic equation, Comput. Math. Appl. 196 (2007), No. 13-16, 1851–1863.10.1016/j.cma.2006.06.014Search in Google Scholar
[45] T. Warburton and J. S. Hesthaven, On the constants in hp-finite element trace inverse inequalities, Comput. Methods Appl. Mech. Engrg. 192 (2003), No. 25, 2765–2773.10.1016/S0045-7825(03)00294-9Search in Google Scholar
A Remarks on the numerical realization
The following three subsections illustrate the numerical realization of the C0IP method with the proposed local parameter selection in established software packages such as FEniCS (using the unified form language), deal.II, and NGSolve.
A.1 Realization using the unified form language (UFL)
The unified form language [2] provides a standardized way of implementing variational formulations for PDEs. It supports discontinuous Galerkin discretizations such as the C0IP method [38].
The first step defines the shape of the element domain and several geometric quantities required for the bilinear form Ah:
cell = triangle
n = FacetNormal(cell)
h = FacetArea(cell)
vol = CellVolume(cell)
Given a polynomial degree k ∈ ℕ with k ⩾ 2, these variables allow to define the penalty parameter σIP,E according to the formula (1.5):
a = 4.0
sigma = 3.0 * a * k * (k − 1)/8.0 * h(“+”)**2 * avg(1/vol)
sigma_boundary = 3.0 * a * k * (k − 1) * h**2/vol
The following lines specify the Lagrange finite element and the trial and test functions:
element = FiniteElement(“Lagrange”, cell, k)
u = TrialFunction(element)
v = TestFunction(element)
Abbreviate the Hessian matrix of finite element functions via
def hess(v):
return grad(grad(v))
The realization of the bilinear form Ah from (1.3) reads
A = inner(hess(u), hess(v)) * dx \
- inner(avg(dot(hess(u) * n, n)), jump(grad(v), n)) * dS \
- inner(dot(grad(u), n), dot(hess(v) * n, n)) * ds \
- inner(jump(grad(u), n), avg(dot(hess(v) * n, n))) * dS \
- inner(dot(hess(u) * n, n), dot(grad(v), n)) * ds \
+ sigma/h(“+”) * inner(jump(grad(u), n), jump(grad(v), n)) * dS \
+ sigma_boundary/h * inner(dot(grad(u), n), dot(grad(v), n)) * ds
Finally, the right-hand side may be defined as follows:
f = Coefficient(element)
F = f * v * dx
Many software packages employ the UFL such as DUNE [7], FEniCS [3], and Firedrake [41]. With slight adaptions of the above UFL code, the C0IP method with automatic penalty selection can be realized within these three software frameworks for arbitrary polynomial degree.
A.2 Realization in deal.II
The C++ software package deal.II [5] (version 9.4) provides a finite element implementation for rectangular grids. An example implementation of the C0IP method in the file examples/step-47.cc employs the following heuristic choice of the penalty parameter (adapted to the notation in this paper) for an edge E ∈ 𝓕 with d± := 2|T±|/hE:
This leads to a slight over-penalization compared to the parameter suggested in (3.2):
The application of the automatic parameter selection with guaranteed stability solely requires the replacement of the lines 512–518 in the file step-47.cc by
const double a = 4.0;
const double sigma_over_h =
2 * a * pow(k − 1, 2.0) * (
1/cell- > extent_in_direction(
GeometryInfo < dim > ::unit_normal_direction[f]) +
1/ncell- > extent_in_direction(
GeometryInfo < dim > ::unit_normal_direction[nf]));
for the parameter on the interior edges and the lines 622–625 by
const double a = 2.0;
const double sigma_over_h =
2 * a * pow(k - 1, 2) /
cell- > extent_in_direction(
GeometryInfo <dim> ::unit_normal_direction[face_no]);
for the parameter on the boundary edges.
A.3 Realization in Netgen/NGSolve
The C++ software package NGSolve [42] with an easy-to-use python interface provides an implementation of a hybridized C0IP method. Straight-forward modifications of the existing code enable the C0IP method presented in this paper. Given an object mesh of the Mesh class representing a regular triangulation into simplices, the Lagrange finite element of order k ∈ ℕ with 2 ⩽ k can be defined as follows:
k = 2
fes = H1(mesh, order = k, dirichlet = “left|right|bottom|top”, dgjumps = True)
This allows to define variables for trial and test functions:
u = fes.TrialFunction()
v = fes.TestFunction()
For an edge E ∈ 𝓕 with adjacent triangles T± ∈ 𝓣, recall the notation d± := 2 |T±|/hE for the altitude of T±. This quantity and the normal vector νE can be defined by special coefficient functions in NGsolve:
d = specialcf.mesh_size
n = specialcf.normal(2)
Thus, the automatic choice of the penalty parameter σIP,E reads, for interior and boundary edges,
a = 4
sigma_over_h = 3 * a * k * (k - 1)/4 * (1/d + 1/ d.Other())
sigma_over_h_bdr = 3 * a * k * (k - 1)/2 * 1/d
Abbreviating the Hessian of a function v
def hesse(v):
return v.Operator(“hesse”)
allows for the following variables for the average of the binormal part of the Hessian in the jump term 𝓙 from (2.3):
mean_d2udn2 = 0.5*InnerProduct(n, (hesse(u)+hesse(u.Other()))*n)
mean_d2vdn2 = 0.5*InnerProduct(n, (hesse(u)+hesse(u.Other()))*n)
This and the definition of the normal jump
def jumpdn(v):
return n*(grad(v)-grad(v.Other()))
lead to the bilinear form Ah from (1.3) via
A = BilinearForm(fes)
A + = SymbolicBFI(InnerProduct (hesse(u), hesse(v)) )
A + = SymbolicBFI(-mean_d2udn2 * jumpdn(v), skeleton = True)
A + = SymbolicBFI(-mean_d2vdn2 * jumpdn(u), skeleton = True)
A + = SymbolicBFI(sigma_over_h * jumpdn(u) * jumpdn(v), VOL, skeleton = True)
A + = SymbolicBFI(sigma_over_h_bdr * jumpdn(u) * jumpdn(v), BND, skeleton = True)
© 2024 Philipp Bringmann, Carsten Carstensen, and Julian Streitberger, published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.
