Projection-Based Finite Elements for Nonlinear Function Spaces

# Projection-Based Finite Elements for Nonlinear Function Spaces

Philipp Grohs Philipp Grohs
Universität Wien
Fakultät für Mathematik
Oskar Morgenstern Platz 1
1090 Wien
Austria
Hanne Hardering Hanne Hardering
Technische Universität Dresden
Institut für Numerische Mathematik
Zellescher Weg 12–14
01069 Dresden
Germany
Oliver Sander Oliver Sander
Technische Universität Dresden
Institut für Numerische Mathematik
Zellescher Weg 12–14
01069 Dresden
Germany
and  Markus Sprecher Markus Sprecher
ETH Zürich
Seminar für Angewandte Mathematik
Rämistraße 101
8092 Zürich
Switzerland
###### Abstract.

We introduce a novel type of approximation spaces for functions with values in a nonlinear manifold. The discrete functions are constructed by piecewise polynomial interpolation in a Euclidean embedding space, and then projecting pointwise onto the manifold. We show optimal interpolation error bounds with respect to Lebesgue and Sobolev norms. Additionally, we show similar bounds for the test functions, i.e., variations of discrete functions. Combining these results with a nonlinear Céa lemma, we prove optimal and discretization error bounds for harmonic maps from a planar domain into a smooth manifold. All these error bounds are also verified numerically.

\pdfsuppresswarningpagegroup

=1

AMS classification: 65N30, 65D05

Keywords: geometric finite elements, projection, nonlinear manifold, interpolation errors, discretization errors, harmonic maps

We investigate the discrete approximation of functions from a Euclidean domain to a closed embedded submanifold of , . Such functions are involved in a variety of partial differential equations (PDEs), from fields like liquid crystal physics [4] and micromagnetics [13]. In these applications, the manifold is , the unit sphere in . In Cosserat-type material models [36, 30, 29] the manifold is , where is the special orthogonal group. Further examples are the investigation of harmonic maps into manifolds [6], signal processing of manifold-valued signals [32], and the denoising of manifold-valued images [5].

We are interested in functions of Sobolev smoothness. By this we mean functions from spaces

 Wk,p(Ω,M)\vbox:={v∈Wk,p(Ω,Rn):v(x)∈Ma.e.},

where we denote by the standard Sobolev space for and . Throughout the paper, and will denote the corresponding Sobolev semi norm and full norm of -valued functions, respectively.

Spaces of approximating functions will be constructed by pointwise projection. Given a finite element grid of , and a set of values at Lagrange points on , we construct nonlinear finite element functions by first interpolating in by piecewise polynomials in , and then projecting pointwise onto . This results in a finite-dimensional set of functions which, as it turns out, is a subspace of for arbitrary . While the approach presented here is based on Lagrangian interpolation in , other linear FE space can be used in principle (see [37] for an example).

The idea to generalize finite elements spaces by a pointwise projection operator has already appeared several times [37, 16, 35]. For functions taking values in the special orthogonal group , Gawlik and Leok have studied -norms of interpolation errors [16]. We will extend these results to general closed submanifolds of , and to interpolation errors in Sobolev norms.

To this end, let be the standard nodal interpolation operator for -valued Lagrangian finite elements, and set the interpolation operator with a pointwise projection. For smooth manifolds, approximation qualities can be inferred from the linear ones of , as we can switch back and forth between discrete functions into and into via the definition of and the identity

 QRn∘QM=QRn.

An alternative proof that uses the Lipschitz continuity of the closest-point projection has been given in [37].

A priori, test functions for manifold-valued settings are defined as variations of particular manifold-valued functions. We show that test functions for functions defined by polynomial interpolation and projection can also be constructed directly, using Euclidean interpolation followed by a projection. We show the same Sobolev interpolation error bounds for these discrete test functions as for the finite element functions themselves.

We then discuss finite element discretizations of PDEs with values in . Prototypically, we focus on harmonic maps from a domain to , which we regard as minimizers of the Dirichlet energy in a suitable Sobolev space. The corresponding discrete solution is defined as a local minimizer of the same energy in , which is well-defined because is suitably conforming.

To estimate we combine a simple nonlinear Céa lemma with the interpolation results for . To show optimal bounds we use the abstract theory of [23], showing that the four criteria stated there are fulfilled by projection-based finite elements. We will also provide inverse estimates. In classical finite element theory they are used in many proofs, e.g., in Nitsche’s method of weighted norms for uniform convergence estimates [11]. In this work we will use them to justify a priori bounds on discrete minimizers of the harmonic map energy. Both interpolation and discretization error bounds are verified numerically in the two final chapters.

There is one alternative construction for conforming finite element spaces for manifold-valued problems, known as geodesic finite elements [33, 34, 19, 20]. To evaluate the relative merits of the two methods we briefly revisit their theoretical relationship, and we repeat all numerical tests using geodesic finite elements. We observe that while geodesic finite elements yield lower errors, projection-based finite elements can be much faster.

## 1. Projection-based finite element spaces

Let be discretized by a finite union of affine-equivalent, regular and quasi-uniform polyhedra , such that the closures intersect in common faces. On we consider scalar-valued Lagrangian finite element spaces with the nodal basis and associated Lagrange points .

We will define the space of projection-based finite elements as the image of an interpolation operator. First we consider the canonical interpolation operator for continuous functions with values in into the space of -valued Lagrangian finite elements.

###### Definition 1.

The interpolation operator corresponding to a set of basis functions and nodes is defined by

 QRnv\vbox:\vbox:=∑i∈Iv(ξi)ϕifor all v∈C(Ω,Rn).

For a manifold embedded in and a function the values of will in general not be on away from the . To get -valued functions we compose pointwise with the closest-point projection

 P:Rn→M,P(q)\vbox:\vbox:=argminr∈M∥r−q∥Rn,

where denotes the Euclidean distance. While the closest-point projection is usually not well defined for all , if is regular enough it is well defined in a neighborhood of [1].

This pointwise projection induces a superposition operator by

 (Pv)(x)\vbox:=P(v(x))

for all and . We then define -valued interpolation by composition of and .

###### Definition 2.

Set

 C(Ω,M;ρ)\vbox:={v∈C(Ω,M):diam(v(Th))<ρ ∀Th∈G},

where the denotes the geodesic diameter of a subset . Provided that is small enough, define the interpolation operator

 QM:C(Ω,M;ρ)→C(Ω,M)byQM\vbox:=P∘QRn.

The space of projection-based finite elements is defined as the range of this interpolation operator.

###### Definition 3.

Let , an embedded submanifold, and the closest-point projection. For a given set of basis functions we define

 (1) Vh(Ω,M)\vbox:={vh:Ω→Ms.t.∃v∈C(Ω,M) and vh=QMv}.

As the operator only uses the values at the Lagrange nodes , we have the equivalent definition

 Vh(Ω,M)\vbox:={vh:Ω→M,∃(ci)i∈I⊂M s.t. vh=P(∑i∈Iciϕi)}.

It has to be noted that while there exist nodal values , for any function , for given values , there exists an interpolating function only if the values are close enough depending on such that .

### 1.1. Conformity

The question of conformity of projection-based elements, i.e., whether holds, can be reduced to the continuity of the superposition operator on and of the operator .

Denote by the differential of the closest-point projection at applied to . Let for given coefficients , and . By the chain rule we have

 ∂∂xjvh(x)=P′(∑i∈Iciϕi(x))[∑i∈Ici∂ϕi∂xj(x)]

for every such that is differentiable at and is differentiable at for all .

If we assume that is a smooth embedded submanifold, there exists a tubular neighborhood such that the closest-point projection is smooth [26, Prop. 6.1.8]. In particular, the pointwise norm of can be estimated in terms of the radius of curvature using explicit calculations in terms of the local parametrization of the manifold [1]. Thus, the -conformity of follows directly from the chain rule and smoothness of the Lagrange basis ,

where denotes the operator norm of the differential .

### 1.2. Relationship to geodesic finite elements

Projection-based finite elements are closely related to the geodesic finite elements proposed in [33, 34, 19] and analyzed in [20, 21]. Geodesic finite elements are constructed by replacing polynomial interpolation of values

 x↦∑i∈Iciϕi(x)=argminq∈Rn∑i∈Iϕi(x)∥ci−q∥2Rn

by the weighted Riemannian center of mass

 (2) x↦argminq∈M∑i∈Iϕi(x)d(ci,q)2,

where is the geodesic distance on . Unlike the construction by pointwise projection, (2) is completely intrinsic, and does not rely on an embedding space. Well-posedness of this definition under suitable conditions on the is shown in [34, 21].

As observed independently by [37] and [16], we recover the projection-based interpolation if we replace the geodesic distance in (2) by the Euclidean distance of the embedding space

 argminq∈M∑i∈Iϕi(x)∥ci−q∥2Rn =argminq∈M(∥q∥2Rn−2⟨q,∑i∈Iϕi(x)ci⟩) =argminq∈M∥∥∥q−∑i∈Iϕi(x)ci∥∥∥2Rn =P(∑i∈Iϕi(x)ci).

This does not mean that projection-based finite elements are equal to geodesic finite elements for embedded manifolds. In general, even if the metric on is obtained by an isometric embedding into Euclidean space, the distance is not the Euclidean distance in the surrounding space. Instead, projection-based interpolation can be interpreted as geodesic finite elements for a general metric space with a non-intrinsic metric. As far as we know, no general existence theory and error estimates exist for this abstract setting.

### 1.3. Preservation of isometries

If an isometry commutes with the projection-based interpolation operator then the finite element space defined in (1) is equivariant under this isometry. In mechanics, this leads to the desirable property that discretizations of objective problems are again objective. Unfortunately, for projection-based finite elements this commutativity only holds under special circumstances.

###### Definition 4.

An isometry (w.r.t. the geodesic distance) is called extendable if there exists an isometry with for all .

Examples for extendable isometries are orthogonal transformations for the sphere and multiplication with special orthogonal matrices for .

###### Theorem 5.

Let be a Riemannian submanifold, the closest-point projection, an extendable isometry and a partition of unity. Then commutes with .

###### Proof.

Let be an extension of to . As an isometry maps closest distances to closest distances, commutes with . By the Mazur–Ulam theorem [27] there exists a linear map with for all , so that obviously commutes with . ∎

One can tell from this proof that only very few isometries are extendable. Indeed, in order to be extendable, needs to be the restriction of a rigid body motion of . In contrast to this rather strong restriction, the geodesic interpolation rule (2) is equivariant under any isometry of by construction.

### 1.4. Discrete test functions and vector field interpolation

The test function space for a function consists of vector fields along that correspond to intrinsic variations within the class of functions considered. For , we call the space of test functions . If we consider as an embedded submanifold, then can be canonically identified with a subset of .

We construct discrete test functions in the same manner, i.e., is a discrete test function for if there exists a variation such that for all , , and [21, 35]. Writing this definition using the coefficients that constitute , the set of all discrete test functions over the discrete function can be defined as

 Wh(Ω,u−1hTM)\vbox:\vbox:={vh∈L2(Ω,Rn):∃(ci)i∈I⊂C((−ϵ,ϵ),M) s.t.vh=ddt∣∣∣t=0P(∑i∈Ici(t)ϕi)}.

Similar to the discrete functions themselves, discrete test functions can be constructed by polynomial interpolation followed by pointwise projection, as by chain rule we have for any and

 vh(x)

and , the differential of the closest-point projection , is again a projection.

###### Proposition 6.

Let be the closest-point projection onto a closed embedded -submanifold . For any , the differential is the orthogonal projection onto the tangent space , with the canonical interpretation of as a subspace of .

###### Proof.

Let . We need to show that for all and

 ⟨P′(y)(ξ)−ξ,ω⟩=0

holds. To see this, we consider the curve defined by , and a vector field along with . As is defined by minimization, the first variation yields at any

 ⟨P(y+tξ)−(y+tξ),w(t)⟩=0.

Differentiating this with respect to yields

 0 =ddt∣∣∣t=0⟨P(y+tξ)−(y+tξ),w(t)⟩ =⟨P′(y)(ξ)−ξ,ω⟩+⟨P(y)−y,w′(0)⟩ =⟨P′(y)(ξ)−ξ,ω⟩.\qed

Thus, the computation of the value at of a test functions along a discrete function corresponds to first interpolating given tangent vectors in , and then projecting the resulting piecewise polynomial function pointwise orthogonally to . Alternatively, by linearity we can first project the orthogonally to , and then interpolate the result in the vector space .

In particular, we can define for an interpolation operator by

 Qu−1hTMv=Pu−1hTM∘QRn(v).

Given some test function along a continuous function , i.e. , we can first interpolate and then , as the interpolation of depends only on the values at the Lagrange nodes , where and agree.

## 2. Interpolation error estimates

In this chapter we will estimate the interpolation errors of and in terms of the mesh width . We also estimate the error of the test vector field interpolation operator .

### 2.1. Properties of Euclidean interpolation

Proving interpolation error bounds for uses several standard results for interpolation in Euclidean spaces. We repeat some of them here for convenience.

Define the usual grid dependent Sobolev norms

 ∥v∥Wl,p(G)\vbox:\vbox:=(∑T⊂G∥v∥pWl,p(T))1p

for functions such that for all . For the rest of this paper, this norm is meant whenever we speak of the -norm of a discrete function, unless explicitly stated otherwise. As we assume shape regularity of the mesh, one can use the Sobolev embedding theorem and elementwise scaling to the reference element to prove that if , is continuous with respect to the grid-dependent -norm, i.e., there exists such that we have

 (3) ∥QRnv∥Wl,p(G)≤C∥v∥Wl,p(Ω) for all v∈Wl,p(Ω).

Note that by the Sobolev embedding theorem, is well-defined for all with . Under these assumptions, we have the following approximation error estimate for [11, 10].

###### Theorem 7.

Let be a bounded Lipschitz domain, a shape-regular, affine-equivalent mesh on , , and Lagrangian nodal basis functions for polynomial order . Then on each element for any with and we have

 |v−QRnv|Wl,p(T,Rn)≤Chmin(m,r+1)−lT|v|Wmin(m,r+1),p(T,Rn)∀v∈Wm,p(T,Rn),

with the constant independent of and .

We will also need the following inverse inequalities.

###### Theorem 8.

Consider a shape-regular, affine-equivalent, quasi-uniform mesh and two pairs and with and such that the space of polynomials up to degree on is a subspace of for each mesh element . Then for all discrete functions of polynomial order

 |QRnv|Wm,q(T)≤Ch−(m−k)−smax{0,1p−1q}|QRnv|Wk,p(T),

where the constant depends on the quasi-uniformity and regularity parameters of the mesh, but not on .

This result here is not as general as it could be. Inverse inequalities with weaker requirements of the mesh appear, e.g., in [12] (no quasi-uniformity), and [17] (no shape regularity). We expect that these generalizations can help to extend the following results on -valued interpolation as well.

### 2.2. M-valued interpolation

We now turn to error bounds for the -valued interpolation operator . Given , we estimate the error , , by observing that , and using the triangle inequality

 ∥QMv−v∥Wl,p ≤∥QMv−QRn(QMv)∥Wl,p+∥QRnv−v∥Wl,p,

(again in the grid-dependent norm). Denoting , both terms on the right can be bounded using Theorem 7, and we obtain

 ∥QMv−v∥Wl,p ≤Ch^m−l(|QMv|W^m,p+|v|W^m,p).

It remains to show estimates of in terms of Sobolev norms of . Unlike in the Euclidean case the Sobolev semi-norm is not by itself sufficient to bound , because lower-order derivatives appear by the chain rule. The proper quantity is the homogeneous norm , known, e.g., from [9]. It replaces the unwieldy smoothness descriptor used in corresponding results for geodesic finite elements [20].

###### Proposition 9.

Let , , and such that the closest-point projection is in in some -neighborhood of . Let be a discrete function from a Lipschitz domain into defined by interpolation of values . Suppose the are contained in a geodesic ball of radius , where is small enough such that . Then

 (4)

where is a constant that depends on the -norm of .

###### Proof.

Let be a multi-index with . By the chain rule, the derivative can be written almost everywhere as a sum of terms of the form

 P(k)(QRnv(x))[D→a1QRnv(x),...,D→akQRnv(x)],

where , and . An expansion of around yields

 (5) ∥P′(QRnv)∥L∞≤1+Lip(P′)d(QRnv,M)≤1+CLip(P′)ρ,

where , and denotes the Lipschitz constant of the map . For

 (6) ∥P(k)(QRnv)∥L∞≤L(P)d(QRnv,M)≤CL(P)ρ,

where the constant depends on the -norm of . Further, we have

 ∥∥P(k)(QRnv(x))[D→a1QRnv(x),...,D→akQRnv(x)]∥∥Lp≤∥∥P(k)(QRnv)∥∥L∞k∏i=1∥∥D→aiQRnv∥∥L^mp|→ai|.

For , this yields

 ∥P′(QRnv(x))[D→aQRnv(x)]∥Lp ≤(1+CLip(P′)ρ)∥QRnv∥W^m,p.

For , we have by the Gagliardo–Nirenberg–Sobolev and Young’s inequalities,

 ∥D→aiQRnv∥L^mp|→ai| ≤C|QRnv||→ai|−1^m−1W^m,p|QRnv|1−|→ai|−1^m−1W1,^mp+C|QRnv|W1,^mp

Combining all of this yields

 ∥D→aQMv∥Lp ≤∑∥∥P(k)(QRnv(x))[D→a1QRnv(x),...,D→akQRnv(x)]∥∥Lp ≤(1+CLip(P′)ρ)∥D→aQRnv∥Lp +CL(P)ρ(|QRnv|W^m,p+max{|QRnv|^mW1,^mp,|QRnv|W1,^mp}).\qed

If , then the highest-order derivatives vanish. In that case, (4) reduces to

 |QMv|W^m,p≤CρL(P)max{|QRnv|^mW1,^mp,|QRnv|W1,^mp}

for the mesh-dependent norm.

We can now state the main theorem.

###### Theorem 10.

Consider the same setting as in Theorem 7. Let be an embedded submanifold, such that the closest-point projection is in , where . Then there exists , depending on , and such that for all and

 (7) |v−QMv|Wl,p≤Ch^m−l[|v|W^m,p+L(P)hα∥v∥W^m,p(∥v∥W^m,p+|v|^mW1,^mp)],

with as in Proposition 9, and the constant depending on the constant in that proposition, the one in Theorem 7, as well as , , and .

###### Proof.

For , we have, using (5) and Theorem 7,

 ∥QMv−v∥Lp ≤∫10∥P′(v+t(QRnv−v))(QRnv−v)∥Lpdt ≤(1+CLip(P′)ρ)∥QRnv−v∥Lp ≤C(1+Lip(P′)ρ)h^m|v|W^m,p.

For , we use , Theorem 7, and Proposition 9 to estimate

 ∥QMv−v∥Wl,p ≤∥QMv−QRn(QMv)∥Wl,p+∥QRnv−v∥Wl,p ≤Ch^m−l(|v|W^m,p+|QMv|W^m,p) ≤Ch^m−l(|v|W^m,p+|QRnv|W^m,p +CρL(P)(|QRnv|W^m,p+max{|QRnv|^mW1,^mp,|QRnv|W1,^mp})).

If denote the Lagrangian interpolation nodes in an element , we have by the Sobolev embedding theorem for some

 ρ≤maxThmaxi,jd(v(ξi,Th),v(ξj,Th))≤C∥v∥C0,α|ξi,Th−ξj,Th|α≤C∥v∥W^m,phα.

By Theorem 7, we can estimate all arising semi-norms of by corresponding semi-norms of . Further, we have by the Sobolev embedding theorem . This yields the assertion. ∎

Note that the constants of our estimates are all independent of . The only dependence on is the factor . However, since appears in the error bound (7) only multiplied with , , it becomes irrelevant for . The bounds are therefore optimal in terms of the mesh width. Extra terms compared to the linear result can be controlled by the closeness parameter of the interpolation nodes, and thus for continuous functions by the mesh width parameter .

We have seen that, due to the chain rule, estimates on obtained from the ones on are always with respect to the homogeneous Sobolev seminorms of the type . As the term does not scale correctly, we cannot expect general inverse estimates in the style of Theorem 8 for . An exception is the special case .

###### Theorem 11.

Let the assumptions of Theorem 8 be fulfilled with and . Then for all projected finite element functions we have

 |QMv|W1,q(T)≤Ch−smax{0,1p−1q}|QMv|W1,p(T).
###### Proof.

The estimate follows directly from Proposition 9 and Theorem 8. We need the condition in order to apply Theorem 7 to estimate

 \qed

### 2.3. Tm-valued interpolation

In Section 1.4, we have defined interpolation of a vector field along a discrete function by . This definition is very similar to that of -valued interpolation, with the difference that the pointwise projection is even linear. This linearity makes proving optimal interpolation error bounds for vector fields along given discrete functions much easier than proving the error bounds for the discrete functions themselves.

###### Theorem 12.

Let be a bounded Lipschitz domain, , such that is bounded independently of , , and . Let the assumptions of Theorem 7 be satisfied. Assume further that is in on . Then there exist constants and such that

 |Qu−1hTMv−v|Wl,p≤Chmin(m,r+1)−l|v|Wmin(m,r+1),p[1+Chα(1+∥uh∥Wm,p+∥uh∥mW1,mp)].
###### Proof.

As , we can estimate for , using (5)

 ∥Qu−1hTMv−v∥Lp =(∫Ω∣∣P′(uh(x))(QRnv(x)−v(x))∣∣pdx)1