An Interface-fitted Mesh Generator and Polytopal Element Methods for - - PowerPoint PPT Presentation

An Interface-fitted Mesh Generator and Polytopal Element Methods for Elliptic Interface Problems

Long Chen

University of California, Irvine chenlong@math.uci.edu Joint work with: Huayi Wei (Xiangtan University), Min Wen(UCI)

Polytopal Element Methods in Mathematics and Engineering Oct 27, 2015, Atlanta, GA

Outline

Elliptic Interface Problems Interface-fitted Mesh Generation 2D Algorithm 3D Algorithm VEM for Elliptic Interface Problems Weak Formulation Virtual Element Method Weak Galerkin Methods Numerical Results

ELLIPTIC INTERFACE PROBLEMS

A Domain with an Interface

Ω + Ω− Γ n

Ω

Figure: A square domain Ω with an interface Γ in it.

Elliptic Interface Problems

Consider −▽ · (β∆u) = f, in Ω\Γ (1) with prescribed jump conditions across the interface Γ: [u]Γ = u+ − u− = q0, [βun]Γ = β+u+

n − β−u− n = q1,

and boundary condition u = g

• n ∂Ω.

Existing Work

Two type of numerical methods:

◮ Numerical methods based on Cartesian meshes. ◮ Numerical methods based on interface-fitted meshes.

Cartesian Mesh Approach

Figure: A Catesian mesh and an interface.

Cartesian Mesh Approach

Pro:

◮ Mesh generation is very simple

Con:

◮ Need to modify the stencil (FDM) or basis function (FEM)

• f the vertex near the interface Γ.

◮ The linear system may be non-symmetric and cause

problems in fast solvers.

◮ Convergence analysis is complicated.

Interface-fitted Mesh

Figure: An interface-fitted mesh in 2D [Wei, Chen, Huang and Zheng,

• SISC. 2014].

Interface-fitted Mesh Approach

Pro:

◮ Can use standard finite element methods to discretize the

interface problem.

◮ Symmetric system can be easily solved by fast solver such

as algebraic multigrid solvers.

◮ Relatively easy error analysis.

Con:

◮ Generate interface-fitted meshes, extremely difficult in 3D!

Our Goal

In the past decade the mesh generation have gotten great

• progress. Here we aim to
• 1. Develop effective and robust interface-fitted (polytopal)

mesh generation algorithm.

• 2. Solve interface problems accurately using conforming FEM

(e.g. Virtual Element method) on polytopal meshes.

• 3. Solve the resulting algebraic system quickly using

(algebraic) multigrid solver.

Outline

Elliptic Interface Problems Interface-fitted Mesh Generation 2D Algorithm 3D Algorithm VEM for Elliptic Interface Problems Weak Formulation Virtual Element Method Weak Galerkin Methods Numerical Results

INTERFACE-FITTED MESH GENERATION

TWO DIMENSIONAL CASE

Delaunay Triangulation

p1 p2 p3 p4 p5 p6 p7

Figure: A Delaunay triangulation of a set of points is a triangulation satisfying the empty circle condition.

Delaunay Triangulation

(a) Non Delaunay. (b) Delaunay.

Figure: Empty circle condition: a circle circumscribing any Delaunay triangle does not contain any other input points in its interior.

2D Interface-fitted Mesh Generation Algorithm

Algorithm 1 2D Interface-fitted Mesh Generation Algorithm INPUT: Grid size: h; Level set function, ϕ(x); Square domain, Ω; OUTOUT: Interface-fitted mesh T ;

1: Find the cut points, the Cartesian mesh points near or on

the interface and some auxiliary points.

2: Construct a Delaunay triangulation on these points. 3: Post processing.

An Example in 2D: Step 1

Figure: Step 1: Find the cut points (red), the mesh points near or on the interface (black) and the auxiliary points (magenta).

An Example in 2D: Step 2

Figure: Step 2: Generate a Delaunay triangulation on these points. In MATLAB, it is simply DT = delaunay(x, y).

An Example in 2D: Step 3

Figure: Step 3: Keep the triangles in the interface elements.

An Example in 2D: Step 4

Figure: An interface-fitted mesh composed by triangles and squares.

Features of 2D Mesh Generator

◮ Simple, efficient and semi-unstructured. ◮ Recovery the interface. ◮ The maximum angle is uniformly bounded (135◦).

Maximum Angle Condition

(a) θmax ≤ 90◦ (b) θmax ≤ 135◦. (c) θmax ≤ 135◦. (d) θmax ≤ 135◦.

Figure: Maximum angles of four types of interface elements. From Semi-Unstructured Grids by C. Pflaum. Computing 2001.

Interface Recovery

Interface Recovery

• We DO NOT code the triangulation case by case. We simply

call DT = delaunay(x, y).

Interface Recovery

• We DO NOT code the triangulation case by case. We simply

call DT = delaunay(x, y). We can PROVE the interface will be preserved in the triangulation generated by delaunay and the maximal angle is minimized.

Lower Convex Hull

(a) Step 1: Lift points to the paraboloid (b) Step 2: Form the lowest convex hull in Rn+1. (c) Step 3: Project the lowest convex hull to Rn.

Figure: Delaunay triangulation is the projection of the lower convex hull of points lifted to the paraboloid f = x2.

Interface Recovery

• The lower convex hull when lift to Rn+1 will always connect the

intersection points except ...

A Degenerate Case

A B C D

Figure: Add an auxiliary point in a rectangle to preserve the interface

Efficiency

Simple, Efficient and Semi-Unstructured.

◮ The greatest advantage: localization. Only call delaunay

for O( √ N) points near the interface. The complexity will be thus O( √ N) in 2-D which can be ignored comparing with the O(N) complexity for assembling the matrix and solving the matrix equation.

◮ No mesh smoothing is needed. The maximal angle is

bounded by 135◦. Good enough for finite element methods.

◮ The mesh is only unstructured near the interface and the

majority of the mesh is structured which leads to nice properties (superconvergence, multigrid etc).

INTERFACE-FITTED MESH GENERATION

THREE DIMENSIONAL CASE

Challenges in Tetrahedra Mesh Generation

Figure: Tetrahedra elements [Shewchuk, 2002].

Most Tetrahedra mesh generation/optimization algorithms cannot guarantee sliver-free for a general 3D domain.

Advertisement: Use ODT mesh smoothing which produces fewer slivers in 3D. L. Chen. Mesh smoothing schemes based on optimal Delaunay triangulations. 2004

3D Interface-fitted Mesh Generation Algorithm

Algorithm 2 3D Interface-fitted Mesh Generation Algorithm

INPUT: Grid size: h; Level set function, ϕ(x); Cube domain, Ω; OUTOUT: Interface-fitted mesh T ;

1: Find the cut points, the Cartesian mesh points near or on the in-

terface and the auxiliary points.

2: Construct a Delaunay triangulation on these points. 3: Post processing: use polytopal meshes.

Polytopal Mesh v.s. Tetrahedron Mesh

Figure: Bad tetrahedron elements will be eliminated and part of their faces will become the boundary of polytopes. As a surface mesh, the triangles are more acceptable.

Features of 3D Mesh Generator

◮ Simple, efficient and semi-unstructured. ◮ Recovery the interface. ◮ The maximum angle of surface mesh is uniformly bounded. ◮ Boundary faces of a polytope is either triangle or square.

Sphere Interface

Figure: Sphere surface with 3, 128 triangles: maximum angle 112.81◦ and minimum angle 11.37◦.

Sphere Interface

Figure: The mesh near the interface. Use 0.338 second to generate a mesh with 70, 169 points).

Heart interface

Figure: Heart surface with 3480 triangles: maximum angle 119.89◦ and minimum angle 4.21◦.

Heart interface

Figure: The mesh near the interface. Use 0.339 second to generate a mesh with 70, 521 points).

Quartics interface

Figure: Quartics surface with 44, 512 triangles: maximum angle 123.9◦ and minimum angle 11.2◦.

Quartics interface

Figure: The mesh near the interface. Use 3.183 second to generate a mesh with 553, 161 points.

Torus interface

Figure: Torus surface with 37, 600 triangles: maximum angle 121.12◦ and minimum angle 11.11◦.

Torus interface

Figure: The mesh near the interface. Use 2.991 second to generate a mesh with 1, 045, 061 points.

Outline

Elliptic Interface Problems Interface-fitted Mesh Generation 2D Algorithm 3D Algorithm VEM for Elliptic Interface Problems Weak Formulation Virtual Element Method Weak Galerkin Methods Numerical Results

Elliptic Interface Problems

Consider −▽ · (β∆u) = f, in Ω\Γ (2) with prescribed jump conditions across the interface Γ: [u]Γ = u+ − u− = q0, [βun]Γ = β+u+

n − β−u− n = q1,

and boundary condition u = g

• n ∂Ω.

Weak Formulation: Case 1

For the elliptic interface problem, if the jump conditions are homogeneous, i.e., q0 = 0 and q1 = 0 on Γ, then the problem is equivalent to that of finding u ∈ H1(Ω) with u = g on ∂Ω such that

• Ω

β▽u · ▽vdx =

• Ω

fv dx, ∀v ∈ H1

0(Ω)

(3)

Weak Formulation: Case 2

If q0 = 0, q1 = 0 on Γ, then the corresponding weak formulation is : (β▽u, ▽v)Ω = (f, v)Ω − q1, vΓ, ∀v ∈ H1

0(Ω).

(4) The jump condition [βun]Γ = q1 holds in H−1/2(Γ) sense.

Weak Formulation: Case 3

If q0 = 0 on Γ, we can find w− : Ω− → R with w− = q0 on ∂Ω− and w− ∈ H1(Ω−). The zero extension of w− ,denote by w satisfies w ∈ H1(Ω− ∪ Ω+) with w = w− on Ω− and w = 0 on Ω+. The choice of w is not unique. The interface problem is equivalent to: find u = p − w with p ∈ H1

0(Ω) such that:

(β▽p, ▽v)Ω = (f, v)Ω − q1, vΓ + (β▽w, ▽v)Ω−, ∀v ∈ H1

0(Ω).

Virtual Element Method

Virtual element method is a generalization of the standard conforming FEM on polygon or polyhedra meshes.

◮ L. Beir˜

ao Da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. Basic principles of virtual element

• methods. Mathematical Models and Methods in Applied

Sciences, 23(1):199-214, 2013.

◮ B. Ahmad, A. Alsaedi, F. Brezzi, L. D. Marini, and A. Russo.

Equivalent projectors for virtual element methods. Computers and Mathematics with Applications, 66(3): 376-391, 2013.

◮ Beir˜

ao da Veiga L., Brezzi F., Marini L.D. and Russo A. The Hitchhiker’s Guide to the Virtual Element Method. Math. Models Meth. Appl. Sci. 24:1541-1573(2014).

Spaces in Virtual Element Methods

We introduce the following space on a simple polygon K Vk(K) := {v ∈ H1(K) : v|∂K ∈ Bk(∂K), ∆v ∈ Pk−2(K)}, where Pk(D) is the space of polynomials of degree ≤ k on D and conventionally P−1(D) = 0, and the boundary space Bk(∂K) := {v ∈ C0(∂K) : v|e ∈ Pk(e) for all edges e ⊂ ∂K}.

VEM k = 1 k = 2 k = 3

Weak Galerkin Methods

Given a polygon mesh Th and an integer k ≥ 1, we introduce Wh = {v = {v0, vb}, v0|K ∈ Pk(K), vb|e ∈ Pk−1(e) ∀e ⊂ K, ∀K ⊂ Th}. Define a modified weak gradient ∇K

w : Wk(K) → ∇Pk(K) as:

given v = {v0, vb} ∈ Wk(K), define ∇K

w v ∈ ∇Pk(K) such that

(∇K

w v, ∇p)K = −(v0, ∆p)K + vb, n · ∇p∂K

for all p ∈ Pk(K).

• L. Mu, J. Wang, and X. Ye. A weak Galerkin finite element method

with polynomial reduction. Journal of Computational and Applied Mathematics, 285:45–58, 2015.

Equivalence

The WG finite element is: find uh ∈ W 0

h such that

aWG(uh, v) = (f, v0) ∀v = {v0, vb} ∈ W 0

h.

where aWG(uh, vh) := (∇wuh, ∇wvh)+

• K⊂Th

χb(ub−u0|∂K)·χb(vb−v0|∂K). The non-conforming VEM is: find uh ∈ V 0

h such that

¯ aVEM

h

(uh, vh) = (f, Π0

kvh)

∀vh ∈ V 0

h ,

where ¯ aVEM

h

(u, v) := (∇Π∇

k u, ∇Π∇ k v) +

• K∈Th

χb(u − Π0

ku) · χb(v − Π0 kv).

Outline

Elliptic Interface Problems Interface-fitted Mesh Generation 2D Algorithm 3D Algorithm VEM for Elliptic Interface Problems Weak Formulation Virtual Element Method Weak Galerkin Methods Numerical Results

NUMERICAL RESULTS

Implementation and Errors

Our numerical test was based on MATLAB package iFEM.

• L. Chen. iFEM: an integrated finite element method package in
• MATLAB. 2009.

We consider the following errors uI − uhA =

• β1/2

h

▽(u−

I − u− h )2 0,Ω−

h + β1/2

h

▽(u+

I − u+ h )2 0,Ω+

h

1/2 uI − uh∞ = max(u−

I − u− h )0,∞,Ω−

h , u+

I − u+ h )0,∞,Ω+

h )

uI − uh0 =

• u−

I − u− h )2 0,Ω−

h + u+

I − u+ h )2 0,Ω+

h

1/2 where uh is the numerical approximation and uI is the nodal interpolation.

A Cube Domain with a Sphere Interface

We chose the following setting φ(x, y, z) =x2 + y2 + z2 − r2 u+ =10(x + y + z) u− =5ex2+y2+z2 + 20 where r = 0.75. The coefficient β is piecewise constant.

A Cube Domain with a Sphere Surface

Numerical Tests: β− = 1 and β+ = 1

#Dof h ||uI − uh||A ||uI − uh||∞ ||uI − uh||0 5489 0.125 0.56436 0.0665034 0.0570874 38481 0.0625 0.237781 0.0204017 0.0147826 285281 0.03125 0.0959572 0.00615071 0.00381557 2189860 0.015625 0.0365602 0.00177001 0.000991932

10

1

10

−2

10

−1

log(1/h) Error

||DuI−Duh|| C1h1.3154 ||uI−uh||∞ C2h1.7634

10

1

10

−3

10

−2

log(1/h) Error

||uI−uh|| C1h1.9494

Numerical Tests: β− = 1 and β+ = 10

#Dof h ||uI − uh||A ||uI − uh||∞ ||uI − uh||0 5489 0.125 0.548404 0.0606638 0.0416347 38481 0.0625 0.224057 0.0197568 0.0135145 285281 0.03125 0.0900665 0.00623164 0.00378079 2189860 0.015625 0.0339467 0.0019323 0.00103507

10

1

10

−2

10

−1

log(1/h) Error

||DuI−Duh|| C1h1.3356 ||uI−uh||∞ C2h1.677

10

1

10

−3

10

−2

log(1/h) Error

||uI−uh|| C1h1.7828

Numerical Tests:β− = 1 and β+ = 100

#Dof h ||uI − uh||A ||uI − uh||∞ ||uI − uh||0 5489 0.125 0.853561 0.0600918 0.0403713 38481 0.0625 0.338113 0.0204745 0.0134808 285281 0.03125 0.122999 0.00646717 0.00380721 2189860 0.015625 0.0450448 0.00203136 0.00104888

10

1

10

−2

10

−1

log(1/h) Error

||DuI−Duh|| C1h1.4191 ||uI−uh||∞ C2h1.6667

10

1

10

−3

10

−2

log(1/h) Error

||uI−uh|| C1h1.7623

Numerical Tests:β− = 1 and β+ = 1000

#Dof h ||uI − uh||A ||uI − uh||∞ ||uI − uh||0 5489 0.125 2.25987 0.0600214 0.0402502 38481 0.0625 0.885157 0.0205811 0.0134799 285281 0.03125 0.298777 0.00649235 0.00380979 2189860 0.015625 0.106891 0.00204206 0.00105082

10

1

10

−2

10

−1

10 log(1/h) Error

||DuI−Duh|| C1h1.4773 ||uI−uh||∞ C2h1.6666

10

1

10

−3

10

−2

log(1/h) Error

||uI−uh|| C1h1.7601

Numerical Tests: β− = 1 and β+ = 10000

#Dof h ||uI − uh||A ||uI − uh||∞ ||uI − uh||0 5489 0.125 6.99301 0.0600142 0.0402382 38481 0.0625 2.73468 0.0205919 0.0134798 285281 0.03125 0.91165 0.00649489 0.00380979 2189860 0.015625 0.324751 0.0020432 0.00105152

10

1

10

−2

10

−1

10 log(1/h) Error

||DuI−Duh|| C1h1.487 ||uI−uh||∞ C2h1.6666

10

1

10

−3

10

−2

log(1/h) Error

||uI−uh|| C1h1.7597

Algebraic Multigrid Solvers

#Dof β+ = 1 β+ = 10 β+ = 100 β+ = 1000 β+ = 10000 5489 8 8 8 8 8 38481 8 8 8 8 8 285281 10 10 10 10 10 2189860 11 13 12 12 11 Table: Iteration steps of MGCG using W-cycle and AMG preconditioner with fixed β− = 1 and increased β+.

Summary

Simple, Efficient and Semi-Unstructured 2-D and 3D Interface-fitted Mesh Generator.

THANK YOU FOR YOUR ATTENTION!

THANK YOU FOR YOUR ATTENTION! Supported by NSF DMS-1418934