Numerical method for solving matrix coefficient elliptic equation with sharp-edged interfaces

doi:10.1016/j.jcp.2010.06.005

Journal of Computational Physics

Volume 229, Issue 19, 20 September 2010, Pages 7162-7179

https://doi.org/10.1016/j.jcp.2010.06.005 Get rights and content

Abstract

Solving elliptic equations with sharp-edged interfaces is a challenging problem for most existing methods, especially when the solution is highly oscillatory. Nonetheless, it has wide applications in engineering and science. An accurate and efficient method is desired. We propose an efficient non-traditional finite element method with non-body-fitting grids to solve the matrix coefficient elliptic equations with sharp-edged interfaces. Extensive numerical experiments show that this method is second order accurate in the L^∞ norm and that it can handle both sharp-edged interface and oscillatory solutions.

Introduction

The importance of elliptic interface problems has been well recognized in a variety of disciplines, such as electromagnetics, material science, fluid dynamics and so on. However, designing highly efficient methods for these problems is a difficult job because of the low global regularity of the solution. Since 1977, after the pioneering work of Peskin [14], much attention has been paid to the numerical solution of elliptic equations with discontinuous coefficients and singular sources on regular Cartesian grids. In many studies, simple Cartesian grids are preferred. In this way, the complicated procedure of generating an unstructured grid can be bypassed, and we can use well developed fast algebraic solvers.

The “immersed boundary” method [14], [15] uses a numerical approximation of the δ-function, which smears out the solution on a thin finite band around the interface Γ. In [16], the “immersed boundary” method was combined with the level-set method, resulting in a first order numerical method that is simple to implement, even in multiple spatial dimensions. However, for both methods, the numerical smearing at the interface forces continuity of the solution at the interface, regardless of the interface condition [u] = a, where a might not be zero.

In [11], [12], a fast method for solving Laplace’s equations on irregular regions with smooth boundaries was introduced. By using Fredholm integral equations of the second kind, solutions can be extended to a rectangular region. Since solutions are harmonic, Fredholm integral equations can be used to capture the jump conditions in the solution and its normal derivative, [u] ≠ 0 and [u_n] = 0. Then these jump conditions are used to evaluate the discrete Laplacian, and then a fast Poisson solver on a regular region can be applied with second or higher order accuracy.

In [5], the “immersed interface” method was presented. This method achieves second order accuracy by incorporating the interface conditions into the finite difference stencil in a way that preserves the interface conditions in both solution and its co-normal derivative, [u] ≠ 0 and [βu_n] ≠ 0. The corresponding linear system is sparse, but not symmetric or positive definite. Various applications and extensions of the “immersed interface” method are discussed in [8].

In [20], a finite element method was developed for solving elliptic problems with the interface conditions [u] = 0 and [βu_n] ≠ 0. Interfaces are aligned with cell boundaries. Second order accuracy was obtained in an energy norm. Nearly second order accuracy was obtained in the L² norm.

In [6], a fast iterative method in conjunction with the “immersed interface” method has been developed for constant coefficient problems with the interface conditions [u] = 0 and [βu_n] ≠ 0. Non-body-fitting Cartesian grids are used, and then associated uniform triangulations are added on. Interfaces are not necessarily aligned with cell boundaries. Numerical evidence shows that this method’s conforming version achieves second order accuracy in the L^∞ norm, and higher than first order for its non-conforming version.

The boundary condition capturing method [9] uses the Ghost fluid method (GFM) [2] to capture the boundary conditions. The GFM is robust and simple to implement, and so is the resulting boundary condition capturing method. The boundary condition capturing method has been sped up by a multi-grid method [17]. The convergence proof for the method is provided in [10]. The boundary condition capturing method can be obtained from discretizing the weak formulation provided in [10]. The convergence proof follows naturally. The method can solve the elliptic equation with interface conditions [u] ≠ 0 and [βu_n] ≠ 0 in two and three dimensions. However, the original version is only first order accurate. In a recent work [25], the method is improved to second order accuracy for smooth interfaces.

In [4], a non-traditional finite element formulation for solving elliptic equations with smooth or sharp-edged interfaces was proposed with non-body-fitting grids for [u] ≠ 0 and [βu_n] ≠ 0. It achieved second order accuracy in the L^∞ norm for smooth interfaces and about 0.8th order for sharp-edged interfaces. In [21], the matched interface and boundary (MIB) method was proposed to solve elliptic equations with smooth interfaces. In [19], the MIB method was generalized to treat sharp-edged interfaces. With an elegant treatment, second order accuracy was achieved in the L^∞ norm. However, for oscillatory solutions, the errors degenerate. Also, there has been a large body of work from the finite volume perspective for developing high order methods for elliptic equations in complex domains, such as [22], [23] for two dimensional problems and [24] for three dimensional problems. Another recent work in this area is a class of kernel-free boundary integral (KFBI) methods for solving elliptic BVPs, presented in [18].

In this paper, inspired by the boundary condition capturing method [9] and the weak formulation derived in [10], we further generalize the method introduced in [4]. We use a finite element formulation for solving elliptic equations with sharp-edged interfaces with β being uniformly elliptic (therefore, positive definite) and lower order terms present. We provided proofs for the generalized version of the theorems in [10], and we proved the resulting linear system is (unsymmetric) positive definite if β is positive definite and lower order terms are not present. We also provided extensive numerical experiments. Compared with the previous work in [4], we improved the order of accuracy for sharp-edged interface from 0.8th to close to second order, see Example 4. Compared with the results in [19], the more oscillatory the solution is, the more advantageous our method is, see Example 1, Example 2, Example 3. The orders of accuracy for different regularity of solutions and different regularity of the interface are listed in Table 11.

Section snippets

Equations and weak formulation

Consider an open bounded domain Ω ⊂ R^d. Let Γ be an interface of co-dimension d − 1, which divides Ω into disjoint open subdomains, Ω⁻ and Ω⁺, hence Ω = Ω⁻ ⋃ Ω⁺ ⋃ Γ. Assume that the boundary ∂Ω and the boundary of each subdomain ∂Ω^± are Lipschitz continuous as submanifolds. Since ∂Ω^± are Lipschitz continuous, so is Γ. A unit normal vector of Γ can be defined a.e. on Γ, see Section 1.5 in [3].

We seek solutions of the variable coefficient elliptic equation away from the interface Γ given by $- ▽ \cdot (β (x) ▽ u (x)) + p ($

Numerical method

For ease of discussion in this section, and for accuracy testing in the next section, we assume that a, b and c are smooth on the closure of Ω, β and f are smooth on Ω⁺ and Ω⁻, but they may be discontinuous across the interface Γ. However, ∂Ω, ∂Ω⁻ and ∂Ω⁺ are kept to be Lipschitz continuous. We assume that there is a Lipschitz continuous and piecewise smooth level-set function ϕ on Ω, where Γ = {ϕ = 0}, Ω⁻ = {ϕ < 0} and Ω⁺ = {ϕ > 0}. A unit vector $n = \frac{▽ ϕ}{| ▽ ϕ |}$ can be obtained on $\bar{Ω}$ , which is a unit normal

Numerical experiments

Consider the problem $\begin{matrix} - \nabla \cdot (β \nabla u) + p \cdot \nabla u + qu = f,in Ω^{\pm}, \\ [u] = a,on Γ, \\ [(β \nabla u) \cdot n] = b,on Γ, \\ u = g,on \partial Ω, \end{matrix}$ on the rectangular domain Ω = (x_min, x_max) × (y_min, y_max). Γ is an interface and prescribed by the zero level-set {(x, y) ∈ Ω ∣ ϕ(x, y) = 0} of a level-set function ϕ(x,y). The unit normal vector of Γ is $n = \frac{\nabla ϕ}{| \nabla ϕ |}$ pointing from Ω⁻ = {(x, y) ∈ Ω ∣ ϕ(x, y) ⩽ 0} to Ω⁺ = {(x, y) ∈ Ω ∣ ϕ(x, y) ⩾ 0}.

In all numerical experiments below, the level-set function ϕ(x, y), the coefficients β^±(x, y), p^±(x, y), q^±(x, y) and the solutions $\begin{matrix} u = u^{+} (x, y), in Ω^{+}, \\ u = u^{-} (x, y), \end{matrix}$

Conclusion

In this paper, we generalized the previous work in [4] to solve matrix coefficient second order elliptic equations with lower order terms present (see Example 11) for interface problems. We provided proofs for the generalized version of theorems in [4]. We also proved that the matrix for the linear system generated by our method is positive definite (but not symmetric). Through numerical experiments, our method achieved second order accuracy in the L^∞ norm, and we can handle the difficulties of

Acknowledgement

This research is partially supported by Louisiana Board of Regents RCS Grant No. LEQSF (2008–11)-RD-A-18.

References (25)

Songming Hou et al.
A numerical method for solving variable coefficient elliptic equations with interfaces
J. Comput. Phys.
(2005)
W.-J. Ying et al.
A kernel-free boundary integral method for elliptic boundary value problems
J. Comput. Phys.
(2007)
Y.C. Zhou et al.
High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources
J. Comput. Phys.
(2006)
Lawrence C. Evans
Partial Differential Equations
(1998)
R. Fedkiw et al.
A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)
J. Comput. Phys.
(1999)
P. Grisvard, Elliptic problems in nonsmooth domains – monographs and studies in mathematics. Pitman Advanced Publishing...
R.J. LeVeque et al.
The immersed interface method for elliptic equations with discontinuous coefficients and singular sources
SIAM J. Numer. Anal.
(1994)
Z. Li
A fast iterative algorithm for elliptic interface problems
SIAM J. Numer. Anal.
(1998)
Z. Li, T. Lin, X. Wu, New cartesian grid methods for interface problems using the finite element formulation,...
Zhilin Li et al.
The Immersed Interface Method: Numerical Solutions of Pdes Involving Interfaces and Irregular Domains
(2006)

Xu-Dong Liu et al.

A boundary condition capturing method for Poisson’s equation on irregular domains

J. Comput. Phys.

(2000)

Xu-Dong Liu et al.

Convergence of the ghost fluid method for elliptic equations with interfaces

Math. Comp.

(2003)

Cited by (82)

A Petrov–Galerkin immersed finite element method for steady Navier–Stokes interface problem with non-homogeneous jump conditions
2024, Journal of Computational and Applied Mathematics
This paper develops a lowest-order Petrov–Galerkin immersed $Q_{1} - Q_{0}$ method for solving Stokes interface problem. We make use of a uniform, interface-unfitted Cartesian mesh. An immersed Petrov–Galerkin formulation is presented, where the test spaces are conventional finite element spaces and the solution spaces satisfying the jump conditions. For the Stokes and Navier–Stokes interface problem, simple stabilized items are introduced. The nonlinear convective term is treated using Picard’s iteration. Extensive numerical experiments validate the feasibility and optimal convergence order for the Petrov–Galerkin immersed $Q_{1} - Q_{0}$ scheme both with homogeneous and non-homogeneous jumps.
A discrete fracture-matrix approach based on Petrov-Galerkin immersed finite element for fractured porous media flow on nonconforming mesh
2024, Journal of Computational Physics
In this paper, the Petrov-Galerkin immersed finite element method is developed for simulating flow in fractured porous media. We consider the hybrid-dimensional fracture model as an elliptic interface problem with multiple physical coupling interface conditions, and embed the information of low-dimensional fractures into the local functions of interface elements. The immersed finite element space is constructed by explicitly representing the local solution functions in the interface element. This approach addresses the limitation of conforming meshes. More importantly, it can capture both the continuity and discontinuity of pressure, making it effective for both conductive and blocking fractures. We provide a detailed illustration of how to construct local functions on triangular and rectangular meshes. Finally, numerical examples and benchmarks are presented to validate the effectiveness of our proposed method.
Error analysis of Petrov-Galerkin immersed finite element methods
2023, Computer Methods in Applied Mechanics and Engineering
This paper designs and analyzes a new and stable Petrov–Galerkin (PG) immersed finite element method (IFEM) for the second-order elliptic interface problems by introducing stabilization terms based on the classical PG-IFEM, which lacks the local positivity. The Petrov–Galerkin immersed finite element method uses the immersed finite element functions for the trial space and the standard finite element functions for the test space. Both the a priori and a posteriori error estimates of the method are analyzed in this paper. We prove the continuity and inf-sup condition and the a priori error estimate of the energy norm. The proposed a posteriori error estimator is proved to be both reliable and efficient, with both reliability and efficiency constants independent of the location of the interface. Extensive numerical results confirm the numerical scheme’s optimal convergence and indicate the robustness with respect to the interface-mesh intersection and the coefficient contrast, despite the robustness of the inf-sup constant with respect to the interface-mesh intersection has yet been theoretically proved.
A reconstructed discontinuous approximation on unfitted meshes to H(curl) and H(div) interface problems
2023, Computer Methods in Applied Mechanics and Engineering
We propose a high-order discontinuous Galerkin finite element method for solving the $H (div)$ - and $H (curl)$ -elliptic interface problems on unfitted meshes. The vector-valued approximation space is constructed by the patch reconstruction with at most $d$ degrees of freedom per element. The $C^{2}$ -smooth interface is allowed to intersect elements in a quite general setup. The patch reconstruction provides the stability near the interface naturally without any additional stabilization strategy. The method is based on the symmetric interior penalty method and the optimal and suboptimal convergence rates under the energy norm and $L^{2}$ norm are derived. Numerical examples in both two and three dimensions are presented to demonstrate the accuracy of the method.
Optimal error bound for immersed weak Galerkin finite element method for elliptic interface problems
2022, Journal of Computational and Applied Mathematics
Citation Excerpt :
Later Adjerid and Lin [43–48] developed the IFE functions with an arbitrary polynomial degree. Recently by imitating IFEM, which combines the idea of immersed methods with the classical FEM, the immersed idea has been combined with non-classical FEM and finite volume method, such as immersed Petrov–Galerkin method [49,50], nonconforming IFEM [51–53], immersed discontinuous Galerkin method [54,55], partially penalized immersed finite element method (PP-IFEM) [47,54,56–58] and immersed finite volume method [59–61]. Lately, the WG-FEM has been presented for elliptic interface problems [64,72–75].
In Mu and Zhang (2019), an immersed weak Galerkin finite element method (IWG-FEM) is developed for solving elliptic interface problems and it is proved that this method has optimal a-priori error estimate in an energy norm under artificial smoothness assumption on the solution. In this study, we prove that IWG-FEM converges optimally in energy norm under natural smoothness assumption on solution. Furthermore, we show that IWG-FEM converges optimally in the $L^{2}$ norm which did not present in Mu and Zhang (2019) because of the artificial $H^{3}$ smoothness requirement. A series of numerical experiments are conducted and reported to verify the theoretical finding.
A new primal-dual weak Galerkin method for elliptic interface problems with low regularity assumptions
2022, Journal of Computational Physics
Citation Excerpt :
Even though successes have been achieved in the endeavor of solving elliptic interface problems, challenges remain in the search of new and efficient numerical algorithms for problems with very complicated interface geometries, and for problems with low-regularity solutions. The low-regularity of the solution is often caused by the geometric singularities of the interfaces and/or the non-smoothness of the interface data [37,26]. The paper is organized as follows.
This article introduces a new primal-dual weak Galerkin (PDWG) finite element method for second order elliptic interface problems with ultra-low regularity assumptions on the exact solution and the interface and boundary data. It is proved that the PDWG method is stable and accurate with optimal order of error estimates in discrete and Sobolev norms. In particular, the error estimates are derived under the low regularity assumption of $u \in H^{δ} (Ω)$ for $δ > \frac{1}{2}$ for the exact solution u. Extensive numerical experiments are conducted to provide numerical solutions that verify the efficiency and accuracy of the new PDWG method.

View all citing articles on Scopus

View full text

Published by Elsevier Inc.

Numerical method for solving matrix coefficient elliptic equation with sharp-edged interfaces

Abstract

Introduction

Section snippets

Equations and weak formulation

Numerical method

Numerical experiments

Conclusion

Acknowledgement

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

Partial Differential Equations

A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)

J. Comput. Phys.

The immersed interface method for elliptic equations with discontinuous coefficients and singular sources

SIAM J. Numer. Anal.

A fast iterative algorithm for elliptic interface problems

SIAM J. Numer. Anal.

The Immersed Interface Method: Numerical Solutions of Pdes Involving Interfaces and Irregular Domains

A boundary condition capturing method for Poisson’s equation on irregular domains

J. Comput. Phys.

Convergence of the ghost fluid method for elliptic equations with interfaces

Math. Comp.