Adaptive finite element methods for the Laplace-Beltrami problem

(1)

Technical report, IDE1224, May 2012

Adaptive finite element methods for the Laplace-Beltrami problem

Master’s Thesis in Computational Science and Engineering David Samvin

School of Information Science, Computer and Electrical Engineering Halmstad University

(2)

(3)

Adaptive finite element methods for the Laplace-Beltrami problem

Master’s Thesis in Computational Science and Engineering

School of Information Science, Computer and Electrical Engineering Halmstad University

Box 823, S-301 18 Halmstad, Sweden

May 2012

(4)

(5)

Acknowledgements

I would like to express my gratitude to my tutor Prof. Bertil Nilsson and Prof.

Peter Hansbo for all the help and support.

(6)

(7)

Abstract

This thesis concerns the Laplace-Beltrami problem with FEM on two types of surfaces.

The first part of the thesis is dedicated to explain the background and the general theory. In the second part I will show how to implement the theory with Matlab, and in the third part we consider a closed surface with no boundary conditions and a continuous surface with boundary conditions. The closed surface will be approximated with triangle elements and the continuous surface with quadratic elements.

(8)

(9)

1. Background

1.1

Finite Element Method, FEM

The Finite Element Method, FEM, also known as Finite Element Analysis, FEA, is a numerical method for finding approximate solutions of partial differential equations by dividing the domain into finite elements and computes each element locally then combines the entire system to find a solution. The method was developed initially, and prospered, as a computer-based simulation method for the analysis of aerospace structures, but found its way into both design and analysis of complex structural systems.

There are a lot of different types of elements when using FEM, in both 2D and 3D problems. This thesis concerns a 2D problem embedded in ℝ³and therefore we will be using triangular and quadrilateral elements shown in Fig.1.

Fig.1 - Triangular and Quadrilateral elements

(12)

2

1.2

Laplace-Beltrami operator

Laplace–Beltrami operator is a linear operator and is defined as the divergence of the gradient. The operator can be extended to operate on tensors as the divergence of a derivative along tangent vectors of a manifold.

Divergence measures the magnitude of a vector field’s source or sink at a given point.

div 𝑓 = ∇ ∙ 𝑓 =

^𝜕𝑈_𝜕𝑥

+

^𝜕𝑉_𝜕𝑦

+

^𝜕𝑊_𝜕𝑧

.

^(1.1)

Gradient of a scalar field is a vector field that points in the direction of the greatest rate of increase of the scalar field.

grad 𝑓 = ∇𝑓 =

^𝜕𝑓_𝜕𝑥

𝐢 +

^𝜕𝑓_𝜕𝑦

𝐣 +

^𝜕𝑓_𝜕𝑧

𝐤 .

^(1.2) where i, j and k stands for the unit vector.

(1.1) and (1.2) put together gives the Laplace –Beltrami operator,

∆𝑓 = div grad 𝑓 .

^(1.3)

1.3

Bezier surface

Bezier surface are mathematical splines used in finite element modeling and computer aided design. The surface is defined by a set of control points. In general, the surface passes through the first and last control

(13)

3

point but not through the center control points, rather, it is “stretched”

towards them.

A two-dimensional Bezier surface can be defined as a parametric surface where the position of a point 𝕡 as a function of the parametric coordinates u, v is given by

𝕡(𝑢, 𝑣) = � � 𝐵

_𝑖^𝑛

𝑚

𝑗=0 𝑛

𝑖=0

(𝑢) 𝐵

_𝑗^𝑚

(𝑣) 𝕜

_𝑖,𝑗

A given surface of order (𝑛, 𝑚) is defined by a set of (𝑛 + 1)(𝑚 + 1) control points 𝕜_𝑖,𝑗, where

𝐵

_𝑖^𝑛

(𝑢) = �𝑛𝑖� 𝑢

^𝑖

(1 − 𝑢)

^𝑛−𝑖 is a Bernstein polynomial, and

�𝑛𝑖� = 𝑛!

𝑖! (𝑛 − 𝑖)!

is the binomial coefficient.

(14)

4

2. Introduction

In this paper we consider the problem

∆

_Γ

𝑢 = 𝑓 𝑜𝑛 Γ

^(2.1)

𝑢 = 0 𝑜𝑛 𝜕Γ

Where Γ a two-dimensional projected surface is embedded in ℝ³ and

∆

_Γ

is the Laplace-Beltrami operator on Γ.

For now we consider the case of a closed surface in which case

𝜕Γ = 0

and there are no boundary conditions, but at the end we will also consider a case of an open surface in which case

𝜕Γ ≠ 0

and boundary conditions are present.

The problem can be written on weak form as seeking 𝑢 ∈ 𝐻₀¹ ( Γ ) such that

(∇

_Γ

𝑢, ∇

_Γ

𝑣) = (𝑓, 𝑣)

_Γ

, ∀ 𝑣 ∈ 𝐻

₀¹

Γ

^(2.2) where

∇

_Γdenotes the surface gradient and (. , . )_Γ the 𝐿₂scalar product over Γ.

(15)

5

3. Finite element approximation

3.1

Finite element formulation

We introduce a triangulation 𝒯_ℎ = { 𝑇 } of Γ into a geometrically conforming, quasiuniform, finite element mesh of Γ, resulting in a discrete approximation Γ_ℎ of Γ. Denote by ℎ_𝑇 the diameter of element 𝑇 and by ℎ = max_{𝑇 ∈ 𝒯}_ℎℎ_𝑇 the global mesh size parameter. We shall use continuous, piecewise polynomial, approximations of 𝑢

𝑉

^ℎ

= {𝑣 ∈ 𝐻

¹

( Γ

_ℎ

) ∶ 𝑣 |

_𝑇

∈ 𝑃

^𝑘

( 𝑇 ) for all 𝑇 ∈ 𝒯

_ℎ

}

or

𝑉

₀^ℎ

= {𝑣 ∈ 𝐻

¹

( Γ

_ℎ

) ∶ 𝑣 |

_𝑇

∈ 𝑃

^𝑘

( 𝑇 ) for all 𝑇 ∈ 𝒯

_ℎ

, 𝑣 = 0 on 𝜕Γ

_ℎ

}

and seek 𝑢^ℎ ∈ 𝑉^ℎ such that

(∇

_Γ_ℎ

𝑢

^ℎ

, ∇

_Γ_ℎ

𝑣)

_Γ_ℎ

= (𝑓, 𝑣)

_Γ_ℎ

, ∀ 𝑣 ∈ 𝑉

^ℎ ^(3.1)

In the following, we shall use finite element technology to compute the approximation

_Γ_ℎ

of

_Γ

. We then first assume that

𝑢^ℎ

can be written as

𝑢

^ℎ

= � 𝑢

_𝑖^ℎ

𝑖

𝜑

_𝑖

(𝜉, 𝜂),

(16)

6

where

(𝜉, 𝜂)

are the coordinates of a reference element,

𝜑

_𝑖

(𝜉, 𝜂)

^{are the} reference shape functions, and

𝑢

_𝑖^ℎ are the nodal values of 𝑢^ℎ. We shall next use this parameterization also to approximate Γ.

3.2

Parameterization

Unlike Dziuk [3] and Demlow [1] we consider an isoperimetric parameterization of the surface (the same idea can be used for arbitrary parameterizations). We first assume that we can write, on each 𝑇, 𝑥 = 𝑥(

𝜉, 𝜂

^{) where}

(𝜉, 𝜂)

denotes the coordinates used to define the reference element and 𝑥 = (𝑥, 𝑦, 𝑧) are the physical coordinates.

For any given parameterization, we can extend it outside the surface by defining

𝑥

_surf

(𝜉, 𝜂) = 𝑥(𝜉, 𝜂),

𝑥

_out

(𝜉, 𝜂) = 𝑥(𝜉, 𝜂) + 𝑡 𝑛(𝜉, 𝜂)

where 𝑛 is the normal and 𝑡 is a fixed parameter. In some models, where the surface is an idealized thin structure, it is natural to think of 𝑡 as a thickness.

For the purpose of the current extension, we may simply set 𝑡 = 1. Next, we think of the extention as a function of

(𝜉, 𝜂, 𝜁)

^{so that}

𝑥

_out

(𝜉, 𝜂, 𝜁) = 𝑥

_surf

(𝜉, 𝜂) + 𝜁 𝑛(𝜉, 𝜂), 0 ≤ 𝜁 ≤ 1.

(17)

7

Using 𝑛 = 𝑥_out − 𝑥_surf the extension of the coordinates can alternatively be written

𝑥(𝜉, 𝜂, 𝜁) = (1 − 𝜁)𝑥

_surf

(𝜉, 𝜂) + 𝜁 𝑥

_out

(𝜉, 𝜂),

which represents an affine map given 𝑥_out. We now specialize to parameterizations of the isoperimetric finite element kind, i.e., approximations of the type

𝑥(𝜉, 𝜂) ≈ 𝑥

^ℎ

(𝜉, 𝜂) = � 𝑥

_𝑖

𝑖

𝜑

_𝑖

(𝜉, 𝜂) ∈ [𝑉

^ℎ

]

³

,

𝑥

_out

(𝜉, 𝜂) ≈ 𝑥

_out^ℎ

(𝜉, 𝜂) = � 𝑥

_𝑖^out

𝑖

𝜑

_𝑖

(𝜉, 𝜂) ∈ [𝑉

^ℎ

]

³

,

where 𝑥_𝑖 are the physical location of the nodes on the surface. 𝑥_out^ℎ will never be used but is needed in the following for the construction of the jacobian of the map. We have

𝑥(𝜉, 𝜂, 𝜁) ≈ �(𝑥

_𝑖

𝑖

𝜑

_𝑖

(𝜉, 𝜂)(1 − 𝜁) + 𝑥

_𝑖^out

𝜑

_𝑖

(𝜉, 𝜂) 𝜁)

At 𝜁 = 0 we find, in the usual way, that the physical derivatives of the shape functions on the surface are approximated, at (𝜉, 𝜂) by

⎣⎢

⎢⎢

⎡𝜕

𝜑

_𝑖

𝜕𝜕𝑥

𝜑

_𝑖

𝜕𝜕𝑦

𝜑

_𝑖

𝜕𝑧 ⎦⎥⎥⎥⎥⎥⎤

= 𝕁^−𝟏(

𝜉, 𝜂, 0

⁾

⎣⎢

⎢⎢

⎢⎡ 𝜕

𝜑

_𝑖

𝜕

𝜉

𝜕

𝜑

_𝑖

𝜕

𝜂

𝑛_𝑧^ℎ

𝜑

_𝑖_⎦^⎥

⎥⎥

⎥⎤

(18)

8

Where the Jacobian

𝕁(𝜉, 𝜂, 0) ∶=

⎣ ⎢

⎢ ⎢

⎡

^𝜕

𝑥

𝜕

𝜉

𝜕

𝑦

𝜕

𝜉

𝜕

𝑧

𝜕

𝜉

𝜕

𝑥

𝜕

𝜂

𝜕

𝑦

𝜕

𝜂

𝜕

𝑧

𝜕

𝜂

𝜕

𝑥

𝜕

𝜁

𝜕

𝑦

𝜕

𝜁

𝜕

𝑧

𝜕

𝜁⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ ,

Finally, the three components of the surface gradient are approximated, elementwise, by

∇

_Γ_ℎ

𝜑

_𝑖

= ( 𝕀 − 𝑛

^ℎ

⨂ 𝑛

^ℎ

) ∇𝜑

_𝑖

We also apply the usual FE approximation over the elements for the normals at nodes, 𝕟^ℎ = ∑_𝑖

𝜑

_𝑖 𝕟_𝑖 followed by normalization. Then

𝕁(𝜉, 𝜂, 0) ∶=

⎣ ⎢

⎢ ⎢

⎢ ⎡

^𝜕

𝑥

^ℎ

𝜕

𝜉

𝜕

𝑦

^ℎ

𝜕

𝜉

𝜕

𝑧

^ℎ

𝜕

𝜉

𝜕

𝑥

^ℎ

𝜕

𝜂

𝜕

𝑦

^ℎ

𝜕

𝜂

𝜕

𝑧

^ℎ

𝜕

𝜂 𝑛

_𝑥^ℎ

𝑛

_𝑦^ℎ

𝑛

_𝑧^ℎ

⎦ ⎥ ⎥ ⎥ ⎥ ⎤

.

(19)

9

4. Adaptivity

As adaptivity criteria we introduce the following error estimation,

‖𝑢 − 𝑈 ‖

_𝐿₂

≤ 𝑐

₁

ℎ ‖𝑓‖

_𝐿₂

‖𝑓‖

_𝐿₂

= �� 𝑓

²

𝑑Ω

Ω

�

1/2

= �� 𝑓

²

𝑑Ω

_𝑒

Ω_𝑒

�

1/2

We only consider according to the error, the highest 30 % of the elements in each adaptivity step.

Element refinement is done by adding a new node between the two existing ones for the elements with highest error as shown in Fig.2.

Fig.2

(20)

10

5. Implementation

5.1

Surface

The surface can be generated by using commercial software like Comsol or Solidworks. From the software we can extract the node coordinates, element topology and the normals at nodes.

In our case we will generate the surface using Bezier surface in Matlab and Mathematica and extract the information needed for the nodes, elements and normals.

From the surface we receive the surface coordinates (𝑥, 𝑦, 𝑧) at each node in

𝑛𝑜𝑑𝑒𝑠 = � 𝑥

₁

. ..

𝑥

_𝑛

𝑦

₁

. ..

𝑦

_𝑛

𝑧

₁

. ..

𝑧

_𝑛

�

and the element topology in

𝑒𝑙𝑚𝑛𝑡𝑠 =

⎣ ⎢

⎢ ⎡ 𝑖

₁

. ..

𝑖

_𝑛

𝑗

₁

. ..

𝑗

_𝑛

𝑘

₁

. ..

𝑘

_𝑛

⎦ ⎥ ⎥ ⎤

where each row is an element, and the normals

𝑛𝑜𝑚𝑎𝑙𝑠 = � 𝑛𝑥

₁

. ..

𝑛𝑥

_𝑛

𝑛𝑦

₁

. ..

𝑛𝑦

_𝑛

𝑛𝑧

₁

. ..

𝑛𝑧

_𝑛

�.

(21)

11

5.2

Parametric domain

The elements used in the parametric domain are in two dimensions and of rectangular or triangular shape as shown in Fig.3.

The linear basis functions 𝜑_𝑖(

𝜉, 𝜂

) in the parametric domain is obtained by

𝜑

₁

= 1 − 𝜉 − 𝜂, 𝜑

₂

= 𝜉, 𝜑

₃

= 𝜂

with derivatives

𝜕𝜑

_𝑖

𝜕𝜉 = � 𝜕𝜑

₁

𝜕𝜉

𝜕𝜑

₂

𝜕𝜉

𝜕𝜑

₃

𝜕𝜉 � = [−1 0 1]

𝜕𝜑

_𝑖

𝜕𝜂 = � 𝜕𝜑

₁

𝜕𝜂

𝜕𝜑

₂

𝜕𝜂

𝜕𝜑

₃

𝜕𝜂 � = [−1 1 0]

Due to linearity of the problem a one point Gauss quadrature is enough, thus points and weights

𝑔𝑝 = [

¹₃

,

¹₃

], 𝑔𝑤 = [

¹₂

]

Fig.3

(22)

12

Numerical integration is applied using the area measure

𝑑Ω

_𝑒

= � 𝜕𝕩

𝜕𝜉 × 𝜕𝕩

𝜕𝜂 �

5.3 Surface projecting

Moving the normals from the nodes to the Gauss points can be done by inserting the Gauss points as 𝜉 and 𝜂 then using the standard FE map

𝑛𝐺 = ∑ 𝜑

_𝑖 _𝑖

(𝜉, 𝜂)𝑛

_𝑖

= [𝜑

1

𝜑

₂

𝜑

₃

] �

𝑛

_𝑥_𝑖

𝑛

_𝑦_𝑖

𝑛

_𝑧_𝑖

𝑛

_𝑥_𝑗

𝑛

_𝑦_𝑗

𝑛

_𝑧_𝑗

𝑛

_𝑥_𝑘

𝑛

_𝑦_𝑘

𝑛

_𝑧_𝑘

� (5.1)

and the projection matrix

ℙ = 𝑛𝐺

^T

𝑛𝐺

(5.2)

and the gradient in three dimensions we obtain

ℂ =

⎣ ⎢

⎢ ⎢

⎡

^𝜕𝜑_𝜕𝑥^𝑖

𝜕𝜑_𝑖

𝜕𝑦

𝜕𝜑_𝑖

𝜕𝑧

⎦ ⎥ ⎥ ⎥ ⎤

(5.3)

Projecting the gradient down on the surface using (5.1), (5.2) and (5.3)

𝔹 = (𝕀 − ℙ) ℂ

where 𝕀 is a 3 by 3 identity matrix.

(23)

13

6. Convergence in a smooth case

6.1

Model 1

We consider a sphere with unit diameter and origin at (1/2, 1/2, 1/2). The exact solution

𝑢 = (𝑥 − 1/2)(𝑦 − 1/2)(𝑧 − 1/2)

inserted into (2.1) then corresponds to the right-hand side

𝑓 = (6 (2 𝑥 − 1) (2 𝑦 − 1) (2 𝑧 − 1))

(3 + 4 (𝑥 − 1) 𝑥 + 4 (𝑦 − 1) 𝑦 + 4 (𝑧 − 1) 𝑧)

Using 𝑓 in the FE assembly we should obtain an approximate solution 𝑈 close to the exact solution 𝑢.

6.2

Model 2

We consider an implicitly defined bubble F(𝑥, 𝑦, 𝑧) ≔ (𝑥 − 𝑧²)²+ 𝑦²+ 𝑧²− 1 = 0 with normal 𝕟 = _|∇F|^∇F. The exact solution

𝑢 = 𝑥 𝑦 ( 𝑧 + 1)

²

inserted into (2.1) using Mathematica then corresponds to the right-hand side

𝑓 = 𝑠𝑒𝑒 𝑎𝑡𝑡𝑎𝑐ℎ𝑚𝑒𝑛𝑡. 1

Using 𝑓 in the FE assembly we should obtain an approximate solution 𝑈 close to the exact solution 𝑢.

(24)

14

6.3

Model 3

A smooth open Bezier surface with the same source term as model 2 and boundary conditions 𝑢 = 0.

6.4

Model 4

A smooth open Bezier surface with the same source term as model 1 and boundary conditions 𝑢 = 0.

(25)

15

7. Solution

7.1

Model 1

A sphere with origin at (1/2, 1/2, 1/2) and no boundary conditions, triangular mesh.

(26)

16

(27)

17

By refining the elements and computing 𝑢, we obtain the exact solution as shown below.

(28)

18

(29)

19

The approximate FE solution U

(30)

20

(31)

21

The difference between the exact solution and the approximated solution i

s

𝑒𝑟𝑟𝑜𝑟 = 𝑢 − 𝑈

(32)

22

(33)

23

7.2

Model 2

A bubble with origin at (0, 0, 0) and no boundary conditions, triangular mesh.

(34)

24

(35)

25

By refining the elements and computing 𝑢, we obtain the exact solution as shown below.

(36)

26

(37)

27

The approximate FE solution U

(38)

28

(39)

29

The difference between the exact solution and the approximated solution i

s

𝑒𝑟𝑟𝑜𝑟 = 𝑢 − 𝑈

(40)

30

(41)

31

7.3

Model 3

Open surface with boundary conditions 𝑢 = 0, triangular mesh.

(42)

32

7.4

Model 4

Open surface with boundary conditions 𝑢 = 0, quadratic mesh.

(43)

33

8. Conclusion

A general adaptive finite element method for the Laplace-Beltrami problem has been developed. The method has been implemented to take care of smooth open and closed surfaces. Comparison with analytical results manifests the method in both cases.

(44)

34

9. References

[1] A. Demlow, Higher order finite element methods and pointwise error estimates for elliptic problems on surfaces, SIAM J. Numer.

Anal. 47 (2009) 805-827.

[2] A. Demlow, G. Dziuk, An adaptive finite element method for the Laplace-Beltrami operator on implicitly defined surfaces, SIAM J. Numer. Anal. 45 (2007) 421-442.

[3] G. Dziuk, Finite elements for the Beltrami operator on arbitrary surfaces, in: Partial Differential Equations and Calculus of Variations (Lecture Notes in Mathematics, Vol. 1357),

Hildebrand S, Leis R (eds). Springer: Berlin, 1988, pp. 142-155.

[4] J. Foley, A. van Dam, S. Feiner, J. Hughes, Computer graphics principles and practice (second edition).

(45)

35

10. Attachment 1

Source term for the bubble, model 2.

𝑓 =

(46)

36

11. Attachment 2 - Main

The code is implemented in Matlab 7.1 R2010a.

(47)

37

(48)

38

(49)

39

12. Attachment 3 – Assemble

Adaptive finite element methods for the Laplace-Beltrami problem