Three-dimensional matched interface and boundary (MIB) method for treating geometric singularities

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

Abstract

This paper reports the three-dimensional (3D) generalization of our previous 2D higher-order matched interface and boundary (MIB) method for solving elliptic equations with discontinuous coefficients and non-smooth interfaces. New MIB algorithms that make use of two sets of interface jump conditions are proposed to remove the critical acute angle constraint of our earlier MIB scheme for treating interfaces with sharp geometric singularities, such as sharp edges, sharp wedges and sharp tips. The resulting 3D MIB schemes are of second-order accuracy for arbitrarily complex interfaces with sharp geometric singularities, of fourth-order accuracy for complex interfaces with moderate geometric singularities, and of sixth-order accuracy for curved smooth interfaces. A systematical procedure is introduced to make the MIB matrix optimally symmetric and banded by appropriately choosing auxiliary grid points. Consequently, the new MIB linear algebraic equations can be solved with fewer number of iterations. The proposed MIB method makes use of Cartesian grids, standard finite difference schemes, lowest order interface jump conditions and fictitious values. The interface jump conditions are enforced at each intersecting point of the interface and mesh lines to overcome the staircase phenomena in finite difference approximation. While a pair of fictitious values are determined along a mesh at a time, an iterative procedure is proposed to determine all the required fictitious values for higher-order schemes by repeatedly using the lowest order jump conditions. A variety of MIB techniques are developed to overcome geometric constraints. The essential strategy of the MIB method is to locally reduce a 2D or a 3D interface problem into 1D-like ones. The proposed MIB method is extensively validated in terms of the order of accuracy, the speed of convergence, the number of iterations and CPU time. Numerical experiments are carried out to complex interfaces, including the molecular surfaces of a protein, a missile interface, and van der Waals surfaces of intersecting spheres.

Introduction

The numerical solution of elliptic equations with discontinuous coefficients and singular sources has attracted much attention in the applied mathematics and scientific computing [5], [6], [7], [13], [16], [21], [26], [27], [28], [29], [31], [32], [34], [36], [41], [42], [43], [46], [61], [62], [63], [65], [66], [68] since Peskin’s pioneer work in 1977 [57]. This class of problems have demonstrated a growing impact to fluid dynamics [14], [18], [19], [30], [53], electromagnetics [23], [24], [33], material science [25], [27] and biologic systems [20], [44], [70], [75]. Due to the lack of regularity in the solution, ordinary numerical schemes, such as finite difference methods, finite element methods and spectral methods, cannot maintain the designed order of convergence and even diverge in the worst scenario. It is very difficult to construct highly accurate and efficient numerical methods for this class of problems. Apart from Peskin’s immersed boundary method (IBM) [22], [35], [57], [58], [59], a number of other elegant methods have been proposed in the literature. Among them, the ghost fluid method (GFM) proposed by Fedkiw, Osher and coworkers [15], [45] is a relatively simple and easy to use sharp interface approach. Finite element formulations [2], [8], [40], [44] are more suitable for weak solutions. An upwinding embedded boundary method was proposed by Cai and Deng [9] for dielectric media. A second-order finite volume based method was proposed Oevermann and Klein [54], which uses bilinear ansatz functions on Cartesian grids and compact stencils. A relevant, while quite distinct approach is the integral equation method for embedding complex geometry [49], [50], [51]. A major advance in the field was due to LeVeque and Li [37], who proposed a remarkable second-order sharp interface scheme, the immersed interface method (IIM). Their method has been improved for computational efficiency and matrix optimization over the past decade [1], [12], [38], [39], [62] and widely applied in practical problems. A new trend in this field is the construction of higher-order (interface) methods [21], [50], [63], [71], [73], [74] that are particularly desirable for problems involving both material interfaces and high frequency waves where conventional local adaptive refinement approaches do not work well. Typical examples are the interaction of turbulence and shock, and high frequency wave propagation in inhomogeneous media. Simple Cartesian grids are preferred in these situations because of the bypass of the mesh generation, better temporal stability and the availability of fast algebraic solvers. Recently, we have proposed a systematic higher-order method, the matched interface and boundary (MIB) method, for treating electromagnetic wave propagation and scattering in dielectric media [71]. In particular, we proposed a way to avoid the use of high order jump conditions, which require large stencils and are numerically unstable in higher-order schemes as analyzed in our work. We found that by repeatedly using only the lowest order jump conditions, an arbitrarily higher-order interface scheme can be constructed in principle [71]. This idea was applied to the fourth-order beam equation with free boundary conditions [72]. More recently, we have generalized the earlier MIB method for solving elliptic equations with curved interfaces [73], [74]. The MIB method makes use of simple Cartesian grids, standard finite difference schemes, lowest order physical jump conditions and fictitious domains. The physical jump conditions are enforced at each intersecting point of the interface and the mesh lines. Because of the use of the subgrid information at the intersecting points, the MIB method overcomes the staircase phenomena of the finite difference approximation for complex geometry and can be of arbitrarily higher-order convergence in principle when the lowest order physical jump conditions are repeatedly enforced. MIB schemes of up to 16th-order accuracy have been constructed for straight interfaces [71], [73]. For curved interfaces, fourth- and sixth-order MIB schemes have been demonstrated [73], [74]. Most recently, we have proposed an interpolation formulation of the MIB method that does not need the use of the fictitious values [74]. We have shown that our interpolation formulation of the MIB method is equivalent to our earlier fictitious domain formulation [74]. We have constructed a three-dimensional (3D) MIB scheme that is of second-order convergence for complex interfaces free of geometric singularities [75].

While much attention was focused on elliptic equations with smooth interfaces, little work has been done on the mathematical treatment of geometric singularities, such as sharp-edged, sharp-wedged and sharp-tipped interfaces. This class of problems are ubiquitous in science and engineering, such as wave-guides analysis [52], electromagnetic wave scattering and propagation [10], [55], [56], friction modeling [64], plasma-surface interaction [47] and turbulent-flow [4]. The presence of geometric singularities, such as cusps and self-intersecting surfaces, is a major obstacle to the accuracy, convergence and stability of the numerical solution of the Poisson–Boltzmann equation for the electrostatic analysis of biomolecules [70]. Technically, it is much more challenging to construct efficient numerical methods for this class of problems due to the lack of regularity in the interface. Because the gradient near the tip of a sharp-edged interface is not well defined in the classical sense, the aforementioned interface methods might not maintain the designed accuracy. A standard approach to these problems is the use of local adaptive approaches [52], [56], which exhibit reduced convergence rates and require dramatic mesh refinements in the vicinity of geometric singularities [11]. These approaches fail to work if the solution is highly oscillatory due to the so called pollution effect [3]. A remarkable result, due to Hou and Liu [26], was of 0.8th-order convergence obtained with a 2D finite element formulation for elliptic interface problems with geometric singularities. Their method is convenient for geometric singularity induced weak solutions. Recently, we have developed a 2D second-order MIB treatment of sharp-edged interfaces [69]. Unlike the Galerkin formulation, our MIB scheme does not specifically handle the geometric singularity induced weak solutions. A standard approach that multiplies the solution with an appropriate polynomial [56] was incorporated in our MIB method to deal with geometry induced weak solutions [69]. The MIB is found to be insensitive to the location of tips. The domain extension on the tip can be computed from any direction by extrapolation without affecting the result [69]. An explanation is that although gradients on the tip extrapolated from different directions contradict each other, different sets of interface conditions from different directions near the tip are consistent with each other for a given domain extension.

Another challenge that is encountered in interface schemes, particularly in higher-order ones, is the asymmetric discretization matrices which render convergence problems in solving the resulting linear algebraic equations. This difficulty is deteriorated by complex interface geometry—the matrices usually become more asymmetric under complex geometric constraints. Moreover, the higher convergence order, the larger matrix bandwidth will be for an interface scheme. Furthermore, discretization matrices are usually not diagonally dominant due to the interface treatment. This causes slow convergence as well. A discrete maximum preserving scheme was proposed to accelerate the convergence of the IIM [12]. In our earlier MIB method, a large number of iterations is required to solve the linear algebraic equations due to the asymmetric matrix.

The objective of the present work is threefold. First, we extend our earlier 2D higher-order MIB schemes [71], [73] to 3D. We report some first known sixth-order 3D results for elliptic problems with smooth interfaces. Such schemes are particularly efficient for problems involving both material interfaces and high frequency waves. Second, we generalize our earlier 2D second-order MIB method for geometric singularities [69] to 3D. Our earlier 3D MIB scheme does not maintain the designed second accuracy for molecular surfaces of proteins due to the presence of non-smooth interfaces [75]. In this work, we demonstrate the second-order convergence for arbitrarily complex interfaces with arbitrary geometric singularities. In fact, there is a critical acute angle restriction in our earlier 2D MIB scheme for sharp-edged interfaces [69]. A new algorithm, which makes use of two interface intersecting points for sharp-edged interfaces, is proposed to reserve this problem. We demonstrate that 3D fourth-order MIB schemes can be constructed for complex interfaces and interfaces with moderate geometric singularities. Finally, due to the constraint of complex interface geometry, it is impossible to create a high order MIB scheme with a symmetric discretization matrix. In this work, we have addressed this problem by making the MIB matrix optimally symmetric and banded. Here, by ‘optimally symmetric and banded’ we mean as symmetric and banded as possible under a given geometric constraint. The matrix optimization is achieved by an appropriate selection of auxiliary grid points. This optimization procedure dramatically improves the convergence property of the present MIB method. An essential feature of the proposed MIB method is that it locally simplifies a 2D or a 3D interface problem into 1D-like ones. As a result, the MIB method is able to systematically achieve higher-order convergence and to efficiently deal with geometric singularities. The proposed MIB algorithm, particularly its matrix acceleration and treatment of singularities, provides the crucial technical foundation to a new generation of Poisson–Boltzmann equation solvers for the electrostatic analysis of biomolecules [70].

In the next two sections, the theoretical formulation and the computational algorithm are given to the 3D MIB method for elliptic interface problems. Detailed consideration of topological variations is given, particularly for 3D geometric singularities. New algorithm is proposed for the sharp-edged interfaces and tips. The MIB matrix is made optimally symmetric and banded. A systematic procedure is introduced to construct 3D higher-order MIB schemes. The proposed MIB method is extensively validated by benchmark tests, such as the interfaces of missile, protein, flower and intersecting spheres. This paper ends with a conclusion to summarize the main points.

Section snippets

General setting of the MIB method for elliptic interface problems

Consider a given interface Γ which divides an open bounded domain ΩR3 into disjoint open subdomains, domain Ω and domain Ω+, i.e., Ω = Ω  Ω+  Γ. We solve the following 3D elliptic equation·(β(x,y,z)u(x,y,z))-κ(x,y,z)u(x,y,z)=q(x,y,z),x,y,zΩwith variable coefficients β(x, y, z) and κ(x, y, z) which may admit jumps at the interface Γ. To make the problem well-posed, this equation should be solved with appropriate boundary conditions and two interface jump conditions[u]=u+-u-,[βun]=β+u+·n-β-u-·n,

Results and discussion

In this section, we examine the convergence order, test the accuracy, validate the matrix optimization and demonstrate the capability of the proposed MIB II schemes. The impact of the proposed algorithm acceleration is studied with the molecular surfaces of a diatom and a protein. Molecular surfaces are generated with the MSMS program [60] with probe radius 1.4 and density 10. Reported CPU time is recorded on an AMD Turion 64™ MT-30 laptop with 1.6 GHz clock speed. The ability of the proposed

Conclusion

This paper introduces the three-dimensional (3D) higher-order matched interface and boundary (MIB) method for solving elliptic equations with discontinuous coefficients and geometric singularities. Our previous higher-order MIB methods in 1D [71] and 2D [73], [74] have been generalized to 3D. The present MIB method is of fourth-order convergence for complex interfaces with moderate geometric singularities, and of six-order convergence for smooth interfaces. Moreover, a recently proposed 2D

Acknowledgements

This work was supported in part by NSF Grants DMS-0616704 and IIS-0430987. The authors thank Shan Zhao, Yongcheng Zhou and Weihua Geng for useful discussions.

References (75)

  • T.Y. Hou et al.

    A hybrid method for moving interface problems with application to the Hele–Shaw flow

    J. Comput. Phys.

    (1997)
  • J.K. Hunter et al.

    Reactive autophobic spreading of drops

    J. Comput. Phys.

    (2002)
  • S. Jin et al.

    Robust numerical simulation of porosity evolution in chemical vapor infiltration II. Two-dimensional anisotropic fronts

    J. Comput. Phys.

    (2002)
  • H. Johansen et al.

    A Cartesian grid embedding boundary method for Poisson’s equation on irregular domains

    J. Comput. Phys.

    (1998)
  • M.C. Lai et al.

    An immersed boundary method with formal second-order accuracy and reduced numerical viscosity

    J. Comput. Phys.

    (2000)
  • M.N. Linnick et al.

    A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains

    J. Comput. Phys.

    (2005)
  • W.K. Liu et al.

    Immersed finite element method and its applications to biological systems

    Comput. Methods Appl. Mech. Eng.

    (2006)
  • X.D. Liu et al.

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

    J. Comput. Phys.

    (2000)
  • A. Mayo et al.

    Fourth order accurate evaluation of integrals in potential theory on exterior 3D regions

    J. Comput. Phys.

    (2007)
  • A. McKenney et al.

    A fast Poisson solver for complex geometries

    J. Comput. Phys.

    (1995)
  • M. Oevermann et al.

    A Cartesian grid finite volume method for elliptic equations with variable coefficients and embedded interfaces

    J. Comput. Phys.

    (2006)
  • C.S. Peskin

    Numerical analysis of blood flow in heart

    J. Comput. Phys.

    (1977)
  • C.S. Peskin et al.

    A 3-dimensional computational method for blood-flow in the heart. 1. Immersed elastic fibers in a viscous incompressible fluid

    J. Comput. Phys.

    (1989)
  • M. Schulz et al.

    Two-dimensional modelling of the river Rhine

    J. Comput. Appl. Math.

    (2002)
  • J.A. Sethian et al.

    Structural boundary design via level set and immersed interface methods

    J. Comput. Phys.

    (2000)
  • A.K. Tornberg et al.

    Numerical approximations of singular source terms in differential equations

    J. Comput. Phys.

    (2004)
  • S.N. Yu et al.

    Matched interface and boundary (MIB) method for elliptic problems with sharp-edged interfaces

    J. Comput. Phys.

    (2007)
  • S. Zhao et al.

    High order FDTD methods via derivative matching for Maxwell’s equations with material interfaces

    J. Comput. Phys.

    (2004)
  • S. Zhao et al.

    DSC analysis of free-edged beams by an iteratively matched boundary method

    J. Sound Vib.

    (2005)
  • Y.C. Zhou et al.

    High order matched interface and boundary (MIB) schemes for elliptic equations with discontinuous coefficients and singular sources

    J. Comput. Phys.

    (2006)
  • Y.C. Zhou et al.

    On the fictitious-domain and interpolation formulations of the matched interface and boundary (MIB) method

    J. Comput. Phys.

    (2006)
  • L. Adams et al.

    The immersed interface/multigrid methods for interface problems

    SIAM J. Sci. Comput.

    (2002)
  • I. Babuška

    The finite element method for elliptic equations with discontinuous coefficients

    Computing

    (1970)
  • I. Babuška et al.

    Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave number?

    SIAM J. Numer. Anal.

    (1997)
  • H. Ben Ameur et al.

    Level set methods for geometric inverse problems in linear elasticity

    Inverse Probl.

    (2004)
  • J. Bramble et al.

    A finite element method for interface problems in domains with smooth boundaries and interfaces

    Adv. Comput. Math.

    (1996)
  • S. Caorsi et al.

    Electromagnetic scattering by a conducting strip with a multilayer elliptic dielectric coating

    IEEE T. Electromag. Compat.

    (1999)
  • Cited by (0)

    View full text