{'fraction_non_alphanumeric': 0.06511556867936055, 'fraction_numerical': 0.022197332247225334, 'mean_word_length': 3.7921444254696266, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Three families of grad-div-conforming finite elements

Several smooth finite element de Rham complexes are constructed in three-dimensional space, which yield three families of grad-div conforming finite elements. The simplest element has only 8 degrees of freedom (DOFs) for a tetrahedron and 14 DOFs for a cuboid. These elements naturally lead to conforming approximations to quad-div problems. Numerical experiments for each family validate the correctness and efficiency of the elements for solving the quad-div problem.


Introduction
We are concerned in this paper with the grad-div or H(grad div) finite elements used for conforming discretizations of problems involving a quad-div operator (sometimes referred to as a fourth-order div operator). The quad-div operator appears in linear elasticity [1,21,22], where the integration of (∇(∇ · u)) 2 represents the shear strain energy with the displacement of the elasticity body u. Moreover, the quad-div operator can be written as are two well-known fundamental fourth-order operators of this formulation, which have been studied extensively. Many numerical methods have been proposed for problems involving those two fourth-order operators. We refer to [3,20,[29][30][31][32] for conforming finite element methods and [8,10,11,19,24,[26][27][28][33][34][35] for other methods. However, unlike the biharmonic operator and the quad-curl operator, very limited work has been done for problems involving the quad-div operator [16]. In this paper, we will propose graddiv conforming elements which can lead to conforming approximations of the quad-div problem.
We apply discrete de Rham complexes to construct three families of grad-div conforming finite elements in three space dimensions (3D). The discrete de Rham complex with appropriate smoothness has been an important and useful tool in designing finite elements and analyzing numerical schemes, c.f., [4-6, 12, 18, 25]. Based on the de Rham complex with minimal smoothness, various well-known finite elements for computational electromagnetism or diffusion problems have been arranged in the finite element periodic table [7]. Motivated by problems in fluid and solid mechanics, there is a growing interest in constructing finite element de Rham complexes with enhanced smoothness, sometimes also referred to as Stokes complexes [13,15,20]. In this paper, for the conforming discretization of the quad-div problem, we will consider another variant of the de Rham complex, i.e., where Ω is a bounded Lipschitz domain in R 3 . The precise definition of the notations used in the above complex is given in the subsequent sections. For simplicity of presentation, throughout this paper, we will assume that Ω is contractible. Then the exactness of (1.1) follows from standard results in, e.g., [4]. The two dimensional (2D) version of the ∇× -H(grad div; Ω) where ∇× = (∂ x 2 u, −∂ x 1 u) T . If we rotate the complex (1.2) by π 2 , we will get the following complex with H(curl 2 ; Ω) := {u ∈ L 2 (Ω) : ∇ × u ∈ L 2 (Ω) and ∇ × ∇ × u ∈ L 2 (Ω)}, which has been studied in [20]. Therefore in this paper, we will only focus on the complex (1.1) to construct 3D grad-div elements.
Our new finite elements fit into a subcomplex of (1.1): The tool of discrete complexes makes it possible to reach the goal of this paper, i.e., constructing grad-curl conforming elements. Moreover, by this tool, we also fit the quaddiv problem and its conforming finite element approximations into the framework of the finite element exterior calculus (FEEC) [4,5] and hence enable the use of various tools from FEEC for the numerical analysis. For instance, we construct interpolation operators that commute with the differential operators. After that, the convergence result follows from a standard argument.
To validate the newly proposed elements, we fit them in the theoretical framework of the conforming finite element method and obtain the approximation property of the numerical solution. Moreover, we carry out a numerical experiment, which validates our theoretical results.
The remaining part of the paper is organized as follows. In Section 2, we present notations and some basic facts from homological algebra. In Section 3, we define shape functions which form local exact sequences and prove their properties. In Section 4, we construct three families of grad-curl conforming finite elements on tetrahedra and develop some theoretical results for them. In Section 5, we introduce the counterparts for cuboids.
In Section 6, we present a regularity estimate for the quad-div problem and apply the new elements to the problem. In Section 7, we provide numerical examples to verify the correctness and efficiency of our methods. Finally, concluding remarks and future work are given in Section 8.

2.1.
Notations. Unless otherwise specified, throughout the paper we assume Ω ∈ R 3 is a contractible Lipschitz domain. We adopt conventional notations for Sobolev spaces such as H m (D) or H m 0 (D) on a contractible sub-domain D ⊂ Ω furnished with the norm · m,D and the semi-norm |·| m,D . In the case of m = 0, the space H 0 (D) coincides with L 2 (D) which is equipped with the inner product (·, ·) D and the norm · D . When D = Ω, we drop the subscript D. We use H m (D) and L 2 (D) to denote the vector-valued Sobolev spaces (H m (D)) 3 and (L 2 (D)) 3 .
In addition to the standard Sobolev spaces, we also define For a subdomain D, a face f , or an edge e, we use P k to represent the space of polynomials with degree at most k, andP k , the space of homogenous polynomials of degree k. The corresponding sets of vector polynomials are denoted as P k = (P k (D)) 3 and P k = (P k (D)) 3 , respectively. We also define We use Q i,j,k (D) to denote the polynomials with three variables (x 1 , x 2 , x 3 ) where the maximal degree is i in x 1 , j in x 2 , and k in x 3 . For simplicity, we drop the subscripts i and j when i = j = k. Similarly, we use Q i,j (f ) to denote such polynomial spaces in 2D.
Let T h be a partition of the domain Ω consisting of tetrahedra or cuboids. We denote h K as the diameter of an element K ∈ T h and h as the mesh size of T h . We adopt the (1, 0, 1), (0, 1, 1), and (1, 1, 1)): where the affine mapping By a simple computation, we have We use C to denote a generic positive h-independent constant.

2.2.
Basic facts from homological algebra. We review some basic facts from homological algebra. For further details, we refer to, for example, [4]. A differential complex is a sequence of spaces V i and operators d i : condition for the exactness of (2.5) is the following dimension condition:

Local spaces and polynomial complexes
To define a finite element space, we must supply, for each element K ∈ T h , • shape functions; • DOFs to guarantee appropriate continuity.
In this section, we will specify shape functions for each space involved in the complex (1.4).
The local complex of function spaces on each K ∈ T h for (1.4) is denoted as follows: In addition to the complex (1.4), we introduce the complex (3.1) with two parameters r and k to specify degrees of spaces, which lead to several versions of complexes.
We let Σ r h (K) be P r (K) for a tetrahedral element or Q r (K) for a cuboid element, and let V r h (K) be R r (K) for a tetrahedral element or Q r−1,r,r (K) × Q r,r−1,r (K) × Q r,r,r−1 (K) for a cuboid element. For a tetrahedral element K, we set 3 is an operator which maps a scalar function to a vector field: As a special case of the Poincaré operators (c.f. [13,18]), p has the following properties: • the null-homotopy identity The above two properties can be obtained via simple calculations. By the null-homotopy identity (3.3), the right hand side of (3.2) is a direct sum.
. This shows that (3.1) is a complex. It remains to show the exactness. It's easy to check the exactness at Σ r h (K) and V r h (K). We now show that, for To prove the exactness at , we only need to show the div operator ∇· : In the following lemma, we show that W r−1,k h (K) contains some polynomial subspaces.
It plays an essential role in analyzing the approximation properties of the finite element Proof. We claim that In fact, by the polynomial preserving property of p, To show (3.4), it remains to show the two sides of (3.4) have the same dimension. By the null-homotopy identity (3.3), the right hand side is a direct sum. Therefore, which is exactly the dimension of P r−1 (K). Combining (3.4) and the fact that P r−2 (K) ⊆ We are now ready to construct grad-div conforming finite elements and complexes. We focus on analyzing tetrahedral elements and only present the definition of cuboid elements in Section 5.

Three families of grad-div conforming elements on tetrahedra
The global discrete complex with specified degree for each space is given by In this section, we construct grad-div conforming finite elements and complexes on tetrahedra. Assigning r = k − 1, k, and k + 1 in (4.1) leads to three versions of grad-div The DOFs for the Lagrange element Σ r h can be given as follows.
• Vertex DOFs M v (u): : uvds for all v ∈ P r−2 (e i ) and for all edges e i .
• Face DOFs M f (u): for all faces f i .
• Interior DOFs M K (u): The DOFs for Σ +,k−1 h can be given similarly, with only one additional interior integration DOF on K to deal with the interior bubble function. We denoteπ h as the H 1 interpolation operator to Σ +,k−1 h by these DOFs.
We choose the space V r h as the first family of Nédélec elements, which has the following DOFs: • Edge DOFs M e (u) (with a unit tangential vector τ i ): and for all edges e i .
• Face DOFs M f (u) (with a unit normal vector n i ): and for all faces f i .
• Interior DOFs M K (u): and for all elements K i .
By counting the number of DOFs, we find the DOF set has the same dimension. Then it suffices to show that if all the DOFs vanish on a function u, then u = 0. To see this, we first observe that ∇ · u = 0 by the unisolvence of the DOFs of Σ +,k−1 h (K). Then . By the face DOFs of W r−1,k h (K), u · n = 0 on faces. By integration by parts, we have which, together with the interior DOFs, leads to (u, q) K = 0 for any q ∈ P r−2 (K). (4.8) Since u · n = 0 on faces,û = (x 1φ1 ,x 2φ2 ,x 3φ3 ) T . Note that u andû are related by (2.1).
Provided u ∈ H 1/2+δ (Ω), and ∇ · u ∈ H 3/2+δ (Ω) with δ > 0 (see Lemma 4.1), we can define an H(grad div) interpolation operator i h whose restriction on K is denoted as i K and defined by Actually, this follows from the exactness of the standard finite element differential forms (e.g., [4] where V, E, F, and K denote the number of vertices, edges, faces, and 3D cells, respectively. The dimension count of the space V r h reads: From the DOFs (4.4) -(4.7), From the above dimension count, we have where we have used Euler's formula V − E + F − K = 1. This completes the proof.
We summarize the interpolations defined in Section 4.1 in the following diagram: Here Σ, V , W are three subspaces of H 1 (Ω), H(curl; Ω), and H(grad div; Ω) in which π h (orπ h ), r h , and i h are well-defined.
Now we show that the interpolations in (4.11) commute with the differential operators.
In addition to Lemma 3.2, this result also plays a key role in the error analysis below for the interpolations. ∇ · i h u =π h ∇ · u for all u ∈ W. (4.14) Proof. The proof of (4.12) can be found in [23,Chapter 5]. A similar trick can be used to prove (4.13) and (4.14). For simplicity of presentation, we omit it.
The following lemma relates the interpolation on K to that onK. the second set of (4.6) are equivalent with those forû onK. In addition, all the DOFs in (4.4), (4.5), and the first set of (4.6) are differed from those forû onK by a factor 1 det(B K ) . This means the DOFs for defining i K u and those for defining iKû are differed by a constant factor. According to Proposition 3.4.7 in [9], we complete the proof.
Next, we establish the approximation property of the interpolation operator.
Theorem 4.2. For r = k − 1, k, or k + 1, if u ∈ H s+(r−k) (Ω) and ∇ · u ∈ H s (Ω), 3/2+δ ≤ s ≤ k with δ > 0, then we have the following error estimates for the interpolation Proof. We only prove the results for integer s to avoid the technical complications. To prove (4.15), we first apply the transformation (2.1) and Lemma 4.5 to derive, for a general From Lemma 3.2, P r−1 (K) ⊆ W r−1,k h (K). Therefore, we have iKp =p whenp ∈ P r−1 (K). We obtain, with the help of Lemma 4.1 and Theorem 5.5 in [23], Mapping back to the general element K leads to Summing the above inequality (4.18) over K ∈ T h , we obtain (4.15).
As for (4.16) and (4.17), we use Lemma 4.4 and Lemma 3.2 to obtain, for i = 0, 1, where we have used the standard approximation property of the Lagrange finite element.

Three families of grad-div conforming elements on cuboids
In this section, we will construct three families of grad-div conforming finite element spaces W r−1,k for a cuboid mesh by taking r = k − 1, k, and k + 1 in (4.1).
Proceeding as in the tetrahedral case, we choose the Lagrange element space of order r on cuboids for Σ r h . The space Σ +,k−1 h is the (k − 1)-th order Lagrange element space, which is enriched with an interior bubble function when k = 2. In addition, we choose the first Nédélec element space of order r for V r h . The DOFs for these spaces can be found in [23,Chapter 6]. In the following, we will define DOFs for W r−1,k . The three grad-div elements with r = k − 1, r = k, and r = k + 1 with k = 2 for a cuboid element are shown in Fig. (5.1) • Vertex DOFs M v (u) at all vertices v i of K: • Edge DOFs M e (u) at all edges e i of K: M e (u) = 1 length(e i ) e i ∇ · uqds, ∀q ∈ P k−3 (e i ), i = 1, 2, · · · , 12 .
• Face DOFs M f (u) at all faces f i of K (with the unit normal vector n i ): • Interior DOFs M K (u):  The same theoretical results as the tetrahedral elements can be obtained by a similar argument. We omit them for this case.

The application to quad-div problems
In this section, we use the H(grad div)-conforming finite elements to solve the quad-div problem which is stated as follows: For f ∈ H(curl; Ω) and g ∈ L 2 (Ω), find u, such that u · n = 0 on ∂Ω, ∇ · u = 0 on ∂Ω.
Here, to make the problem consistent, g = ∇ × f . By taking curl on both sides of the first equation of (6.1), we see that the condition ∇ × u = g holds automatically.
The variational formulation is to seek u ∈ H 0 (grad div; Ω) such that It follows from Lax-Milgram Lemma that (6.2) is well-posedness. Taking v = ∇ × ψ with ψ ∈ H 0 (curl; Ω) in (6.2) leads to which implies ∇ × u = g holds in the sense of L 2 (Ω). Since the regularity of the solution plays a crucial role in the error analysis, we will first derive a regularity result for the quad-div problem before proceeding further. In particular, when Ω is convex, α 0 is at least 1.
Remark 6.1. With Theorem 6.1, the interpolation i h u is well-defind.
We now present the finite element scheme. We define the finite element space with vanishing boundary conditions The H(grad div)-conforming finite element method reads: seek u h ∈ W 0 h , such that The following approximation property of u h follows immediately from Céa's lemma and the duality argument.
Theorem 6.2. For r = k − 1, r = k, or r = k + 1, if u ∈ H s+(r−k) (Ω) and ∇ · u ∈ H s (Ω), 3/2 + δ ≤ s ≤ k with δ > 0, then we have the following error estimates for the numerical When Ω is a Lipschitz polyhedron, even with the lowest regularity, the scheme still has a convergence order 1/2 in H(grad div) norm.

Numerical Experiments
We now turn to a concrete example to test our new elements. We consider the problem (6.1) on a unit cube Ω = (0, 1) × (0, 1) × (0, 1) with an exact solution Then, by a simple calculation, the source term f can be derived. Note that in this case g = 0. We denote the finite element solution as u h . To measure the error between the exact solution and the finite element solution, we also denote Example 1. In this example, we test the tetrahedral elements. To this end, we partition the unit cube into N 3 small cubes and then partition each small cube into 6 congruent tetrahedra. We use the lowest-order elements in three families to solve the problem (6.1) on the uniform tetrahedral mesh. We observe from the tables that the numerical solution converges to the exact solution with a convergence order 1 for the family r = k − 1, 2 for the family r = k, and 2 for the family r = k + 1 in the sense of L 2 -norm. In addition, the three families have the same convergence order 2 in the H(div)-norm and 1 in the H(grad div)-norm, respectively. All the results coincide with Theorem 6.2, which confirms the correctness of the elements and their properties. as the weights and nodes of Legendre-Gauss quadrature rule of an order p. We also denote {w l n } p n=1 and {l n } p n=1 as the weights and nodes of Legendre-Gauss-Lobbato quadrature rule of an order p. For u = (u 1 , u 2 , u 3 ) T , we define three discrete norms. and where (x K c , y K c , z K c ) is the center of element K and 2h K x , 2h K y , 2h K z are the lengths of edges parallel to x, y, z axes, respectively. Tables 7.4, 7.5, and 7.6 shows errors measured in various norms for the lowest-order cuboid elements in the three familes. We also depict error curves with a log-log scale in Fig.   7.1. From Fig. 7.1 (A), we can observe superconvergence phenomena that |||∇(∇ · e h )||| W and |||e h ||| V converge to 0 with one order higher than ∇(∇ · e h ) and e h . In addition, from Fig. 7.1 (B)(C), we can observe superconvergence of |||∇(∇ · e h )||| W .
We can not observe any superconvergence of |||e h ||| V for r = 2, 3 and |||∇ · e h ||| U for all the 3 families when k = 2. To further investigate the superconvergence of ∇ · e h , we test the third-order (k = 3) element in the family of r = k − 1. The results are shown in Table   7.7 and Fig. 7.1(D). In this case, we can observe superconvergence of ∇ · e h .
Using these superconvergent results, together with some recovery techniques, we can construct a solution with higher accuracy if needed, which is one of the reasons that we explore the superconvergence of cuboid elements.
We conclude this section by pointing out that the three families of elements bear their own advantages. The family of r = k − 1 can be the best choice if we pursue a low computational cost, while the family with r = k + 1 stands out for its higher accuracy in the L 2 -norm without any recovery techniques.   1/22 6.5478e-08 6.5476e-08 6.6200e-07 4.8351e-07 5.5322e-05 6.5117e-06

Conclusion
In this paper, we constructed conforming finite element de Rham complexes with enhanced smoothness. This naturally leads to grad-div conforming elements in 3D which can be utilized to solve the quad-div problem. The simplest elements in our construction have only 8 and 14 DOFs for a tetrahedron and a cuboid respectively, which makes commercial adoption of the elements feasible.
Since the quad-div problem has not been extensively studied both in mathematical theory and numerical methods, there are still some mysteries about the regularity and the discretization of this problem. However, by the de Rham complex, we relate the grad-div conforming elements to the FEEC. This allows further systematic developments of the new elements and the quad-div problem. We believe that the framework may shed new light into studying this problem.
In the future, we will apply the newly proposed elements to solve practical problems and further investigate the superconvergence phenomena. Moreover, we will construct simple finite element subcomplexes for the following variant of de Rham complex: where Φ(Ω) := {v ∈ L 2 (Ω), ∇ × v ∈ H 1 (Ω)}.