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 − = q 0 , [ βu n ] Γ = β + u + n − β − u − n = q 1 , and boundary condition u = g on ∂ Ω .

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) of 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 p 6 p 7 p 5 p 3 p 1 p 4 p 2 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 (b) Step 2: Form the (c) Step 3: Project to the paraboloid lowest convex hull in the lowest convex hull R n +1 . to R n . Figure: Delaunay triangulation is the projection of the lower convex hull of points lifted to the paraboloid f = � x � 2 .

Interface Recovery ����������������� ����� ��� The lower convex hull when lift to R n +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 − = q 0 , [ βu n ] Γ = β + u + n − β − u − n = q 1 , and boundary condition u = g on ∂ Ω .

Weak Formulation: Case 1 For the elliptic interface problem, if the jump conditions are homogeneous, i.e., q 0 = 0 and q 1 = 0 on Γ, then the problem is equivalent to that of finding u ∈ H 1 (Ω) with u = g on ∂ Ω such that � � fv dx, ∀ v ∈ H 1 β ▽ u · ▽ vdx = 0 (Ω) (3) Ω Ω

Weak Formulation: Case 2 If q 0 = 0, q 1 � = 0 on Γ, then the corresponding weak formulation is : ∀ v ∈ H 1 ( β ▽ u, ▽ v ) Ω = ( f, v ) Ω − � q 1 , v � Γ , 0 (Ω) . (4) The jump condition [ βu n ] Γ = q 1 holds in H − 1 / 2 (Γ) sense.