Research Article  Open Access
Liqun Wang, Liwei Shi, "Numerical Method for Solving Matrix Coefficient Elliptic Equation on Irregular Domains with SharpEdged Boundaries", Advances in Numerical Analysis, vol. 2013, Article ID 967342, 10 pages, 2013. https://doi.org/10.1155/2013/967342
Numerical Method for Solving Matrix Coefficient Elliptic Equation on Irregular Domains with SharpEdged Boundaries
Abstract
We present a new secondorder accurate numerical method for solving matrix coefficient elliptic equation on irregular domains with sharpedged boundaries. Nontraditional finite element method with nonbodyfitting grids is implemented on a fictitious domain in which the irregular domains are embedded. First we set the function and coefficient in the fictitious part, and the nonsmooth boundary is then treated as an interface. The emphasis is on the construction of jump conditions on the interface; a special position for the ghost point is chosen so that the method is more accurate. The test function basis is chosen to be the standard finite element basis independent of the interface, and the solution basis is chosen to be piecewise linear satisfying the jump conditions across the interface. This is an efficient method for dealing with elliptic equations in irregular domains with nonsmooth boundaries, and it is able to treat the general case of matrix coefficient. The complexity and computational expense in mesh generation is highly decreased, especially for moving boundaries, while robustness, efficiency, and accuracy are promised. Extensive numerical experiments indicate that this method is secondorder accurate in the norm in both two and three dimensions and numerically very stable.
1. Introduction
Let () be an openbounded domain with a Lipschitz continuous boundary .
We consider the variable coefficient elliptic equation where refers to the spatial variable, is the gradient operator, and the righthand side is assumed to lie in . The coefficient is a matrix that is uniformly elliptic, and its entries are continuously differentiable on .
For a given function on the boundary , the Dirichlet boundary condition is prescribed as
Elliptic partial differential equations are often used to construct models of the most basic theories underlying physics and engineering, such as electromagnetism, material science, and fluid dynamics. Different kinds of boundary conditions arise with the wide range of applications, such as Dirichlet boundary condition, Neumann boundary condition, and Robin boundary condition.
Elliptic equation on irregular domains has been studied by many researchers and several techniques have been developed. Finite element methods use a mesh triangulation to capture the boundary [1–4]. However, in many situations, such as when the boundary is moving, the mesh generation may be both computational expensive and challenging. A more preferred method is to combine the Cartesian grid method with levelset approach [5–7] to capture the boundary. Higherorder accuracy can be achieved by modifying standard difference formulas; examples are ShortleyWeller discretization [8] and the fast solution proposed by Mayo in [9] for solving Poisson or biharmonic equation on irregular regions with smooth boundary.
Since the pioneer work of Peskin [10] in 1977, much attention has been paid to the numerical solution of elliptic equations with discontinuous coefficients and singular sources on regular Cartesian grids. A large class of finite difference methods have been proposed. The main idea is to use difference scheme and stencils carefully near the interface to incorporate jump conditions and achieve highorder local truncation error using Taylor expansion. The immersed interface method [10, 11] incorporates the interface conditions into the finite difference scheme near the interface to achieve secondorder accuracy based on a Taylor expansion in a local coordinate system. The boundary condition capturing method [12] uses the ghost fluid method [13] to capture the boundary conditions. The method extends the solution from one side across the interface using the jump conditions. In [14], the matched interface and boundary method was proposed to solve elliptic equations with smooth interfaces. In [15], the matched interface and boundary method was generalized to treat sharpedged interfaces. With an elegant treatment, secondorder accuracy was achieved in the norm. However, for oscillatory solutions, the errors degenerate.
The existing finite element schemes on nonbodyfitted meshes are usually designed by modifying the finite element basis near the interface. The extended finite element method [16–18] was proposed by Ted Belytschko and collaborators to ease difficulties in solving problems with localized features that are not efficiently resolved by mesh refinement. In immersed finite element method [19], a Lagrangian solid mesh moves on top of a background Eulerian fluid mesh which spans the entire computational domain. The penalty finite element method [20, 21] modifies the bilinear form near the interface by penalizing the jump of the solution value (with no general flux jump) across the interface.
Also, there has been a large body of work from the finite volume perspective for developing highorder methods for elliptic equations in complex domains, such as [22, 23] for twodimensional problems and [24] for threedimensional problems. Another recent work in this area is a class of kernelfree boundary integral (KFBI) methods for solving elliptic BVPs, presented in [25].
In this paper, we implemented the nontraditional finite element method into matrix coefficient elliptic equation on irregular domains. The idea of applying interface problem solvers into noninterface problems is completely new. The nontraditional finite element method is originally developed in [26–32] for solving elliptic or elasticity equations with sharpedged interfaces. This method uses nonbodyfitting Cartesian grids and uses different basis for the solution and test function; its linear system is independent of jump condition. In this way, it is easy to implement, especially for moving interface. It is capable of treating the more general case of variable matrix coefficient and can deliver (almost) secondorder accuracy in norm. Since the method is based on the assumption of discontinuous solution, in this work, we embed the research domain with irregular boundary into a fictitious rectangular/cube. After setting the function and coefficient in the fictitious part, we construct the jump conditions and calculate the integration on the boundary of the research domain (we call it interface in the following discussion). By implementing the nontraditional finite element method on the fictitious domain containing the irregular domains which are under research, we overcome the difficulty of capturing the sharpedged boundary with high accuracy. Both two and threedimensional models are studied. Extensive numerical experiments show the efficiency of this method.
The rest of the paper is organized as follows. The next section presents the discretization of two and threedimensional irregular domains. In Section 3, we construct the jump conditions and develop the weak formulation of the problem. In Section 4, six extensive numerical experiments are presented to show the secondorder accuracy of our method. Finally some conclusions are drawn in Section 5.
2. Discretization of the Domain
For simplicity, we first embed the irregular domain into a fictitious rectangular/cube ; the boundary of is denoted by , the domain between and is denoted by , and let . For ease of discussion, from now on we call the interface and the boundary of . In this paper, we restrict ourselves to a fictitious rectangular domain in twodimensional space and a fictitious cube domain in threedimensional space; see Figure 1. We assume that there is a Lipschitz continuous and piecewise smooth levelset function on , where , , and . A unit vector can be obtained on , which is a unit normal vector of at pointing from to .
(a) The rectangular domain
(b) The cube domain
2.1. Discretization of TwoDimensional Problems
In the twodimensional case, is a matrix coefficient that is uniformly elliptic in . Given two positive integers and , set and . We define a uniform Cartesian grid with grid points given by for and . A grid point is called a boundary point if or ; otherwise, it is called an interior point. The grid size is defined as .
Two sets of grid functions are needed and they are denoted by
To triangulate the domain into a set of uniform triangular cells, we divide each rectangular region into two right triangular cells: one is bounded by , and and the other by , and . Collecting all these triangular cells, we obtain a uniform triangulation of ; see Figure 2(a). Note that we can also choose the hypotenuse to be and get a different uniform triangulation from the same Cartesian grid. There is no conceptual difference for our method on these two triangulations.
(a) A schematic illustration of triangulation
(b) Interface cells
If , the grid point is counted as in ; otherwise, it is counted as in .
2.2. Discretization of ThreeDimensional Problems
For the threedimensional problem, is a matrix coefficient that is uniformly elliptic in . Given positive integers , , and , set , , and .
We define a uniform Cartesian grid with grid points given by for , , and . A grid point is called a boundary point if , , or ; otherwise it is called an interior point. The grid size is defined as .
The following two sets of grid functions are needed:
Each cube cell region is cut into six tetrahedron cells. Upon collecting all these cells, we obtain a uniform territorialization of the cube cell region ; see Figure 3.
(a) Cube cells of threedimensional problems
(b) Territorialization of threedimensional problems
If , the grid point is counted as in ; otherwise it is counted as in .
In both two and threedimensional cases, a cell belongs to one of the following two different sets:
If a cell belongs to , we call it a regular cell; otherwise we call it an interface cell, written as . and are separated by the interface segment, denoted by ; see Figures 2(b) and 4.
(a) Case 1: The interface segment is a triangle
(b) Case 2: The interface segment is a polygon
Two extension operators are needed. The first one is . For any , is a standard continuous piecewise linear function, which is a linear function in every cell and matches on grid points. Clearly such a function set, denoted by , is a finite dimensional subspace of . The second extension operator is constructed as follows. For any with at boundary points, is a piecewise linear function and matches on grid points. It is a linear function in each regular cell, just like the first extension operator in a regular cell. In each interface cell, it consists of two pieces of linear functions: one is on and the other is on .
3. Construction of Jump Conditions and Weak Formulation
In [26–32], jump conditions and are known and are used to derive the linear combination of the values of the points across the interface . However, in this paper, we need to construct and for each interface cell from the points with known value across the interface and on and then calculate the linear integral
Here is smooth on and is kept to be Lipschitz continuous. is the unit vector defined in Section 2.
For different and we have different and . For simplicity, set as 0 and on ; then . For , since , the difficult part is to generate . For accuracy consideration, we choose the position of the ghost point as the midpoint of the grid points on . Let be the value of point ; then , where is the value of the grid point / on .
For twodimensional problem, cell is cut into two parts, and ; see Figure 2(b). If part , then is chosen to be point 1, ; if part , then is chosen to be the midpoint of points and , . Since the values of points 4 and 5 are known as and , the value of can be determined by
Let
Then
For threedimensional problem, there are two kinds of interface cells.
Case 1. The interface cell is cut into a triangular pyramid and a pentahedron ; see Figure 4(a). If part , then is chosen to be point 1, ; if part , then is chosen to be the midpoint of points 2, 3, and 4, .
Case 2. The interface cell is cut into a tetrahedron and a pentahedron . We can always rotate the cell to make , such that is the midpoint of points 3 and 4, .
In Figure 4(a), the values of points 5, 6, and 7 are known as , , and ; therefore, the value of can be calculated by
In Figure 4(b), the values of points 5, 6, 7, and 8 are known as , , , and ; therefore, the value of can be calculated by
Let
Then
Lemma 1. All coefficients in equations (8)–(12) are finite and independent of and .
Proof. See Lemma 3.1 in [26].
Upon above discussion, the system described in (1)(2) is equivalent to the following system:
Since is known, we only care about the value on . We generalize the weak formulation for the elliptic equation with matrix coefficient. The usual Sobolev space is used. For , instead of the usual inner product, we choose one which is better suited to our problem:
In this way, we have the following definition.
Definition 2. A function with and on is a weak solution of system (1)(2) if it satisfies
A classical solution of system (1)(2), , is necessarily a weak solution.
Theorem 3. If , then there exists a unique weak solution of equations (1)(2).
Proof. See Theorem 2.1 in [33].
Theorem 4. For all , can be constructed uniquely, provided and are given.
Proof. See Theorem 3.1 in [26].
Method 1. Find a discrete function such that on boundary points and so that for all , we have
Theorem 5. If is positive definite, then the matrix for the linear system generated by Method 1 is positive definite.
Proof. For any vector ,
where and are basis functions for the solution and the test functions, respectively. According to Lemma 1, matrix is independent of ; choose such that everywhere. Let
Considering the positive definiteness of , we have
Therefore, is positive definite.
In our implementation, the integrals are computed with Gaussian quadrature rule. For each cell, we denote the midpoint of each edge as . In numerical computation, we apply the average of all the in each cell.
4. Numerical Experiments
Of all the numerical experiments below, in , the levelset function , the coefficients , and the solutions are given. Hence the source term and the Dirichlet boundary condition can be directly evaluated and the jump conditions and can be constructed. The norm of the error in solution is measured over the domain . The uniform Cartesian grid is restricted to the whole domain .
4.1. Numerical Examples of TwoDimensional Problems
Example 1. This example has a “star” interface. The levelset function , the coefficients , and the solution are given as follows: with , , , and
The numerical result with our method using a uniform Cartesian grid is shown in Figure 5. Table 1 shows the error on different grids.

Example 2. This example has a “happy face” interface. The levelset function , the coefficients , and the solution are given as follows:
The numerical result with our method using a uniform Cartesian grid is shown in Figure 6. Table 2 shows the error on different grids.

Example 3. This example has a “target" interface. The levelset function , the coefficients , and the solution are given as follows:
The numerical result with our method using a uniform Cartesian grid is shown in Figure 7. Table 3 shows the error on different grids.

4.2. Numerical Examples of ThreeDimensional Problems
Example 4. The levelset function , the coefficients , and the solution are given as follows:
Figure 8 shows the interface and boundary condition using a uniform Cartesian grid. Table 4 shows the error on different grids.

Example 5. The levelset function , the coefficients , and the solution are given as follows:
Figure 9 shows the interface and boundary condition using a uniform Cartesian grid. Table 5 shows the error on different grids.

Example 6. The levelset function , the coefficients , and the solution are given as follows:
Figure 10 shows the interface and boundary condition using a uniform Cartesian grid. Table 6 shows the error on different grids.

5. Conclusion
We developed a novel numerical method for solving matrix valued coefficient elliptic equation on irregular domains. The research domain with sharpedged boundary is embedded into a fictitious rectangular/cube and ghost point is used to construct the jump conditions on the interface. Nontraditional finite element method, which employs different basis for the solution and the test function, is implemented in the whole domain. Both two and threedimensional problems are discussed in this paper. Extensive numerical experiments show that this method is secondorder accurate in the norm in two and three dimensions. The work shows the robustness and efficiency of the method and indicates that this method can be used as an alternative way to solve elliptic equations; it lays a foundation for interface problems on irregular domains with various kinds of boundary conditions; as such, the variety of problems that are solvable is largely expanded.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors are grateful to the editor and referees for their valuable comments and suggestions. L. Wang’s research is supported by Science Foundation of China University of Petroleum Beijing (no. YJRC201348).
References
 A. Quarteroni, Numerical Models for Differential Problems, vol. 2, Springer, Milan, Italy, 2009. View at: Publisher Site  MathSciNet
 A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, vol. 23, Springer, Berlin, Germany, 1994. View at: MathSciNet
 R. H. Nochetto, M. Paolini, and C. Verdi, “An adaptive finite element method for twophase Stefan problems in two space dimensions. II. Implementation and numerical experiments,” SIAM Journal on Scientific and Statistical Computing, vol. 12, no. 5, pp. 1207–1244, 1991. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 A. Schmidt, “Computation of three dimensional dendrites with finite elements,” Journal of Computational Physics, vol. 125, pp. 293–312, 1996. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Sussman, P. Smereka, and S. Osher, “A level set approach for computing solutions to incompressible 2phase flow,” Journal of Computational Physics, vol. 114, pp. 146–159, 1994. View at: Publisher Site  Google Scholar
 S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, vol. 153, Springer, New York, NY, USA, 2003. View at: MathSciNet
 A. du Chéné, C. Min, and F. Gibou, “Secondorder accurate computation of curvatures in a level set framework using novel highorder reinitialization schemes,” Journal of Scientific Computing, vol. 35, no. 23, pp. 114–131, 2008. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 G. H. Shortley and R. Weller, “The numerical solution of Laplaces equation,” Journal of Applied Physics, vol. 9, pp. 334–348, 1938. View at: Publisher Site  Google Scholar
 A. Mayo, “The fast solution of Poisson's and the biharmonic equations on irregular regions,” SIAM Journal on Numerical Analysis, vol. 21, no. 2, pp. 285–299, 1984. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 C. S. Peskin, “Numerical analysis of blood flow in the heart,” Journal of Computational Physics, vol. 25, no. 3, pp. 220–252, 1977. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 C. S. Peskin and B. F. Printz, “Improved volume conservation in the computation of flows with immersed elastic boundaries,” Journal of Computational Physics, vol. 105, no. 1, pp. 33–46, 1993. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 X.D. Liu, R. P. Fedkiw, and M. Kang, “A boundary condition capturing method for Poisson's equation on irregular domains,” Journal of Computational Physics, vol. 160, no. 1, pp. 151–178, 2000. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. P. Fedkiw, T. Aslam, B. Merriman, and S. Osher, “A nonoscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method),” Journal of Computational Physics, vol. 152, no. 2, pp. 457–492, 1999. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 S. Yu, Y. Zhou, and G. W. Wei, “Matched interface and boundary (MIB) method for elliptic problems with sharpedged interfaces,” Journal of Computational Physics, vol. 224, no. 2, pp. 729–756, 2007. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 Y. C. Zhou, S. Zhao, M. Feig, and G. W. Wei, “High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources,” Journal of Computational Physics, vol. 213, no. 1, pp. 1–30, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 N. Moes, J. Dolbow, and T. Belytschko, “A finite element method for crack growth without remeshing,” International Journal for Numerical Methods in Engineering, vol. 46, no. 1, pp. 131–150, 1999. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 G. J. Wagner, N. Moës, W. K. Liu, and T. Belytschko, “The extended finite element method for rigid particles in Stokes flow,” International Journal for Numerical Methods in Engineering, vol. 51, no. 3, pp. 293–313, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 J. Chessa and T. Belytschko, “An extended finite element method for twophase fluids: flow simulation and modeling.,” Journal of Applied Mechanics, vol. 70, no. 1, pp. 10–17, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 Y. Gong, B. Li, and Z. Li, “Immersedinterface finiteelement methods for elliptic interface problems with nonhomogeneous jump conditions,” SIAM Journal on Numerical Analysis, vol. 46, no. 1, pp. 472–495, 2008. View at: Publisher Site  Google Scholar  MathSciNet
 I. Babuška, “The finite element method for elliptic equations with discontinuous coefficients,” Computing, vol. 5, pp. 207–213, 1970. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 A. Hansbo and P. Hansbo, “An unfitted finite element method, based on Nitsche's method, for elliptic interface problems,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 4748, pp. 5537–5552, 2002. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 P. Colella and H. Johansen, “A Cartesian grid embedded boundary method for Poissons equation on irregular domains,” Journal of Computational Physics, vol. 60, pp. 85–147, 1998. View at: Google Scholar
 M. Oevermann and R. Klein, “A Cartesian grid finite volume method for elliptic equations with variable coefficients and embedded interfaces,” Journal of Computational Physics, vol. 219, no. 2, pp. 749–769, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 M. Oevermann, C. Scharfenberg, and R. Klein, “A sharp interface finite volume method for elliptic equations on Cartesian grids,” Journal of Computational Physics, vol. 228, no. 14, pp. 5184–5206, 2009. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 W. Ying and C. S. Henriquez, “A kernelfree boundary integral method for elliptic boundary value problems,” Journal of Computational Physics, vol. 227, no. 2, pp. 1046–1074, 2007. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 S. Hou, W. Wang, and L. Wang, “Numerical method for solving matrix coefficient elliptic equation with sharpedged interfaces,” Journal of Computational Physics, vol. 229, no. 19, pp. 7162–7179, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 S. Hou, Z. Li, L. Wang, and W. Wang, “A numerical method for solving elasticity equations with interfaces,” Communications in Computational Physics, vol. 12, no. 2, pp. 595–612, 2012. View at: Publisher Site  Google Scholar  MathSciNet
 S. Hou, L. Wang, and W. Wang, “A numerical method for solving the elliptic interface problems with multidomains and triple junction points,” Journal of Computational Mathematics, vol. 30, no. 5, pp. 504–516, 2012. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 L. Wang, S. Hou, and L. Shi, “A numerical method for solving 3D elasticity equations with sharpedged interfaces,” International Journal of Partial Differential Equations, vol. 2013, Article ID 476873, 13 pages, 2013. View at: Publisher Site  Google Scholar
 L. Wang, S. Hou, and L. Shi, “An improved nontraditional finite element formulation for solving the elliptic interface problems,” Journal of Computational Mathematics. In press. View at: Google Scholar
 S. Hou, P. Song, L. Wang, and H. Zhao, “A weak formulation for solving elliptic interface problems without body fitted grid,” Journal of Computational Physics, vol. 249, pp. 80–95, 2013. View at: Publisher Site  Google Scholar
 L. Wang, A numerical method for solving the elliptic and elasticity problems [Ph.D. Dissertation], Louisiana Tech University, 2011.
 S. Hou and X.D. Liu, “A numerical method for solving variable coefficient elliptic equation with interfaces,” Journal of Computational Physics, vol. 202, no. 2, pp. 411–445, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
Copyright
Copyright © 2013 Liqun Wang and Liwei Shi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.