Accessible Published by De Gruyter November 7, 2020

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

Daniel Arndt ORCID logo and Guido Kanschat ORCID logo


Finite elements of higher continuity, say conforming in H 2 instead of H 1 , require a mapping from reference cells to mesh cells which is continuously differentiable across cell interfaces. In this article, we propose an algorithm to obtain such mappings given a topologically regular mesh in the standard format of vertex coordinates and a description of the boundary. A variant of the algorithm with orthogonal edges in each vertex is proposed. We introduce necessary modifications in the case of adaptive mesh refinement with nonconforming edges. Furthermore, we discuss efficient storage of the necessary data.

MSC 2010: 65N50

1 Introduction

The Bogner–Fox–Schmit element [4] is usually referred to as an element on rectangles. While it is particularly simple there, already Petera and Pittman [8] presented an isoparametric version in 1994. Their mapping to the reference cell requires a global, smooth mapping from a domain with Cartesian cells to a deformed domain. Furthermore, they compute the derivative degrees of freedom by a global smoothing process. Here, we present a simple procedure generating such a mapping based on local information only involving cells adjacent to the same vertex. It extends to nonconforming meshes due to adaptive refinement.

It seems that smooth approximations and mappings have not attracted much attention in the finite element community until recently, most likely inspired by the advent of isogeometric methods. We refer in particular to Sangalli, Takacs, and their coworkers, discussing smooth piecewise spline approximation on irregular triangular and quadrilateral meshes in [6, 7]. Since their analysis is based on Argyris-type degrees of freedom, they require at least polynomials of degree five. Here, we use degree three, which becomes possible by a restriction to simpler topology.

An obvious purpose of a C 1 -mapping is the H 2 -conforming discretization with an isoparametric Bogner–Fox–Schmit element. This can be achieved in a fairly simple way by ensuring that the 2 d edges leaving an interior vertex of a d-dimensional mesh are aligned pairwise.

When the mesh describes a smooth manifold in a higher-dimensional space, the consistency of normal vector of adjacent cells is always an issue, which has originated fixes like averaging of these vectors. Such an averaging improves the consistency error, but does not produce a conforming method, while a C 1 -mapping does.

More complicated is the case of the element introduced by Austin, Manteuffel, and McCormick [2] on rectangles. We have introduced commuting interpolation operators for this element in [9], which at the end of each edge involve the normal derivative of the normal vector component. On the other hand, due to the choice of the anisotropic polynomial space, the Hermitian/Lagrangian interpolation involved requires the tangential derivative along the orthogonal edge, which only works on rectangles, where the two edges are actually orthogonal. This property can be maintained if we impose additionally on the mapping that all edges leaving a vertex are either orthogonal or aligned. Therefore, below we describe two versions of the mapping, without and with orthogonalization.

The mapping without additional orthogonalization transfers easily to adaptive meshes with hanging nodes. There, position and derivative values are obtained by extension of the standard interpolation known from multilinear mapping. With additional orthogonalization, a d-cubic space for the mapping is not anymore sufficient and an enrichment is necessary.

The purpose of this article is twofold. On the one hand, we describe the algorithmic construction of a piecewise polynomial C 1 -mapping and the conditions that have to be met. On the other hand, we discuss an efficient layout of data for such a mapping. In Section 2, we recall the notion of smooth mappings and then focus on the properties of polynomial mappings on a single cell. In Section 3, we discuss the construction of such mappings based on given vertex coordinates. We close with remarks on the implementation in Section 4.

2 Mapped Polynomials and Degrees of Freedom

We assume that the domain Ω d is subdivided into a mesh 𝕄 of mesh cells K. Our approach is based on tensor products and works in arbitrary dimension. The examples though are restricted to the most relevant dimensions d = 2 , 3 . While we will discuss the exact shape of the cells below, we assume that the positions of the vertices of the mesh and their connectivity are given. Furthermore, let the boundary Ω be smooth and its tangent plane be known in every point. Each mesh cell K will be represented by a smooth bijection Φ K : K ^ K from the reference cell K ^ = [ - 1 , 1 ] d into the domain d . The whole mesh is then described by the function Φ on the Cartesian product of the reference cell, such that for each mesh cell K there is a factor of this product with Φ K ( K ^ ) = K , thus

Φ : K 𝕄 K ^ d ,

where the range of Φ approximates Ω. Two cells K i and K j are called neighbors of each other if the face F i j = K i K j has dimension d - 1 . We call the mesh conforming if each interior face F i j is a whole facet of each of the adjacent cells. The discussion of nonconforming meshes is deferred to Section 3.2. Let P = [ - 2 , 2 ] d and let T map the 2 d -fold Cartesian product K ^ × × K ^ to P such that the reference cells K ^ are mapped into one orthant of P by translation by the vectors ( ± 1 , , ± 1 ) T with all possible combinations of signs and possibly rotation around their centers by multiples of π 2 . A mesh is called regular if for every interior vertex v, there are 2 d mesh cells K v i and a continuous mapping Φ v such that


The set P together with the mappings Φ v forms an atlas for Ω. Therefore, we define that the mapping Φ is C 1 if each of the mappings Φ v is in C 1 ( P ) .

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

u ( 𝐱 ) = u ^ ( 𝐱 ^ ) , u ( 𝐱 ) = ^ Φ K ( 𝐱 ^ ) - T ^ u ^ ( 𝐱 ^ ) ,

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:


𝐭 K ( 𝐱 ) = ^ Φ ( 𝐱 ^ ) 𝐭 ^ K ^ ( 𝐱 ^ ) ,
𝐧 K ( 𝐱 ) = ^ Φ ( 𝐱 ^ ) - T 𝐧 ^ K ^ ( 𝐱 ^ ) .

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

( ^ Φ 1 ( Φ 1 - 1 𝐱 ) ) - T ^ u ^ 1 ( 𝐱 ^ ) = ( ^ Φ 2 ( Φ 2 - 1 𝐱 ) ) - T ^ u ^ 2 ( 𝐱 ^ ) .

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

𝒩 0 0 ( p ) = p ( - 1 ) , 𝒩 0 1 ( p ) = p ( - 1 ) ,
𝒩 1 0 ( p ) = p ( 1 ) , 𝒩 1 1 ( p ) = p ( 1 ) .

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


4 ψ 0 0 ( ξ ) = ξ 3 - 3 ξ + 2 , 4 ψ 0 1 ( ξ ) = ξ 3 - ξ 2 - ξ + 1 ,
4 ψ 1 0 ( ξ ) = - ξ 3 + 3 ξ + 2 , 4 ψ 1 1 ( ξ ) = ξ 3 + ξ 2 - ξ + 1 .

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


𝒩 i 000 ( p ) = p ( 𝐯 i ) ,
𝒩 i 100 ( p ) = x ^ p ( 𝐯 i ) , 𝒩 i 010 ( p ) = y ^ p ( 𝐯 i ) , 𝒩 i 001 ( p ) = z ^ p ( 𝐯 i ) ,
𝒩 i 110 ( p ) = x ^ y ^ p ( 𝐯 i ) , 𝒩 i 101 ( p ) = x ^ z ^ p ( 𝐯 i ) , 𝒩 i 011 ( p ) = y ^ z ^ p ( 𝐯 i ) ,
𝒩 i 111 ( p ) = x ^ y ^ z ^ p ( 𝐯 i ) .

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

ψ 0 κ λ μ ( x ^ , y ^ , z ^ ) = ψ 0 κ ( x ^ ) ψ 0 λ ( y ^ ) ψ 0 μ ( z ^ ) .

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.5) p ( 𝐱 ^ ) = 𝐯 i 𝒱 κ α i κ j ψ i κ j ( x ^ j ) .

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:

Φ ( ξ , - 1 ) = ( x ^ 𝐩 ( ξ ) , 𝐪 ( ξ ) ) ,

with the vector-valued polynomials

𝐩 ( x ^ ) = Φ ( 𝐯 0 ) ψ 0 0 ( x ^ ) + x ^ Φ ( 𝐯 0 ) ψ 0 1 ( x ^ ) + Φ ( 𝐯 1 ) ψ 1 0 ( x ^ ) + x ^ Φ ( 𝐯 1 ) ψ 1 1 ( x ^ ) , 𝐪 ( x ^ ) = y ^ Φ ( 𝐯 0 ) ψ 0 0 ( x ^ ) + x ^ y ^ Φ ( 𝐯 0 ) ψ 0 1 ( x ^ ) + y ^ Φ ( 𝐯 1 ) ψ 1 0 ( x ^ ) + x ^ y ^ Φ ( 𝐯 1 ) ψ 1 1 ( x ^ ) .

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

𝒩 i α ( Φ K ) = 0 for all  | α | > 1 .

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

𝐪 ( x ^ ) = y ^ Φ ( 𝐯 0 ) ψ 0 0 ( x ^ ) + y ^ Φ ( 𝐯 1 ) ψ 1 0 ( x ^ ) .

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 Globally Differentiable Mapping of Mesh Cells

The main challenge for the construction of C 1 -conforming elements on arbitrary quadrilaterals and hexahedra consists in finding a C 1 -conforming mapping Φ K from the reference cell K ^ to the grid cell K, namely continuity of Φ and Φ across interfaces. In order to be able to use polynomial orders as low as d-cubic, we confine ourselves to the case of meshes with regular vertices in the sense of equation (2.1) with possibly adaptive refinement. Irregular vertices require mapping with higher polynomial degrees, see for instance [7].

We follow the isoparametric approach, choosing a d-cubic mapping defined by Hermitian-like interpolation. Clearly, the values of Φ in the vertices 𝐱 ^ i of K ^ must be the vertices 𝐱 i of T. Furthermore, continuity of the tangential derivatives of Φ along interfaces is implied by continuity of Φ itself. It remains to ensure continuity of the normal derivative.

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.

Figure 1

Determining the direction of edges leaving a vertex. Simple averaging (left) and orthogonalized (right). Auxiliary lines dotted, final directions in red. Dashed lines on the right are final directions on the left.

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.

Figure 2

A triangular ring with bicubic mapping and different scaling of derivatives. Average distance from vertex to its neighbors (left) and distance of the two neighboring vertices (right).

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.

Figure 3

A one-irregular edge with bilinear C 0 -mapping (left) and bicubic C 1 -mapping (center) and augmented mapping with orthogonality at center point (right). The red dot indicates the location of the hanging vertex with linear interpolation on the edge.

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

ψ E ( ξ ) = 16 ξ 2 ( 1 - ξ ) 2 ,

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

ψ E , 1 ( x ^ , y ^ ) = ψ 0 1 ( x ^ ) ψ E ( y ^ ) ,

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

ψ E , 1 ( - 1 , y ^ ) = 0 , x ^ ψ E , 1 ( - 1 , ± 1 ) = 0 ,
y ^ ψ E , 1 ( - 1 , y ^ ) = 0 , x ^ ψ E , 1 ( - 1 , 0 ) = 1 .

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

x ^ Φ = 𝜶 0 10 x ^ ψ 0 10 + 𝜶 0 11 x ^ ψ 0 11 + 𝜶 2 10 x ^ ψ 2 10 + 𝜶 2 11 x ^ ψ 2 11 + 𝜶 E , 1 = : 𝜶 x + 𝜶 E , 1 ,
y ^ Φ = 𝜶 0 01 y ^ ψ 0 01 + 𝜶 0 11 y ^ ψ 0 11 + 𝜶 2 01 y ^ ψ 2 01 + 𝜶 2 11 y ^ ψ 2 11 = : 𝜶 y .

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

(3.1) x ^ Φ = 𝜶 x - 𝜶 x 𝜶 y 𝜶 y 2 𝜶 y 𝜶 E , 1 = - 𝜶 x 𝜶 y 𝜶 y 2 𝜶 y .

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

ψ E / F , ( x ^ , y ^ , z ^ ) = ψ 0 1 ( x ^ ) ψ E ( y ^ ) ψ E ( z ^ ) ,

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.

Figure 4

An increasingly anisotropic wedge shaped mesh. The algorithm of the previous section is applied to the two center vertices. The bottom cells have heights 1, 2, 8, and 16, respectively (from left to right). Horizontal positions are constant.

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.

Figure 5

An increasingly anisotropic wedge shaped mesh with corrected mapping. The bottom cells have heights 2, 4, 8, and 16, respectively (from left to right). Horizontal positions are constant

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 Implementation

The C 1 -conforming meshes described above are implemented in the open-source finite element library deal.II [1, 3]. In particular, the mapping is implemented in terms of the MappingFEField class that allows to describe the geometry by a finite element vector.

The mesh data structure in d dimensions requires d + 1 vector-valued degrees of freedom in each non-hanging vertex, one of order zero for the position of the vertex and d of order one for the derivatives. While Hermite interpolation uses even higher order degrees of freedom, these are set to zero for mapping purposes. Using the local numbering of degrees of freedom described in (2.4), these are assigned by the following algorithm. Note that the algorithm uses additional intermediate data in each vertex, which can be deleted after completion in order to free space for numerical computation.

4.1 The Algorithm

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

  1. (1)

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

  2. (1)

    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. (1)

    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. (1)

    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. (1)

    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.

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


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

                     Step (3): Averaging derivatives.


Step (3): Averaging derivatives.

                     Step (5): Computing hanging vertices.


Step (5): Computing hanging vertices.

                     Orthogonalization with bicubic polynomials.


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.

                     Trilinear continuous geometry.


Trilinear continuous geometry.



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.


We have constructed a global C 1 -mapping based on tensor products of cubic Hermite interpolation for quadrilateral and hexahedral meshes. Such a mapping is possible under the condition that all vertices are topologically regular. The construction easily extends to hanging vertices due to adaptive mesh refinement and to anisotropic cells in boundary layers. The algorithm can be implemented in two sweeps over the mesh using only local information in each of them.

Funding source: Deutsche Forschungsgemeinschaft

Award Identifier / Grant number: 1648 SPPEXA

Funding statement: The authors gratefully acknowledge support of DFG through priority program 1648 SPPEXA. Guido Kanschat gratefully acknowledges support by the special program “Numerical Analysis of Complex PDE Models in the Sciences” of the Erwin Schrödinger International Institute for Mathematics and Physics (ESI) at University of Vienna.


[1] G. Alzetta, D. Arndt, W. Bangerth, V. Boddu, B. Brands, D. Davydov, R. Gassmoeller, T. Heister, L. Heltai, K. Kormann, M. Kronbichler, M. Maier, J.-P. Pelteret, B. Turcksin and D. Wells, The deal.II library, Version 9.0, J. Numer. Math. 26 (2018), 10.1515/jnma-2018-0054. Search in Google Scholar

[2] T. M. Austin, T. A. Manteuffel and S. McCormick, A robust multilevel approach for minimizing 𝐇 ( d i v ) -dominated functionals in an 𝐇 1 -conforming finite element space, Numer. Linear Algebra Appl. 11 (2004), no. 2–3, 115–140. Search in Google Scholar

[3] W. Bangerth, R. Hartmann and G. Kanschat, deal.II – A general purpose object oriented finite element library, ACM Trans. Math. Softw. 33 (2007), 10.1145/1268776.1268779. Search in Google Scholar

[4] F. K. Bogner, R. L. Fox and L. A. Schmit, The generation of interelement compatible stiffness and mass matrices by the use of interpolation formulae, Proceedings of Conference on Matrix Methods in Structural Mechanics, Air Force Institute of Technology, Wright Patterson (1965), 397–444. Search in Google Scholar

[5] F. Bonizzoni and G. Kanschat, H1-conforming finite element cochain complexes and commuting quasi-interpolation operators on Cartesian meshes, preprint (2020), Search in Google Scholar

[6] A. Collin, G. Sangalli and T. Takacs, Analysis-suitable G 1 multi-patch parametrizations for C 1 isogeometric spaces, Comput. Aided Geom. Design 47 (2016), 93–113. Search in Google Scholar

[7] M. Kapl, G. Sangalli and T. Takacs, The Argyris isogeometric space on unstructured multi-patch planar domains, preprint (2017), Search in Google Scholar

[8] J. Petera and J. F. T. Pittman, Isoparametric Hermite elements, Internat. J. Numer. Methods Engrg. 37 (1994), 3489–3519. Search in Google Scholar

[9] N. Sharma and G. Kanschat, A contraction property of an adaptive divergence-conforming discontinuous Galerkin method for the Stokes problem, J. Numer. Math. 26 (2018), no. 4, 209–232. Search in Google Scholar

Received: 2020-10-07
Revised: 2020-10-20
Accepted: 2020-10-27
Published Online: 2020-11-07
Published in Print: 2021-01-01

© 2021 Walter de Gruyter GmbH, Berlin/Boston