Optimal strategies to control human behaviors during a catastrophic event

1.   Introduction

In practical applications, partial differential equations (PDEs) can characterize and describe various physical phenomena and engineering problems effectively. In many engineering fields, elliptic PDEs with discontinuous coefficients are an important class of mathematical and physical equations, which can be used to describe various physical phenomena and engineering problems, such as bi-material interface cracks, multi-phase flows, fluid-structure interaction, as well as crystal growth. These situations involve multiple materials or fluids with distinct properties on the entire physical domain. These physical problems and natural phenomena modeled by various types of elliptic PDEs, referred to as "interface problems", may exhibit some non-smooth characteristics, such as discontinuities, high gradients, and singularities. The non-smooth features on the interfaces will significantly affect the regularity of the solution for model problem ^[1,2]. As a result, these non-smooth features on the interface pose difficulties in accurately representing and simulating the interactions between different materials or fluids.

Due to their flexibility in mesh generation and their ability to handle problems with irregular geometries, finite element methods (FEMs) are now viewed as powerful tools for solving PDEs that occur in various engineering disciplines. When addressing interface problems in a finite element framework, the discrete background mesh can be categorized as fitted mesh and unfitted mesh, based on the alignment between the mesh and the interface. Classical FEMs often exhibit poor convergence rates when solving interface problems under the unfitted mesh framework, while satisfactory accuracy can be achieved if fitted meshes are used, as shown in the left figure of Figure 1. For a triangulation ${\mathcal T}_h = \{K\}$ of a bounded domain $\Omega \subset {\mathbb R}^2$, the chunkiness parameter $c_K$, defined as the ratio of the circumradius ${h_K}$ of triangle $K$ to the inradius ${\rho_K}$ of that triangle, measures the minimum angle of triangle $K$ and is a significant indicator of mesh quality. A large chunkiness parameter indicates that the triangle $K$ is extremely elongated, which leads to poor mesh quality. The right figure of Figure 1 illustrates the mesh quality of the fitted mesh mentioned earlier and clearly shows yellow and orange triangles with larger chunkiness parameters. A triangulation ${\mathcal T}_h = \{K\}$ based on a bounded domain $\Omega \subset {\mathbb R}^2$ is called shape-regular if $c_K = \frac{h_K}{\rho_K}\leq \beta$ for every element $K\in{\mathcal T}_h$, where the positive constant $\beta$ does not depend on the given mesh parameter ^[3]. It is extremely difficult to construct high-quality meshes that are tailored to match different types of interfaces in cases where the computational domain or interface involve complex shapes. Moreover, a significant amount of time and computational resources would be required to create an interface-fitted mesh with high quality, especially when dealing with time dependent interface problems. For example, when calculating the solution to problems involving a moving interface, it is necessary to perform mesh partitioning at different temporal layers since the numerical solution at the previous temporal layer may be projected onto the new mesh at the current temporal layer ^[4].

Although FEMs based on a unfitted mesh suffer from low precision and a poor rate of convergence in various norms when handling interface problems, unfitted mesh methods based on fixed quasi-uniform meshes have always attracted extensive attention in the field of numerical solutions for PDEs due to their advantages, such as flexibility, adaptability to dynamic physical processes, and handling of irregular domains ^[5]. Therefore, numerical methods with unfitted meshes, such as the immersed FEM ^[6,7,8,9], the immersed finite volume element method (FVEM) ^{[10,11,12,13,14]}, and the generalized or extended FEM (GFEM/XFEM) ^{[15,16,17,18,19,20]}, have made significant advancements for solving interface problems since the beginning of this century. Given the fundamental similarity in the ideas of GFEM and XFEM, they will be hereafter collectively referred to as GFEM.

In the context of immersed FEM and immersed FVEM, the basis functions close to the interface require some local modifications such that the given jump conditions associated with the interface problem are preserved as accurately as possible, while standard basis functions are employed on elements located away from the interface, as usual. Hence these two methods maintain the same number of basis functions as their classic counterparts. Both conforming and non-conforming immersed FEMs have been presented in ^[6,7,8]. It was shown in ^[6] that the conforming immersed FEM yields satisfactory accuracy and convergence in terms of the $L^{\infty}$ norm, while the non-conforming immersed FEM fails to do so in some cases of significantly large ratios of jump coefficients. Additionally, numerical results in ^[9] have already shown that the classic immersed FEM exhibits comparatively larger point-wise errors around the interface and leads to poor convergence behavior in both $L^2$ and $L^{\infty}$ norms when the mesh is more refined. To address these issues, the authors of ^[9] developed a partial-penalty immersed FEM, in which the penalty strategies were adopted to maintain the desired convergence rates in the $L^2$ and $L^{\infty}$ norms for elliptic interface problems. In order to locally preserve the physical conservation laws, the authors of ^[10] presented an immersed FVEM to solve elliptic interface problems in the late 1990s; however, the numerical results sometimes exhibited oscillating features when there was a large ratio of jump coefficients. In ^[11,12,13], various techniques such as the source removal technique and adding stability terms were employed to improve the traditional immersed FVEM for solving elliptic interface problems. At the same time, the convergence of approximate solutions in the sense of the $L^{\infty}$ norm was investigated through extensive numerical examples. In ^[14], a new immersed FVEM without stability terms was proposed and the authors explored the convergence of the $L^{\infty}$ norm through numerical examples.

The GFEM is an extended version of the traditional FEM, incorporating non-polynomial shape functions into the standard FEM approximation space. However, the GFEM may cause severe ill-conditioning of the linear equations ^[18,21,22]. Recently, a stable GFEM has grown to become a very efficient and robust numerical method for addressing the issue of ill-conditioning in GFEM ^{[21,22,23,24]}, as it satisfies these excellent properties: (a) It can achieve an optimal convergence rate (first-order accuracy for linear finite element shape functions) in the sense of the energy norm; (b) the growing speed of the scaled stiffness matrix conditioning exhibits a trend almost identical to that of traditional FEM; (c) it exhibits robustness as the relative distance between the interface and mesh boundaries decreases. In ^[21,22,23], the SGFEM was utilized to solve elliptic interface problems, and a detailed convergence analysis of the energy norm was provided in ^[23]. In ^[24], a higher-order SGFEM was developed, and it has been used to solve elliptic interface problems ^[25] and elliptic eigenvalue problems ^[26], where it achieved optimal higher-order convergence rates. In ^[27,28], a strongly SGFEM with singular and distance functions as enrichment functions was introduced to solve specific non-smooth interface problems with corner interfaces. To reduce computational complexity, an improved SGFEM was investigated for interface problems in ^[29], which did not require the evaluation of distance function gradients. Additionally, a proof of convergence analysis of the energy norm was provided. The SGFEM has also been extensively applied to other mathematical models and physical problems with interfaces, such as crack problems ^[30,31], parabolic interface problems ^[4,32], fracture mechanics ^[33,34,35], and so on.

In the SGFEM studies mentioned above, the study of convergence mostly focuses on the energy norm or the $H^1$ semi-norm; and what is more, most of the interfaces associated with model problems are trivial in the numerical experiments. However, as far as we know, there has been limited studies on the numerical investigation of $L^{\infty}$ norm convergence using the SGFEM for the elliptic interface problems with complicated interfaces. In this paper, a modified SGFEM is presented to solve elliptic interface problems with complicated but smooth interfaces, where a one-side distance function or level set function serves as the enrichment function. The modified SGFEM is based on an unfitted mesh, and there are no stability terms or penalty parameters in the weak form of the model problem, that is to say, it is a conforming FEM. The numerical convergence behavior of the modified SGFEM is verified with respect to the $L^2$ and $L^{\infty}$ norms through several non-trivial numerical examples, and at the same time, its stability and robustness are evaluated.

The remaining sections are structured as follows. The elliptic interface problem and its corresponding weak form are presented in Section 2. In Section 3, the methodologies of standard SGFEM and proposed SGFEM with a one-side enrichment function are explained in detail. In Section 4, several non-trivial numerical examples are presented, in which the relevant interfaces have different shapes to illustrate the convergence, stability, and robustness of the presented SGFEM. The conclusions are drawn in the last section.

2.   Elliptic interface problem

In this paper, we will consider a classical second-order elliptic interface problem defined on a bounded and convex domain $\Omega \subseteq {\mathbb R}^2$. The model problem with a Neumann boundary condition is stated as follows:

where the solution domain $\Omega$ is divided into the two mutually disjoint subdomains $\Omega_0$ and $\Omega_1$ by a smooth interface $\Gamma = \overline{\Omega}_0\cap\overline{\Omega}_1$, as shown in Figure 2, and $\vec{n}_b$ is the unit outward normal vector of the boundary $\partial{\Omega}$. The coefficient function $a(\mathit{\boldsymbol{x}}) = a_i$ for $\mathit{\boldsymbol{x}} \in {\Omega}_i, i = 0, 1$, in which $a_0$ and $a_1$ are known positive constants, is discontinuous on the interface $\Gamma$. The jump conditions of the solution and its gradient are specified as

where $\vec{n}$ stands for the unit outward normal vector of the given smooth interface $\Gamma$, and the notation $\left[u \right]_\Gamma = u|_{{\Omega}_0}-u|_{{\Omega}_1}$ represents the difference of the function $u$ restricted to the given smooth interface $\Gamma$.

The weak form of the classical second-order elliptic interface problem (2.1) is described as follows: Find $u\in \mathcal{E}(\Omega)$, such that

where the bilinear form $B(\cdot, \cdot)$ and linear functional $F(\cdot)$ are defined as

and the energy space is expressed as

The compatibility condition $F(1) = 0$ holds true for the given functions $f$, $g$, and $q$. For the elliptic interface problem (2.1) with Neumann boundary condition, the weak solution $u$ is uniquely determined in the energy space $\mathcal{E}(\Omega)$, up to a constant term. The imposition of an essential constraint on the boundary $\partial \Omega$ allows the determination of a unique solution.

Remark 2.1. The focus of the discussion now shifts to the clarification of the interface jump conditions (2.3). In many physics problems, it is common to come across weakly or strongly discontinuous interface problems ^[19,20]. Weak discontinuities are expressed through kinks in the solution function $u$, where the function is continuous on the interface (i.e., $\left[u \right]_\Gamma = 0$), and a discontinuity is present in the gradient (the jump in the flux may be non-homogeneous, i.e., $\left[a \frac{\partial u}{\partial \vec{n}} \right]_\Gamma = q \neq 0$). On the other hand, strong discontinuities refer to jumps in the solution function, where the solution function $u$ satisfies $\left[u \right]_\Gamma \neq 0$ on the interface. In this paper, we only consider the case of weak discontinuity.

3.   The modified SGFEM

This section begins with an intensive review of GFEM, followed by an introduction to the stable GFEM (SGFEM) with a one-side distance function.

For a bounded and convex domain $\Omega$, let ${\mathcal T}_h = \{K\}$ be a regular and uniform finite element mesh with size $h$. The mesh consists of closed quadrilateral elements $K$, such that $\Omega = \cup_{K \in {\mathcal T}_h} K$. The nodes related to the finite element mesh ${\mathcal T}_h = \{K\}$ are represented as $\{\mathit{\boldsymbol{x}}_i\}_{i \in I_h}$, where $I_h$ is the index set. For each node $\mathit{\boldsymbol{x}}_i$, the classical bilinear finite element basis function is defined as $\phi_i$, and its support is represented as $supp\{{\phi}_i\} = \bar{\omega}_i$, where ${\omega}_i$ is the patch related to $\mathit{\boldsymbol{x}}_i$. These finite element basis functions $\{\phi_i\}_{i\in I_h}$ possess the property of partition of unity, i.e., $\sum\limits_{i\in I_h} \phi_i(\mathit{\boldsymbol{x}}) = 1$.

For the weak form (2.4), the approximation solution $u_h$ for the standard FEM is defined as follows:

where $\Bbb{S}_{FEM}$ is the standard FEM approximation space, and $u_i$ is the nodal degrees of freedom (DOFs) associated with the basis function $\phi_i$. It is widely acknowledged that the standard FEM can obtain satisfactory numerical results only in the case of smooth problems, whereas standard FEMs with unfitted meshes have poor approximation properties when addressing interface problems ^[22]. To improve approximation accuracy, some special (non-polynomial) basis functions are added into the initial FEM space to establish a good approximation space with the capability of mimicking the non-smooth features on the interface. Consequently, the standard FEM approximation space $\Bbb{S}_{FEM}$ is converted into the GFEM approximation space $\Bbb{S}_{GFEM}$. The approximate solution ${u}_h$ for the interface problem, which belongs to the approximation space $\Bbb{S}_{GFEM}$, can be written as

where the set $I_{h, enr} \subset I_h$ represents the index set of enriched nodes; the function $\Pi_i$, referred to as enrichment function, usually depends on the specific problem; the coefficients $c_i$ and $b_i$ are standard and enriched nodal DOFs related to the components $\phi_i$ and $\phi_i\Pi_i$, respectively; and $\Bbb{S}_{ENR}$ stands for the global enrichment space.

For smooth interface problems with weak discontinuities, a so-called distance function

or the absolute value of the level set function

which possesses a weak discontinuous property is commonly used as the enrichment function ^[18,22,36]. It can be seen that ${D}(\mathit{\boldsymbol{x}})$ is equal to zero for any point $\mathit{\boldsymbol{x}}\in\Gamma$, and the continuity is maintained in ${D}(\mathit{\boldsymbol{x}})$ but not in its derivative. Unfortunately, the system matrix could be badly conditioned and bring about a large condition number when the GFEM is applied to the interface problem, leading to the GFEM approximation solution (3.2) having poor accuracy for interface problems ^{[21,22,23,24,32]}. To improve the accuracy and conditioning of the GFEM, a stable GFEM with a local modification of the enrichment function was developed in ^[21] and further studied in ^[22,23,24]. Therefore, the approximation solution ${u}_h$, belonging to the SGFEM approximation space $\Bbb{S}_{SGFEM}$ can be corrected as

and

where $\mathcal{I}_h$ is an interpolation operator with the finite element hat function as the interpolation function. In ^[32], an SGFEM with a one-side enrichment function $\tilde{D}$ (see (3.7)) was proposed to solve parabolic interface problems, where the corresponding approximation function was written as

where the one-side enrichment function is given by

We note that the modified enrichment functions $D - \mathcal{I}_h D$ and $\tilde{D} - \mathcal{I}_h \tilde{D}$ are equal to zero at the nodes and not equal to zero in the interior of the interface elements. The SGFEM approximation solution obtained from (3.5) or (3.6) satisfies these excellent properties: (a) It can achieve an optimal convergence rate (first-order accuracy for linear finite element shape functions) in the sense of the energy norm; (b) the growth speed of the scaled stiffness matrix conditioning exhibits a trend almost identical to that of traditional FEM; (c) it remains robust as the relative distance between the interface and mesh boundaries decreases. Certainly, it is possible to replace the one-side enrichment function with

For convenience, in this paper, SGFEMs with a one-side enrichment function $(3.7)$ (or $(3.8)$) is referred to as $\mbox{SGFEM0}$ (or $\mbox{SGFEM1}$). The enrichment functions are problem-dependent, and different enrichment strategies are used to approximate different types of interface problems. For interface problems with strong discontinuities, it can be observed that there is a jump in the solution on the interface. As a result, the Heaviside step or sign function is commonly chosen as an enrichment function to model the strong discontinuity of the model problem on the interface ^[19,20]. For the non-smooth interface problem with a corner on the interface, its solution always has singularity at the corner, thus the modified singular function and modified distance function are required to approximate the solution more accurately and to ensure convergence ^[27,28].

Using the Ritz-Galerkin method, the variational problem (2.4) can be discretized in the SGFEM approximation space $\Bbb{S}_{SGFEM}$ spanned by the finite element basis functions $\phi_i$ and non-polynomial basis functions $\phi_i(\tilde{D} - \mathcal{I}_h \tilde{D})$. Let $\Bbb{S}_{SGFEM}\subset \mathcal{E}(\Omega)$ be a finite dimensional space, and the SGFEM solution $u_h\in\Bbb{S}_{SGFEM}$ satisfies

In light of the Galerkin orthogonality, we have the following error estimate:

which implies that the approximate solution $u_h$ belonging to the finite dimensional space $\Bbb{S}_{SGFEM}$ is the best approximation in the energy space ${\mathcal E}(\Omega)$. Hence, a space with good approximation properties is crucial in ensuring that the SGFEM produces accurate approximations for non-smooth problems. One noteworthy feature of the suggested SGFEM is that the bilinear form $B(\cdot, \cdot)$ does not include additional stabilization terms in the discrete weak form (3.9) and the computational expense of assembling the stiffness matrix and load term is only marginally higher than that of conventional FEM. Furthermore, the stiffness matrix derived from the discrete weak form (3.9) remains unchanged regardless of variations in the gradient jump conditions.

The stiffness matrix ${\mathbf A}$ associated with the aforementioned weak form (3.9) is a symmetric definite matrix. The condition number of the stiffness matrix ${\mathbf A}$ is a pivotal indicator when solving linear equation systems. Its scaled condition number (SCN) is now discussed in detail. The SCN of matrix ${\mathbf A}$ is measured by the condition number of scaled matrix ${\mathbf{\hat{A}}} = {\mathbf D}{\mathbf A}{\mathbf D}$, where ${\mathbf D}_{ii} = {\mathbf A}_{ii}^{-1/2}$ is a diagonal matrix. Then, the SCN $\mathcal{K}$ of stiffness matrix ${\mathbf A}$ is given by

where $\kappa_2({\mathbf {\hat{A}}})$ refers to the spectral condition number of the invertible matrix ${\mathbf {\hat{A}}}$.

Remark 3.1. We mention that the GFEM with the enriched node set $I_{h, enr}^{\Gamma}$ is called topological GFEM, and the GFEM with the enriched nodes within a specified region adjacent to the interface is referred to as geometric GFEM ^[22,32]. These two methods have some shortcomings in solving interface problems. Topological GFEM cannot obtain the first-order accuracy in the energy norm for linear elements, while geometric GFEM brings more enriched DOFs and leads to bad SCN of stiffness matrix ^[22]. It is obvious that the SGFEM and topological GFEM have a minimum of enriched DOFs when compared with the geometric GFEM, and their extra computational costs are minimal among these GFEMs (including SGFEM). Therefore, in the following numerical experiments, topological and geometric GFEMs are not taken into consideration.

4.   Numerical experiments

In the experiment section, we consider the standard FEM along with various SGFEMs. These SGFEMs include SGFEM, SGFEM0, and SGFEM1. The convergence, stability, and robustness of the proposed SGFEM are verified through various types of numerical examples. In the numerical experiments, we discretize the bounded domain $\Omega$ using an $N \times N$ rectangular mesh. Assuming $a_0 = 1$ or $a_1 = 1$, several numerical tests with different contrasts $\rho$ (defined by $\rho = \frac{a_1}{a_0}$) are conducted. In certain specific cases, numerical tests about variable coefficients are also implemented.

For comparison, the numerical relative errors of the approximation solution $u_h$ in the sense of $L^2$ and $L^{\infty}$ norms are given by

where

denote the $L^2$ and infinity norms of the function $u$ over the whole domain $\Omega$. To calculate $\|u\|_{L^{\infty}(K)}$, we choose $10\times10$ uniformly distributed points on the element $K$ and compute the maximum absolute value of $u$ among these sampling points. Similarly, we compare the SCN of the stiffness matrix ${\mathbf A}$ measured by (3.11).

4.1.   Convergence and SCN

In this subsection, the convergence and stability of various SGFEMs are verified by some numerical examples with straight or curved interfaces. For a straight or curved interface given by the set $\Gamma = \{\mathit{\boldsymbol{x}} = (x, y) \in \Omega : \gamma(\mathit{\boldsymbol{x}}) = 0 \}$, the domain $\Omega \cap \{\mathit{\boldsymbol{x}}: \gamma(\mathit{\boldsymbol{x}}) < 0\}$ is denoted as $\Omega_0$ and $\Omega \cap \{\mathit{\boldsymbol{x}}: \gamma(\mathit{\boldsymbol{x}}) > 0\}$ is denoted as $\Omega_1$.

Example 4.1. In this example, a straight interface problem defined on the solution domain $\Omega = (0, 1)^2$ is considered, where the interface $\Gamma$ is defined by a straight line segment $\{\mathit{\boldsymbol{x}} \in \Omega : y-\tan(\theta_0)(x-1-d_0)-1 = 0 \}$, and its analytical solution $u$ is represented by

with $\alpha_0 = 2$, $\theta_0 = \frac{\pi}{6}$, and $d_0 = \frac{1}{\pi}$, where the point $A(1+d_0, 1)$ outside the solution domain is the center of the polar coordinate system $(r, \theta)$ with the polar axis $\{\mathit{\boldsymbol{x}}: x > 1+d_0, y = 1\}$. This numerical example was introduced in ^[37].

The mesh size is set as $h = \frac{1}{2^{k+1}}$, where $k = 1, 2, \cdots, 7$, in this example. The relative errors associated with various SGFEMs and FEM are presented in Figure 3 for different contrasts $\rho$. It is noted that the relative errors of the approximation solution $u_h$ in both $L^2$ and $L^{\infty}$ norms are slightly larger for $\rho = 100$ compared to $\rho = 10$. As might be expected, it is clear form Figure 3 that the standard FEM only has poor convergence rates in both $L^2$ and $L^{\infty}$ norms, since the mesh is not an interface-fitted mesh. It can be observed that various SGFEMs (including SGFEM, SGFEM0, and SGFEM1) exhibit second-order accuracy in the sense of $L^2$ and $L^{\infty}$ norms. The SCNs ${\mathcal K}_{FEM}$, ${\mathcal K}_{SGFEM}$, ${\mathcal K}_{SGFEM0}$, and ${\mathcal K}_{SGFEM1}$ associated with the stiffness matrices are plotted against $h$ in Figure 4, where it can be observed that they grow at the same rate, although the magnitude of the SCNs obtained by various SGFEMs is slightly larger compared to those of the standard FEM. These numerical results show that the SGFEMs possess excellent convergence and are indeed stable.

Example 4.2. In this example, a circular interface problem with different flux conditions is considered, where the circular interface is defined as $\Gamma = \{\mathit{\boldsymbol{x}} \in \Omega : (x-x_0)^2+ (y-y_0)^2 -r_0^2 = 0 \}$ centered at $O(x_0, y_0)$ with radius $r_0$. The classical solution $u$ of the interface problem, with a homogeneous jump condition in flux $\left[a \frac{\partial u}{\partial \vec{n}} \right]_\Gamma = 0$, is chosen as

and the classical solution $u$ of the interface problem, with a non-homogeneous jump condition in flux $\left[a \frac{\partial u}{\partial \vec{n}} \right]_\Gamma = (a_1-a_0)\sin(xy)(y\cos\theta+x\sin\theta) \neq 0$, is given by

where the point $O(x_0, y_0)$ inside the given domain $\Omega = (0, 1)^2$ is the center of the polar coordinate system $(r, \theta)$, and the polar axis is $\{\mathit{\boldsymbol{x}}: x > x_0, y = y_0\}$. For the numerical example, we will take $x_0 = \frac{1}{\sqrt{5}}, y_0 = \frac{1}{\sqrt{3}}$, and $r_0 = \frac{1}{\sqrt{10}}$. This numerical example with the homogeneous jump condition in flux was introduced in ^[37].

The mesh size is set as $h = \frac{1}{2^{k+1}}, k = 1, 2, \cdots, 7$. The relative errors with homogeneous and non-homogeneous jump conditions in flux, associated with the different SGFEMs and FEM are presented in Figures 5–7 for different contrasts $\rho$. It is evident that all SGFEMs yield second-order accuracy in terms of various norms, whereas the standard FEM only achieves first-order accuracy. The approximate solution obtained by various SGFEMs in the $L^2$ norm is slightly larger than the second-order accuracy for both cases of smaller and larger contrasts ($\rho = 1/1000$ and $\rho = 1000$). The SCNs of the stiffness matrices obtained by different SGFEMs and FEM are plotted with respect to $h$ in Figure 8 for different contrasts $\rho$. The figure shows that various SGFEMs (including SGFEM, SGFEM0, and SGFEM1) are well conditioned, that is to say, their SCNs grow at the rate of $O(h^{-2})$ with further mesh refinement for different contrasts $\rho$. Thus, the proposed SGFEM is indeed stable.

Example 4.3. This example focuses on the interface problem with variable coefficients, with the computational domain and interface matching those in Example 4.2. The variable coefficients are assumed to be $a_0(\mathit{\boldsymbol{x}}) = \frac{x^2+y^2+7}{7}$ and $a_1(\mathit{\boldsymbol{x}}) = \frac{xy+50}{5}$. The classical solution $u$ of the interface problem, with a non-homogeneous jump condition in flux $\left[a \frac{\partial u}{\partial \vec{n}} \right]_\Gamma = \left(a_0(\mathit{\boldsymbol{x}})-a_1(\mathit{\boldsymbol{x}})\right)\left(\frac{\partial a_1}{\partial x}\cos\theta+\frac{\partial a_1}{\partial y}\sin\theta\right)r_0^5 \neq 0$, is given by

The relative errors and SCNs against mesh size $h = \frac{1}{2^{k+1}}, k = 1, 2, \cdots, 7$, are reported in Figures 9 and 10, respectively. It is observed that various SGFEMs attain the optimal convergence order $O(h^{2})$ in various norms, while the FEM only attains the suboptimal convergence order $O(h)$. Figure 10 shows that the SCNs ${\mathcal K}_{A}$ of various SGFEMs grow at the expected growth rate $O(h^{-2})$.

In the following, three interface problems with more complicated interfaces are considered, where the curvature of the interfaces is not constant.

Example 4.4. This numerical example was presented in ^[38]. The interface is a curved interface $\Gamma = \{\mathit{\boldsymbol{x}} \in \Omega : y-3x(x-0.3)(x-0.8)-0.38 = 0 \}$ which is defined on the solution domain $\Omega = (-1, 1)^2$, and the analytical solution $u$ is chosen as

The relative errors and the SCNs against mesh size $h = \frac{2}{2^{k+1}+1}, k = 1, 2, \cdots, 7$, are shown in Figures 11 and 12, respectively. It can be seen from Figure 11 that various SGFEMs achieve the optimal convergence order $O(h^{2})$ in both $L^2$ and $L^{\infty}$ norms, while the FEM only achieves the suboptimal convergence order $O(h)$ in both $L^2$ and $L^{\infty}$ norms. Figure 12 shows that the SCNs ${\mathcal K}_{A}$ of different SGFEMs grow at the expected rate $O(h^{-2})$ with further mesh refinement for different contrasts $\rho$, which indicates that they are stable.

Example 4.5. In this case, the test example has an ellipsoid interface $\Gamma = \{\mathit{\boldsymbol{x}} \in \Omega : \frac{x^2}{\alpha^2}+\frac{y^2}{\beta^2}-1 = 0 \}$ with a semimajor axis $\alpha = \frac{\pi}{4.71}$ and a semiminor axis $\beta = \frac{\pi}{6.28}$, as illustrated in Figure 13. The solution $u$ defined on the computational domain $\Omega = (-1, 1)^2$ is given as follows:

The relative errors and SCNs are plotted against the mesh size $h = \frac{2}{2^{k+1}+1}, k = 1, 2, \cdots, 7$, in Figures 14 and 15. As can be seen from Figure 14, the approximate solutions obtained by different SGFEMs yield second-order accuracy in terms of various norms, while the FEM has poor convergence. It is evident from Figure 15 that the SCNs ${\mathcal K}_{A}$ obtained by different SGFEMs grow at the almost identical rate $O(h^{-2})$ for different contrasts $\rho$.

Example 4.6. Next, a more complicated interface problem from ^[8] is investigated. The problem is defined on the domain $\Omega = (-1, 1)^2$ and the interface is a flower-like curve defined as $\Gamma = \{\mathit{\boldsymbol{x}} \in \Omega : r^4\big(1+0.4\sin(6\theta)\big)-0.3 = 0 \}$, as shown in Figure 13. Its analytical solution $u$ is given by

where the point $(0, 0)$ is the center of the polar coordinate system $(r, \theta)$, and the polar axis is the right half of the x-axis.

The relative errors and SCNs are plotted against the mesh size $h = \frac{2}{2^{k+1}+1}, k = 1, 2, \cdots, 7$, and the similar numerical behaviors are observed in Figures 16 and 17. Examination of the figures shows that the relative errors measured by different SGFEMs diminish approximately quadratically in the sense of $L^2$ and $L^{\infty}$ norms. Figure 17 clearly shows that the SCNs of various SGFEMs grow at a rate of $O(h^{-2})$ with further mesh refinement for different contrasts $\rho$.

4.2.   Robustness test

To study the robustness of the suggested SGFEM, a mesh of $8\times8$ elements is discretized on the unit square domain $\Omega = [0, 1]^2$. Additionally, the interface $\Gamma$ is defined by a horizontal line segment $\{\mathit{\boldsymbol{x}} \in \Omega : 0\leq x\leq1, y = \alpha+\delta\}$ with the discontinuous coefficient

where $\alpha = 0.25\, \, \, \mbox{or}\, \, \, \alpha = 0.50, \delta_j = 0.06\times 2^{-j+1}, j = 1, 2, \cdots, 20$. It is evident that some mesh elements are cut by the interface. As the parameter $\delta$ decreases, the horizontal straight line segment $\Gamma$ approaches the edges of the fixed mesh. Figure 18 exhibits the changes in the SCN of the stiffness matrices for different SGFEMs and standard FEM. It is clear that the SCNs of SGFEMs are independent of the relative distance $\frac{\delta}{h}$ in both interface scenarios, which indicates all SGFEMs are robust as the relative distance between the interface and mesh boundaries decreases.

5.   Conclusions

In this study, a modified SGFEM was proposed to resolve the second-order elliptic interface problem involving non-trivial interfaces. The proposed SGFEM utilizes a one-side enrichment function, which is a straightforward extension of the SGFEM with a two-side distance function. Several non-trivial numerical examples on SGFEM with a one-side enrichment function were investigated, showing that the modified SGFEM can achieve the optimal convergence order $O(h^{2})$ in both $L^2$ and $L^{\infty}$ norms for linear elements. The numerical experiments also show that the proposed SGFEM is stable and robust, and applicable to any smooth interface, regardless of its concavity or shape. A study on the higher-order SGFEM for parabolic interface problems will be investigated in our forthcoming paper.

Author contributions

Pengfei Zhu: Conceptualization, Methodology, Investigation, Software, Data curation, Validation, Funding acquisition, Project administration, Supervision, Writing—original draft. Kai Liu: Conceptualization, Methodology, Software, Writing—review and editing. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work was supported by the Growth Project for Young Scientific and Technological Talents of Ordinary Colleges and Universities in Guizhou Province (No. KY[2022]317), the Basic Research Project for Natural Science of Guizhou Province (No. ZK[2022]218), and the Science and Technology Planning Project of Guizhou Province (No. [2023]159).

Conflict of interest

The authors declare that they have no relevant competing interests.