A Differentiable Mapping of Mesh Cells Based on Finite Elements on Quadrilateral and Hexahedral Meshes

2.1 Mapped Shape Functions

Shape functions on a single cell K are defined by pull-back of functions on K ^ such that for 𝐱 = Φ K ⁢ ( 𝐱 ^ ) there holds

where ( ∇ ^ ⁢ Φ K ) - T is the transpose of the inverse of ∇ ^ ⁢ Φ K . For later use, we introduce the Jacobi determinant J K = det ⁡ ∇ ^ ⁢ Φ K . The implicit assumption here is the invertibility of ∇ ^ ⁢ Φ K . More so, in order to obtain uniform approximation results, we require uniform boundedness of the singular values of ∇ ^ ⁢ Φ K from below and above. A straightforward geometrical condition based on inscribed circles and diameter exists only for affinely mapped simplices. Already for bilinear quadrilaterals, we have to combine a maximum angle condition with relative edge lengths; for higher order mappings, the situation is even harder. While we do not provide proofs for such uniform bounds, the constructions in Sections 3.2 and 3.3 serve to avoid deterioration of such bounds. We refer to Figure 4 for mappings which violate such a condition. Furthermore, we will abbreviate Φ K by Φ if no confusion arises. We denote the components of Φ by φ i . Tangential and normal vectors to the boundary of K ^ are transformed to the corresponding ones on K by contravariant and covariant transformation, respectively:

(2.2)

Continuity of function values can be achieved by standard techniques of mapped finite elements. Here, we have to obtain continuity of first derivatives by additional measures. First investigating a vertex at position 𝐱 , we see that for every two cells K 1 and K 2 sharing a vertex, we obtain the additional condition

This condition can be used to establish conforming degrees of freedom on a mesh 𝕄 . A simple choice is to require ∇ ^ ⁢ Φ 1 = ∇ ^ ⁢ Φ 2 and thus decouple C 1 -continuity of the mapped element into C 1 -continuity of the mapping and C 1 -continuity of the finite element with respect to the reference geometry.

For a face of codimension one shared by cell K 1 and K 2 , all tangential derivatives of u 1 and u 2 coincide just by continuity of the function itself and the mapping. Thus, we have to ensure continuity of the normal derivative along this face. By Taylor expansion, a necessary condition is the continuity of ∂ 𝐭 ⁡ ∂ 𝐧 ⁡ u ⁢ ( p ) for any tangential vector 𝐭 of the face. By continuity, this implies that u must be C 2 at every vertex.

2.2 Hermitian Tensor Product Elements in Arbitrary Dimension

For the standard d-linear mapping, the normal derivative of Φ at an interface is determined by the vertices not adjacent to this interface, namely vertices not shared by the two cells at the interface. Therefore, continuity of this normal derivative can only be achieved using nonlocal information on a cell. Therefore, we use Hermitian interpolation and assign suitable values to the normal derivative in the vertices. The choice of such values can be made differently, depending on the structure of the mesh.

The shape function space we are going to use with respect to the reference cell can be described as follows. On the interval [ - 1 , 1 ] with vertices v 0 = - 1 and v 1 = 1 , a cubic polynomial p ⁢ ( x ^ ) is uniquely determined by the values of the Hermite node functionals

All degrees of freedom in one dimension are at the end of the interval and the node functionals are dual to the well known basis functions

(2.3)

Higher-dimensional versions of this element are obtained by simple tensor products, see for instance [5]. In this case, the degrees of freedom are all located in the vertices 𝐯 i ∈ 𝒱 , and 𝒱 is the set of vertices of the reference cell [ - 1 , 1 ] d . In three dimensions for instance, we obtain in each vertex the node functionals

(2.4)

We refer to the number of ones in the upper index as the order of the node functional that coincides with the order of the derivative involved. The associated basis functions for instance for v 0 = ( - 1 , - 1 , - 1 ) T are

This may suffice as an example. The general case can easily be constructed using binary vertex indices.

In arbitrary dimension, we use a binary multi-index κ = ( κ 1 , … , κ d ) and obtain the representation

2.3 Basic Properties of Hermitian-Like Mapping

Consider mapping functions Φ K : [ - 1 , 1 ] d → K , which are vector-valued Hermitian-like tensor product elements themselves. That is, each component of the vector Φ K ⁢ ( 𝐱 ^ ) is of the form in equation (2.5).

This choice has the immediate consequence, that the values of Φ K and its Jacobian on any face of the cell are determined by degrees of freedom in vertices adjacent to this face only. Thus, continuity of the Jacobian can be achieved by choosing them consistently on the cells sharing a vertex. Consider for instance a horizontal edge of the reference square, bounded by vertices say 𝐯 0 = ( - 1 , - 1 ) T and 𝐯 1 = ( 1 , - 1 ) T . Then the position of a point 𝐱 = Φ ⁢ ( x ^ , - 1 ) on the image of this edge is a vector-valued, cubic polynomial, determined by the four degrees of freedom Φ ⁢ ( ± 1 , - 1 ) and ∂ x ^ ⁡ Φ ⁢ ( ± 1 , - 1 ) . The tangential derivative ∂ x ^ ⁡ Φ ⁢ ( x ^ , - 1 ) is uniquely determined by the same degrees of freedom. The normal derivative ∂ y ^ ⁡ Φ ⁢ ( x ^ , - 1 ) is a cubic polynomial determined by the degrees of freedom ∂ y ^ ⁡ Φ ⁢ ( ± 1 , - 1 ) and ∂ x ^ ⁢ y ^ ⁡ Φ ⁢ ( ± 1 , - 1 ) . Accordingly, ∇ ⁡ Φ on this edge has the following form:

with the vector-valued polynomials

It remains to choose the actual node values for this transformation. Obviously, Φ ⁢ ( 𝐯 i ) must be the position of the corresponding vertex of the actual cell. The degrees of freedom involving derivatives on the other hand are not predetermined by a standard mesh geometry, with exception of the domain boundary. As a consequence of (2.2), first order degrees of freedom ∂ j ⁡ Φ ⁢ ( 𝐯 i ) determine up to a factor the directions of the edges attached to this vertex. At this point, we observe that the degrees of freedom of higher order and their basis functions only distort the geometry interior to each mesh cell without changing their boundaries, therefore, we already constrain mappings Φ K by

Thus, the polynomial 𝐪 in the gradient of Φ reduces to

We also could have used the second derivatives to simplify 𝐪 , for instance, making it linear. This is impractical though, since these degrees of freedom influence two edges at a time. Thus, the derivative on one edge would depend on an opposing vertex.

3.1 Shape Regular, Conforming Meshes

Most easily, this can be achieved by choosing the vector connecting the two adjacent vertices. In the situation of Figure 1 on the left, the direction vectors for the two edges leaving the center vertex are chosen parallel to the lines connecting the two vertices on these lines. This procedure is not local to each cell anymore, but requires averaging over the cells adjacent to a vertex.

The effect of such a mapping on a ring consisting of three trapezoidal cells can be seen in Figure 2. The shape regularity of the cell, indicated by distortions of the subdivision, is not much worse than for the corresponding bilinear mapping.

As can be seen in Figure 1, this method produces directions in a vertex which may not be orthogonal. While this is not a problem in general, it produces meshes not suitable for higher regularity divergence-conforming elements, see [2, 9]. Due to the anisotropic polynomial spaces used there, the normal derivative with respect to one edge in a corner point must be the tangent to the other edge. This requires that the edges meeting at one vertex are mutually orthogonal.

This can be achieved for instance by the Gram–Schmidt algorithm. Here, we propose the method illustrated in Figure 1 on the right. First, we compute the lines dividing the angles between the two original directions in half. The directions chosen are obtained by rotating those by π 2 . This way, both directions are close to the ones computed in the first step.

After determining the direction, with or without orthogonalization, the norm of the derivative is chosen as the average length of the two edges it represents, as in Figure 2 on the left. Even simpler, it can be chosen as half the distance of the two points used to compute the direction, as on the right. As the figure shows, simplicity comes at the cost of more strongly differing gradients.

3.2 Hanging Nodes

When dealing with locally refined meshes with nonconforming edges/faces, we follow the accepted approach that the degrees of freedom on the refined side must be constrained such that the mappings on both side coincide. This principle now applies to positions as well as derivatives.

Note that the vertex in the center of the patch in Figure 3 is artificial, hence the name hanging vertex. Thus, its position can be chosen to match the requirements of the mapping. While located at the average of its two neighbors’ positions for standard mappings, its location is now determined by the cubic polynomial mapping on the coarse cell, see Figure 3 left and center. This approach has the advantage of being dimension independent, by taking the values of the ( d - 1 ) -cubic polynomials on the face of the coarse cell at the bisecting points. It can be easily extended to more levels of irregularity or to nonuniform subdivision.

The degrees of freedom of the mapping on one of the small cells in the hanging vertex are completely determined by the two vertices left and right of it. This includes not only the position and first derivatives, but also the mixed second derivative, which might be nonzero here.

We point out here, that the parameter ξ on the large cell traverses the edge only once, while it does so twice on the two small cells. Therefore, the derivative degrees of freedom on the refined cells must be divided by two to maintain consistency with the coarse cell. Since adaptive algorithms typically generate a mesh hierarchy by refinement, this can be obtained by taking refinement level information into account. For details, refer to the data structures discussed in Section 4.

In order to be C 1 -conforming, the direction of the edge going upward from the hanging vertex in Figure 3 is determined by the normal derivative of the mapping on the coarse cell below, which in turn is uniquely determined by its values at the vertices. In particular, this edge may or may not be perpendicular to the other one. If we aim at vertices with orthogonal edges, we cannot avoid adding a shape function on the coarse side. Since this implies a fourth-order polynomial along the edge, the same augmentation is required on the refined cells. Such an augmentation can be obtained by a hierarchical approach: first, the bicubic mapping on the large cell is computed according to the previous section. In a second step, degrees of freedom for additional shape functions are determined such that the vectors in ∇ ^ ⁢ Φ are orthogonal in the centers of the edges (and the faces in 3D).

In two dimensions, we augment the mapping by one vector-valued function for each edge, which allows adjusting the normal derivative. First, add to the one-dimensional basis in equation (2.3) the fourth-order polynomial

which has vanishing degrees of freedom 𝒩 i j for i , j ∈ { 0 , 1 } and the value 1 at ξ = 1 2 . Then, for each edge, add the mapping basis function, which is the product of ψ E in tangential direction with the basis function for the derivative degree of freedom in orthogonal direction. For example, we add the basis function

on the left edge E 1 = { - 1 } × [ - 1 , 1 ] which satisfies

For the derivatives at x = - 1 , y = 0 , this choice yields

A straightforward approach is to choose 𝜶 E , 1 such that ∂ x ^ ⁡ Φ is the projection of 𝜶 x onto a vector orthogonal to 𝜶 y . By using the Gram–Schmidt procedure, this implies

The construction extends by tensor product structure to three-dimensional cells, where we have to consider face and edge augmentations of the form

one on each face, two on each edge. The augmentation on the edges is chosen such that for each face the tangential vector not parallel to the edge is orthogonal to the edge. Then these vectors form a two-dimensional system which can be orthogonalized according to Figure 1, obtaining an orthogonal system of three vectors, which form the directions of the edges of the refined cells. This concept generalizes again to higher dimensions due to its tensor product structure.

After this – note that we have modified the tangentials of the faces – the coefficient for the face augmentation follows (3.1), but orthogonalizing with respect to both tangent vectors.

3.3 Boundary Layers

When discretizing problems exhibiting boundary layers, typically an adjusted mesh with anisotropic mesh cells is used. The procedure described in the previous section averages the gradients of all cells adjacent to a vertex. Thus, it can be expected that the approach results in self-intersecting cells if the mesh is not locally quasi-uniform, for instance, if a cell has a neighbor with high aspect ratio. This effect can indeed be observed in Figure 4, where the meshes have been elongated to show the behavior of the mapping on the small cell more clearly. The original mesh on the left consists of two cells of equal height, forming a wedge to highlight the effect discussed. In the subsequent meshes, the lower cell is stretched towards the bottom more and more. While stretching by a factor 2 hardly affects the quality of the mapping in the smaller cell, the mapping gets close to singular at a factor 8, where the self-intersection could still be avoided by changing the derivatives at the top edge. At a stretching of 16, the mapping is clearly self-intersecting.

An obvious fix would be a scaling of the derivative of the mapping to accommodate for the small cell. But, such a change would scale the vertical derivative in the large cell, thus introducing a variation of the Jacobi determinant J depending on the small cell, which would in turn yield bad constants in local interpolation estimates.

The solution is the choice of a mapping, which is not continuously differentiable, but still yields such finite element functions by transformation. Indeed, we do a similar rescaling of derivatives as in the case of hanging vertices, but this time for the normal derivative instead of the tangential one.

This is achieved by defining on each of the two cells K 1 and K 2 sharing a face F a “length orthogonal to F” called h i ; F . Several definitions of such a quantity can be devised. Here, we suggest the following: for each vertex 𝐯 of F, compute the hyperplane obtained from 𝐯 and all other vertices of F sharing an edge with 𝐯 . Then compute the distance of the (unique) vertex of K i sharing an edge with 𝐯 and not in F to this hyperplane. Finally, take the average over all vertices 𝐯 of F. Figure 5 shows cells with the same parameters as in Figure 4, but with rescaled derivatives on the small cell. Clearly, the quality of the mapping does not deteriorate anymore when the large cell is stretched.

The implementation of such a mapping can be achieved by an additional vector-valued scaling factor stored for each cell and multiplied with the degree of freedom for the derivative. However, note that the scaling factor in one direction must be shared by all cells joined by an edge in that direction. Therefore, this scheme is suitable for boundary layers, where we have this situation, but not for anisotropies that may change side along a line.

4.1 The Algorithm

Let a mesh geometry be given by vertex coordinates as well as the tangent planes in boundary vertices.

1. Assign the coordinates of each vertex to the zero order degree of freedom.

2. For each edge between a vertex and a neighboring vertex, record the connecting vector 𝐞 and a scaling factor δ E taking into account anisotropic boundary cells or locally refined cells (see Section 4.3). Ignore hanging vertices in this process. For instance, this can be achieved cell-wise and the first vertex 𝒗 0 on a cell in two dimensions gets the direction vectors

   𝐞 0 0 = 𝒗 1 - 𝒗 0 ,

   𝐞 0 1 = 𝒗 2 - 𝒗 0 ,

   assuming a lexicographic ordering of the vertices of the cell. Note that this choice yields the standard d-linear C 0 -mapping with discontinuous derivatives.

3. For each vertex and each pair of edges e 1 and e 2 with respective scaling factors δ 1 and δ 2 reaching the vertex from opposite sides, compute the corresponding first order degree of freedom as the weighted mean

   𝜶 v * = 1 δ 1 ⁢ 𝐞 1 + 1 δ 2 ⁢ 𝐞 2 1 δ 1 + 1 δ 2 .

   For a shape regular mesh without hanging vertices, this is simply the average. If the vertex is on the boundary, project the obtained vector to the tangent plane.

4. Optionally: orthogonalize the derivative degrees of freedom in each vertex. If the mesh has hanging vertices, compute the additional degrees of freedom according to equation (3.1).

5. Compute the positions and derivatives for hanging vertices according to Section 3.2.

4.2 Examples

The steps of the algorithm as performed by the deal.II finite element library are illustrated in Figure 6 (a)–(c). Subfigure (d) there illustrates, that it is impossible to satisfy both orthogonality constraints and hanging node constraints using cubic tensor product polynomials. There are not enough degrees of freedoms on the coarse cell left to orthogonalize the gradients at hanging nodes. This problem can be solved by enriching the ansatz space by ansatz functions that have zero value and unit derivative both in x- and y-direction at all mid points of faces.

Figure 6

The algorithm as performed by the deal.II library. Subfigures (a) to (c) illustrate the mesh obtained after steps (2), (3), and (5) without orthogonalization, respectively. Subfigure (d) shows the attempt of orthogonalization without enrichment, which is bound to fail.

(a)

Step (1) and (2): The bilinear mapping.

(b)

Step (3): Averaging derivatives.

(c)

Step (5): Computing hanging vertices.

(d)

Orthogonalization with bicubic polynomials.

While the presentation above often used the two-dimensional case as an example, the approach does not depend on the dimension. An example can be seen in Figure 7.

Figure 7

The result of the algorithm in three dimensions computed by the deal.II library.

(a)

Trilinear continuous geometry.

(b)

Tricubic C 1 -geometry.

The last step of this algorithm could be deferred to the point, where a finite element computation is executed on the mesh and coordinates are needed on the smaller cells. In this case, the necessary values are computed on the fly and no data is attached to the hanging vertices. If the data is precomputed as described above, and as implemented in deal.II, the higher order degrees of freedom in these vertices may be nonzero due to interpolation and must be stored as well.

After the degrees of freedom of the mapping in each vertex have been computed using this algorithm, they can be used locally to determine the geometry of each cell. On a plain, shape regular and conforming mesh, this can be done in a straightforward way. If on the other hand the mappings are modified for hanging vertices or anisotropic cells, derivative degrees of freedom must be rescaled.

4.3 The Scaling Factors δ E

Rescaling on shape-regular mesh hierarchies with hanging vertices can be achieved without storing additional data by choosing for each edge δ E = 2 - ℓ , where ℓ is the refinement level of the cells adjacent to E. This level is unique, since hanging vertices are ignored in the computation. Furthermore, this length is characteristic to the cells themselves, such that the “local” coefficients in the expansion (2.5) are obtained from the “global” coefficients computed by the algorithm by multiplying by 2 - ℓ ⁢ o , where o is the number of derivatives related to the node functional.

The same principle can be applied to anisotropic cells in boundary layers according to Section 3.3, but here the anisotropic rescaling must be stored for each cell. Both methods can be combined.