MULTIGRID METHODS FOR DISCRETE ELLIPTIC PROBLEMS ON TRIANGULAR SURFACES

MULTIGRID METHODS FOR DISCRETE ELLIPTIC PROBLEMS

ON TRIANGULAR SURF ACES

RALF KORNHUBER AND HARRY YSERENTANT

Dedicated to Wolfgang Hackbusch on the occasion of his60th birthday.

Abstract.We construct and analyze multigrid methods for discretized self-

adjoint elliptic problems on triangular surfaces in R3.The methods involve

the same weights for restriction and prolongation as in the case of planar trian-

gulations and therefore are easy to implement.We prove logarithmic bounds

of the convergence rates with constants solely depending on the ellipticity,the

smoothers and on the regularity of the triangles forming the triangular surface.

Our theoretical results are illustrated by numerical computations.

1.Introduction

Geometric di?erential equations play a crucial role in many applications ranging from material science to image processing or numerical relativity[9,12].Numerical discretizations of such problems typically lead to large algebraic systems which become ill–conditioned with decreasing mesh size.For example,the approximation of an evolving surface driven by mean curvature requires the solution of a second order elliptic problem on a triangular surface in each time step[5].

A straightforward approach to construct multigrid methods on triangular sur-faces is to simply use the same weights for restriction and prolongation as for pla-nar triangulations.Such algorithms have been implemented in the software package MC and applied successfully to the numerical solution of the Einstein equations[10]. Special coarsening strategies for unstructured triangular surface meshes have been considered in[1].In spite of the simplicity of the straightforward approach,nu-merical experiments indicate multigrid convergence speed.However,there seems to be no theoretical justi?cation yet.Biorthogonal wavelet bases on manifolds pro-vide an essential step towards mesh independent preconditioners[4].However,the construction and thus the resulting algorithms involve piecewise smooth parame-terizations of the underlying manifold which might cause problems,if the manifold itself has been computed numerically.

In this paper,we provide a convergence analysis for a class of multigrid methods for discretized self-adjoint elliptic problems on a triangular surface M j.As these multigrid methods involve the same weights for restriction and prolongation as for planar triangulations,our results can regarded as a theoretical justi?cation of the above-mentioned straightforward approach.The main di?culty in the analysis is

2KORNHUBER AND YSERENTANT

resulting from the fact that an underlying sequence M k,k=1,...,j,of coarser tri-angular surfaces does no longer generate a sequence of nested?nite element spaces. As a consequence,the existing convergence theory for subspace correction meth-ods[13,14]cannot be applied directly.The main idea of this paper is to generate suitable decompositions of functions on M j by decomposing associated functions on a re?ned reference con?guration M′j.Assuming that M′j is resulting from suc-cessive planar re?nement of coarse triangles forming a reference con?guration M′0, nested?nite element spaces on M′j can be obtained,e.g.,by standard nodal in-

terpolation.Existing estimates for this hierarchy provide the desired estimates for an associated hierarchy on M j.In such a way,we are able to derive logarithmic bounds for the convergence rates.The constants solely depend on the ellipticity,the smoothers and on the regularity of the triangles forming the triangular surface M j. Moreover,we obtain exactly the same weights as in the planar case for restriction and prolongation.

The paper is organized as follows.After stating the problem,we introduce the concept of logically nested triangular surfaces,clarifying the connection of M j and M′j.Section4is devoted to a hierarchical decomposition by generalized interpo-lation.In Section5we prove logarithmic upper bounds for the convergence rates of a corresponding hierarchical basis multigrid method and discuss some possible variants and improvements.A concluding numerical experiment illustrates our theoretical?ndings.

2.Discrete elliptic problems on triangular surfaces

Let M j?R3denote a surface consisting of planar triangles t∈T j.To?x the ideas,M j can be regarded as an approximation of some continuous surface M with j denoting the number of re?nement steps.The space S j,

S j={v∈C(M j)|v|t is linear?t∈T j},

of linear?nite elements on M j is equipped with the usual Sobolev norms

v 20,M

j

= M j v(x)2dx,|v|21,M j= ?M j v 20,M j, v 21,M j= v 20,M j+|v|21,M j,

which are de?ned piecewise here,that is,triangle by triangle.The tangential gradient?M j v:t→R3of a function v:R3→R is the tangential part

?M j v=?v?(ν·?v)ν

of the gradient of v,whereνdenotes the normal of the actual triangle t∈T j.It depends only on the values of v on t and therefore can be evaluated in an obvious manner for functions de?ned only on t.Correspondingly,the parts of the given norms depend only on the relative position of the vertices of the single triangles to each other and can be evaluated as in the planar case.

We consider the discrete variational problem

(2.1)u j∈S?j:a(u j,v)=?(v)?v∈S?j.

Here,S?j?S j is a subspace of S j,?is a linear functional on S?j,and a(·,·)is a symmetric,S?j–elliptic bilinear form.More precisely,

(2.2)α v 21,M

j ≤a(v,v)≤β v 21,M

j

?v∈S?j

MULTIGRID METHODS FOR DISCRETE ELLIPTIC PROBLEMS ON TRIANGULAR SURFACES3

holds with positive constantsα,βso that the energy norm · ,de?ned by

v 2=a(v,v)

is equivalent to · 21,M

j .Usually,both the bilinear form a(·,·)and the right hand

side?depend on the triangular surface M j and therefore on j.For example,the bilinear form

(2.3)a(v,w)= M j?M j v(x)·?M j w(x)dx

is generated by the Laplace–Beltrami operator.It satis?es(2.2),if the boundary of M j is non-empty and homogeneous Dirichlet boundary conditions are prescribed. Moreover,the constantsα,βare independent of j,if M j converges to a su?ciently smooth surface M in a suitable way.On these conditions,it is well-known that u j converges to the solution of the continuous analogue of(2.1)with the same convergence rates as in the planar case[7].

For ease of presentation,we assume from now on that S?j=S j.

3.Logically nested triangular surfaces

Let M′0be a conceptionally coarse,triangular surface consisting of non-degene-rate triangles t′∈T′0and let

φ:M′0→M

denote a parametrization of a continuous surface M over M0.We assume that M′0is conforming in the sense that the intersection of two triangles t,t′∈T′0is either a common edge,a common vertex or empty.Self-intersections of M′0are not excluded and M′0may have a boundary or not.An example is shown in the left picture of Figure3.1.The only condition that we impose on the mappingφis later given implicitly,in De?nition3.2below.

Let T′0,T′1,...,T′j be a sequence of nested triangulations of M′0as resulting from standard red/green re?nement of T′0(see,e.g.,[2,6],or[11,p.66]).A triangle t with the vertices p i∈R3is denoted by t=t(p1,p2,p3).For each k=0,...,j,we identify each t′∈T′k with an associated triangle t?R3according to

(3.1)T′k?t′=t′(p′1,p′2,p′3)?t=t(p1,p2,p3),p i=φ(p′i).

Note that di?erent triangles t associated with di?erent triangles t′∈T′k may overlap or occupy the same place in space.Moreover,it is not excluded at this point that certain associated triangles t degenerate.

De?nition3.1.A sequence of triangular surfaces M0,M1,...,M j formed by triangulations T0,T1,...,T j is called logically nested,if there exists an associated sequence of nested triangulations T′0,T′1,...,T′j such that

(3.2)T k={t?R3|?t′∈T′k:t′?t},k=0,...,j.

is valid.

As an example,let M0be some triangular approximation of M such that all vertices of all triangles t∈T0are located on M.Then a(uniformly re?ned)logically nested sequence M0,M1,...,M j is obtained by inductively bisecting the edges of each t∈T k,shifting the midpoints to M in a suitable way and then connecting the shifted midpoints.In this case,we can simply chose M′0=M0(see Figure3.1).

4KORNHUBER AND YSERENTANT

Figure3.1.Reference con?guration M′0with re?nement M′1and

associated logically nested triangular surface M1

As the coarse approximations M0,M1,...,M j?1do not play any role in the rest of this paper,we will simply say from now on that M j is logically nested. The nested triangulations T′k give rise to nested?nite element spaces

S′0?S′1?···?S′j.

Each?nite element space S′k is spanned by the nodal basis functionsλ(k)p′associated with the nodes p′∈N′k.Due to self-intersection,two or more di?erent nodes might occupy the same point in space.We introduce the mapping F:S j→S′j by (3.3)S j?v→v′=F v∈S′j,(F v)(p)=v(φ(p))?p∈N′j,

which transforms each function v∈S j into an associated function v′∈S′j with the same nodal values.If certain triangles t∈T j degenerate,then this mapping is not one-to-one.This motivates the following de?nition.

De?nition3.2.A logically nested surface M j is called regular,if the norm equiv-alence

(3.4)γ v l,t≤ F v l,t′≤Γ v l,t?t∈T j?v∈S j

holds for l=0,1with positive constantsγ,Γ.The ratioΓ/γquanti?es the regu-larity of M j.

The estimates(3.4)describe how much the size and shape of the triangles t∈T j and t′∈T′j associated to each other can di?er.The constantsγandΓare intention-ally independent of the given re?nement level j and?x indirectly the requirements on the parametrizationφof the surface M.For example,these conditions are sat-is?ed,if M is a piecewise smooth surface such that the restrictionsφ|T,T∈T′0, are su?ciently smooth,regular functions.

We are ready to state the main result of this section.

Proposition3.1.Assume that M j is logically nested and regular.Then the map-ping F:S j→S′j de?ned by(3.3)is linear and bijective and has the properties

(3.5)γ v l,M

j ≤ F v l,M′

≤Γ v l,M

j

?v∈S j,l=0,1.

MULTIGRID METHODS FOR DISCRETE ELLIPTIC PROBLEMS ON TRIANGULAR SURFACES 5

In particular,F maps subsets of linear independent functions of S j onto subsets of linear independent functions of S ′j and vice versa.As a consequence,φ:N ′k →

N k =φ(N ′k )mapping the nodes N ′k of S ′k onto the nodes N k of the subspace F ?1(S ′k )?S j is one-to-one for k =0,...,j .

Exploiting Proposition 3.1,the given problem (2.1)can be easily reformulated as an equivalent problem on the re?ned reference con?guration M ′j .Then,multigrid

algorithms could be derived and analyzed by working on M ′j ,where the nested

sequence S ′0?···?S ′j is available.However,such a strategy would involve the

transformation F and thus the parametrization φexplicitly .In order to avoid that,we prefer to work directly on M j .

4.Hierarchical decomposition

We assume that M j is logically nested and regular as explained in the preced-ing section.As a starting point for multilevel methods we now provide a stable decomposition of S j into nested subspaces V 0?V 1?···?V j ?S j .A correspond-ing decomposition of the reference space S ′j is easily obtained,e.g.,by standard nodal interpolation.The main idea is to transform this decomposition together with well-known stability estimates from S ′j to S j .

Let I ′k :S ′j →S ′k denote the standard nodal interpolation.Utilizing Proposi-

tion 3.1,we introduce the generalized interpolation operators

(4.1)I k :S j →S j ,I k =F ?1I ′k F,k =0,...,j.

Observe that I j v =v and I k v is de?ned on M j and not on M k for k

V 0=I 0S j ,V k =(I k ?I k ?1)S j ,k =1,...,j.Proposition 4.1.The hierarchical decomposition (4.2)is stable in the sense that

c 0 v 21,M j ≤ I 0v 21,M j +j k =1

4k (I k ?I k ?1)v 20,M j ≤c 1(j +1)2 v 21,M j

holds for all v ∈S j with positive constants c 0,c 1depending only on the regularity of M j expressed in terms of the ratio Γ/γfrom De?nition 3.2.

Proof.Let v ∈S j .Denoting v ′=F v ,v ′0=I ′0v ′,and v ′k =(I ′k ?I ′k ?1)v ′for

k =1,...,j ,we get from the triangle inequality and the Cauchy-Schwarz inequality v ′ 0,M ′j ≤ j k =0 v ′k 0,M ′j ≤( j k =04?k )1/2( j k =04k v ′k 20,M ′j )1/2≤23( v ′0 21,M ′j + j k =14k v ′k 20,M ′j

)1/2.A local version of the well-known strengthened Cauchy-Schwarz inequality (see,e.g.,the proof of [14,Lemma 2.7]or [15,Theorem 3.4])and an inverse inequality yield

|v ′|21,t ′≤c ′( v ′0 21,t ′+

j k =1

4k v ′k 20,t ′)?t ′∈T ′0.The constant c ′depends only on the shape and size of the triangles in T ′0.Now the left inequality follows from Proposition 3.1and from the ellipticity (2.2).

6KORNHUBER AND YSERENTANT

It is well–known that for all v′∈S′j and all t′∈T′0the estimates

I′0v′ 21,t′≤C′0(j+1)2 v′ 21,t′

and

4k v′?I′k v′ 20,t′≤C′1(j?k+1) v′ 21,t′,k=1,...,j,

hold with constants C′0,C′1>0depending only on the shape and size of t′.We refer,e.g.,to[15,Theorem3.1and Theorem3.2]or[16]and the references cited therein.Now the assertion follows from the triangle inequality,Proposition3.1and the ellipticity(2.2).

5.Multilevel methods

Utilizing the general framework of subspace correction schemes[13,16],the hierarchical splitting

S j=V0⊕V1⊕···⊕V j

provided in(4.2)gives rise to an hierarchical basis multigrid method for the discrete variational problem(2.1)with the bilinear form a(·,·)and the right hand side?: Starting with the initial residual r j=??a(uνj,·)of the given iterate uνj,we compute corrections v k from the approximate defect problems

v k∈V k:b k(v k,v)=r k(v)?v∈V k,k=j,j?1, 0

For k

r k?1=(r k?a(v k,·))|V

k?1

,k=j,j?1, (1)

Finally,we collect the corrections from all levels to obtain the new iterate

(5.1)uν+1

j =uνj+

j

k=0

v k.

The additive version of(5.1)provides an hierarchical basis preconditioner associated with the bilinear form

(5.2)b(v,w)=

j

k=0

b k(v k,w k),v=

j

k=0

v k,w=

j

k=0

w k,v k,w k∈V k.

This preconditioner is evaluated by simply skipping the update of the residual r k in the corresponding multigrid algorithm.

The symmetric,positive de?nite bilinear forms b k(·,·)on V k are usually called smoothers.We chose

b0(·,·)=a(·,·)

which means that the defect problems on the coarsest space V0are solved exactly. We assume that the remaining smoothers b k(·,·),k=1,...,j,ful?ll the conditions (5.3)a(v,v)≤ωb k(v,v)?v∈V k,ω<2,

and

(5.4)c24k v 20,M

j ≤b k(v,v)≤c34k v 20,M

j

?v∈V k.

Standard smoothers like Jacobi-or symmetric Gau?-Seidel iterations have these properties.

MULTIGRID METHODS FOR DISCRETE ELLIPTIC PROBLEMS ON TRIANGULAR SURFACES7 Theorem5.1.The hierarchical basis multigrid method(5.1)satis?es the error estimate

2≤(1?c(j+1)?3) u j?uνj 2,ν=0,1,....

(5.5) u j?uν+1

j

The hierarchical basis preconditioner b(·,·)satis?es

(5.6)C1(j+1)?2b(v,v)≤a(v,v)≤C2b(v,v)?v∈S j.

The constants c,C1,C2depend only on the ellipticity of a(·,·),on the smoothers b k(·,·),and on the regularity of M j expressed in terms of the ratioΓ/γfrom De?-nition3.2.

Proof.The proof follows from general convergence results for subspace correction methods(cf.Xu[13]or also Yserentant[16]).More precisely,utilizing(5.4), Proposition4.1,and the ellipticity(2.2),we get the left inequality in(5.6)with C1=α?1c1c3or,equivalently,the condition(5.2)in[16]with K1≤C1(j+1)2.

The ellipticity(2.2),the left inequality in Proposition4.1,and5.3immediately yield the right inequality in(5.6)with C2=c?12c?10βor,equivalently the condition (5.7)in[16]with K2≤C1.Now we can apply Theorem5.4in[16]to prove (5.5).

According to Theorem5.1,the convergence rateρj of(5.1)and the condition numberκj of(5.2)behave like

(5.7)ρj=1?O (j+1)?3 ,κj=O (j+1)2 ,

respectively.The order(j+1)2of the condition number is optimal for hierarchi-cal basis preconditioners.The slightly suboptimal bound forρj is caused by the loss of certain orthogonality properties and thus of strengthened Cauchy-Schwarz inequalities for the subspaces V k as de?ned on triangular surfaces.

The asymptotic bounds in(5.2)can be further improved by selecting larger subspaces V k associated with the new nodes p∈N k\N k?1and their neighbors.In the special case of uniform re?nement this means

V k=F?1S′k,k=1,...,j,

providing a generalization of classical multigrid methods(cf.Hackbusch[8]).The results stated in Theorem5.1directly extend to these methods.In addition,the upper bound for K1appearing in the proof now can be improved on additional assumptions on M j,exploiting that the underlying reference con?guration M′0 is not unique.For example,M′0depicted in the left picture of Figure3.1could be replaced by the planar reference con?guration M0′as shown in Figure5.1. Our assumption is that M j allows for such a planar reference con?guration M′0. Then we immediately get an analogue of Proposition4.1with an upper bound depending only on(j+1)and therefore the condition(5.2)in[16]with K1= O(j+1).The proof relies on generalized L2-projection Q k=F?1Q′k F instead of interpolation I https://www.wendangku.net/doc/0d17356613.html,ing estimates involving the K-functional as proposed by Bornemann and Yserentant[3]on the re?ned reference triangulation M′j,we can even achieve condition(5.2)in[16]with K1independent of j.This leads to

ρj=1?O (j+1)?1 ,κj=O(1).

In many cases,a planar reference con?guration is not available,e.g.,for closed sur-faces.Therefore one might relax this assumption by claiming the existence of an

8KORNHUBER AND YSERENTANT

Figure5.1.A planar reference con?guration M0′

atlas of local planar reference con?gurations.Mesh-independent bounds of conver-gence rates seem to require a di?erent type of analysis referring to the continuous surface M and not to a reference con?guration M′0.

6.Implementation

The implementation of the multigrid methods as derived in the previous section requires a reformulation in terms of vectors,sti?ness matrices,restrictions and prolongations.To this end,we de?ne the generalized nodal basis functions

μ(k)p=F?1λ(k)p′,p=φ(p′)∈N k,k=0,...,j.

Recall thatλ(k)

p′,p′∈N′k,is the standard nodal basis of S′k and that the mapping

φ:N′k→N k is one-to-one as a consequence of Proposition3.1.The generalized nodal basis functionsμ(k)p,p∈N k,are a basis of the spaces

V k=F?1S′k=span{μ(k)p|p∈N k},k=1,...,j,

providing the generalized classical multigrid method.The interpolation operators de?ned in(4.1)have the representation

I k v= p∈N k v(p)μ(k)p,v∈S j,

so that

V0=span{μ(0)p|p∈N0},V k=span{μ(k)p|p∈N k\N k?1},k=1,...,j, are basis representations of the subspaces V k generating the hierarchical basis multi-grid method(5.1).

The recursive formula

μ(k?1)

p

= q∈N kλ(k?1)p′(q′)μ(k)q

reveals that for both multigrid methods exactly the same weightsλ(k?1)

p′(q′)as in the

planar case are used for restriction and prolongation.These weights are available from the logical re?nement structure alone.They do not involve any information about a reference con?guration M′0or a parametrizationφwhich appear only in the analysis.Roughly speaking,multigrid can be carried out on triangular surfaces

in the same way as on planar triangulations.

MULTIGRID METHODS FOR DISCRETE ELLIPTIC PROBLEMS ON TRIANGULAR SURFACES

9 Figure7.1.Final surface M8and corresponding solution u8

Figure7.2.Maximal aspect ratiosσj of M j and convergence

ratesρj over the re?nement levels j

7.Numerical experiments

As a numerical example,we consider the bilinear form(2.3)generated by the

Laplace-Beltrami operator on approximations M j of a jelly?sh-like section of the unit sphere.Coarse approximations for j=0,1are depicted in Figure3.1and

the?nal approximation M8with328193unknowns is shown in the left picture

in Figure7.1.We impose homogeneous Dirichlet conditions on the boundary and chose the right hand side?(v)= M j v(x)dx.The right picture in Figure7.1 illustrates the corresponding solution u j for j=8.

According to our analysis in Section5,the shape regularity of the triangles of M j plays a crucial role for the performance of multilevel methods on M j.The left picture in Figure7.2shows the shape regularity of M j as expressed by the maximal aspect ratioσj of the triangles t∈T j over the re?nement level j.The aspect ratios saturate with increasing re?nement.The right picture in Figure7.2shows the convergence rates of a V(1,0)cycle of the generalized classical multigrid method with symmetric Gau?-Seidel smoother.The convergence rates seem to saturate at

10KORNHUBER AND YSERENTANT

aboutρj≈0.5which is very similar to the performance of this algorithm in the planar case.

References

[1]B Aksoylu,A.Khodakovsky,and P.Schr¨o der.Multilevel solvers for unstruc-

tured surface meshes.SIAM https://www.wendangku.net/doc/0d17356613.html,put.,26:1146–1165,2005.

[2]R.E.Bank,T.F.Dupont,and H.Yserentant.The hierarchical basis multigrid

method.Numer.Math.,52:387–404,1988.

[3]F.A.Bornemann and H.Yserentant.A basic norm equivalence in the theory

of multilevel methods.Numer.Math.,64:455–476,1993.

[4]W.Dahmen and R.Schneider.Wavelets on manifolds.I:Construction and

domain decomposition.SIAM J.Numer.Anal.,31:184–230,1999.

[5]K.Deckelnick,G.Dziuk,and https://www.wendangku.net/doc/0d17356613.html,putation of geometric partial

di?erential equations and mean curvature?ow.Acta Numer.,14:139–232, 2005.

[6]P.Deu?hard,P.Leinen,and H.Yserentant.Concepts of an adaptive hierar-

chical?nite element code.IMPACT Comput.Sci.Engrg.,1:3–35,1989.

[7]G.Dziuk.Finite elements for the Beltrami operator on arbitrary surfaces.In

S.Hildebrandt and R.Leis,editors,Partial di?erential equations and calculus of variations,volume1357of Lecture Notes in Mathematics,pages483–490.

Springer-Verlag,1988.

[8]W.Hackbusch.Multi-Grid Methods and Applications.Springer,Berlin,1985.

[9]S.Hildebrandt and H.Karcher,editors.Geometric analysis and nonlinear

partial di?erential equations.Springer,2003.

[10]M.J.Holst.Adaptive numerical treatment of elliptic systems on manifolds.

https://www.wendangku.net/doc/0d17356613.html,put.Math.,15:139–191,2001.

[11]R.Kornhuber.Adaptive Monotone Multigrid Methods for Nonlinear Varia-

tional Problems.Teubner,Stuttgart,1997.

[12]G.Sapiro.Geometric partial di?erential equations and image analysis.Cam-

bridge University Press,Cambridge,2006.

[13]J.Xu.Iterative methods by space decomposition and subspace correction.

SIAM Review,34(4):581–613,December1992.

[14]H.Yserentant.On the multi-level splitting of?nite element spaces.Numer.

Math.,49:379–412,1986.

[15]H.Yserentant.Two preconditioners based on the multilevel splitting of?nite

element spaces.Numer.Math.,58:163–184,1990.

[16]H.Yserentant.Old and new convergence proofs for multigrid methods.Acta

Numer.,pages285–326,1993.

Prof.Dr.Ralf Kornhuber,Freie Universit¨a t Berlin,Institut f¨u r Mathematik II, Arnimallee2-6,D-14195Berlin,Germany

E-mail address:kornhuber@math.fu-berlin.de

Harry Yserentant,Technische Universit¨a t Berlin,Institut f¨u r Mathematik,Stra?des 17.Juni136,D-10165Berlin,Germany

E-mail address:yserentant@math.tu-berlin.de