Elliptic Interface Problems in Heterogeneous Media

Updated 28 December 2025

Elliptic interface problems are partial differential equations with discontinuous coefficients across internal interfaces, essential for modeling composite materials and multiphysics phenomena.
Specialized numerical strategies such as immersed finite elements, weak Galerkin, and discontinuous Galerkin methods accurately enforce jump conditions on non-body-fitted meshes.
Recent research incorporates machine learning and adaptive algorithms to tackle complex geometries and improve convergence and simulation accuracy.

Elliptic interface problems are partial differential equations (PDEs) with piecewise coefficients and interfacial jump conditions, modeling highly heterogeneous media or multiphysics phenomena. The defining feature is the presence of an internal interface Γ across which the solution and/or its flux may be discontinuous due to abrupt changes in material properties. Accurate numerical treatment of such problems is central to numerous applications, including composite materials, subsurface flows, electrostatics with dielectric interfaces, and biological modeling. Recent advances address both algorithmic and theoretical challenges associated with discretizing solutions on meshes that do not align with the geometric interface, combining techniques spanning finite element, finite volume, discontinuous Galerkin, and machine learning paradigms.

1. Mathematical Formulation and Fundamental Properties

Let Ω ⊂ ℝ^d (d = 2, 3) be a bounded domain split by a smooth interface Γ into subdomains Ω⁻ and Ω⁺. The canonical elliptic interface problem seeks $u: \Omega \to \mathbb{R}$ such that

(PDE in each subdomain)

$-\nabla \cdot (\beta(x) \nabla u(x)) = f(x) \quad \text{in } \Omega^- \cup \Omega^+,$

with piecewise-constant $\beta(x)$ .

(Interface jump conditions on $\Gamma$ )

$[u] := u^+ - u^- = g_D, \qquad [\beta \partial_n u] := \beta^+ \partial_n u^+ - \beta^- \partial_n u^- = g_N,$

with $n$ denoting the normal vector to Γ.

(Boundary condition)

$u = g \quad \text{on } \partial\Omega.$

The interface discontinuities in either material coefficients or source terms mandate the enforcement of both value and flux jumps at Γ; these can be homogeneous or inhomogeneous, depending on physical context. Weak formulations are cast in broken Sobolev spaces (e.g., $H^1(\Omega^- \cup \Omega^+)$ ), typically involving integration by parts over Ω⁻ and Ω⁺ individually and explicit inclusion of jump terms in variational identities. Well-posedness, norm equivalence, coercivity, and a priori estimates have been systematically established for both primal and mixed (first-order) formulations (Mu et al., 2018, Yang et al., 13 Sep 2024, Cao et al., 2022, Hu et al., 30 Jan 2024).

2. Immersed and Unfitted Finite Element Approaches

A distinguishing challenge is the treatment of interfaces that do not match mesh lines. Immersed finite element (IFE) methods construct basis functions that incorporate the jump conditions at the element level, allowing for interface-independent (unfitted) grids. On interface elements, polynomial basis functions are made piecewise, with their coefficients determined by Lagrange interpolation, strong or weak jump enforcement, and flux continuity (Guo et al., 2016, Guo et al., 2019, Lin et al., 2015).

IFE construction: For triangular or quadrilateral elements, local IFE shape functions are constructed to satisfy nodal interpolation conditions and interfacial value and flux continuity—sometimes at a selected point (Sherman-Morrison rank-one modification) or in a weak (integral) sense (Guo et al., 2016). Systematic constructions extend to 3D trilinear, bilinear, and high-order polynomial spaces (Guo et al., 2019, Adjerid et al., 2023).
Approximation properties: When geometric and regularity conditions are met (e.g., maximum interface curvature times mesh size is small, suitable interface intersection rules), optimal-order convergence is attained: $O(h)$ in energy norm, $O(h^2)$ in $L^2$ norm, with constants explicit in curvature and coefficient contrast (Guo et al., 2016, Adjerid et al., 2023, Yang et al., 13 Sep 2024, Mu et al., 2018).
Mixed and conservative variants: Augmentation with local piecewise constants achieves local mass conservation (EIFEM) (Jo et al., 2021), while mixed FE formulations (e.g., hybrid or BDM-based) yield $H(\mathrm{div})$ conforming fluxes and robust inf-sup stability, even for unfitted meshes (Cao et al., 2022, Hu et al., 30 Jan 2024).
Stability enhancements: Penalty-based stabilization, particularly at interface edges (PPIFE, SIPFEM, HPFEM), is critical to controlling the jump in IFE functions and restoring coercivity under mesh refinement (Lin et al., 2015, Cao et al., 2015, Wu et al., 2010).
High-order and DG extensions: Geometry-conforming approaches leveraging local Frenet frames allow arbitrary-degree polynomial IFE construction, with exact interface mapping and enforcement of all necessary jump conditions, yielding robust $p$ -convergence (Adjerid et al., 2023).

3. Weak Galerkin, Hybrid, and Unfitted Discretizations

Weak Galerkin (WG) and hybrid continuous/discontinuous approaches further decouple interface enforcement from mesh geometry.

Immersed Weak Galerkin (IWG): The solution is approximated by pairs $(v_0, v_b)$ , with $v_0$ in local IFE or polynomial space and $v_b$ (its trace) on element boundaries, possibly discontinuous. Weak derivatives (gradients and divergences) are defined through local integration by parts. Stabilization via penalties on $(v_0 - v_b)$ secures well-posedness independent of interface location (Mu et al., 2018, Yang et al., 13 Sep 2024).
Hybrid CG/WG/IFE: Away from the interface, standard continuous elements are used; on interface elements, immersed WG or IFE spaces are employed. This design achieves optimal $H^1$ and $L^2$ convergence with minimal overhead, demonstrated in complex geometries and coefficient regimes (Yang et al., 13 Sep 2024).
DG and unfitted penalty methods: Discontinuous Galerkin formulations (e.g., cut-FEM, unfitted IP-DG) enforce jump conditions weakly via Nitsche-type penalty and consistency terms. Local mesh merging strategies handle stability and conditioning in the presence of small interface cuts, ensuring the inf-sup constant and condition number are independent of interface-mesh alignment (Han et al., 2023, Capatina et al., 4 Jul 2025, Wu et al., 2010).

4. Finite Difference and Volume Methods for Complex Interfaces

High-fidelity finite difference and finite volume methods are formulated to maintain accuracy and mass conservation on structured grids with complex or moving interfaces.

Matched Interface and Boundary (MIB) FVM: The MIB-FVM approach uses Cartesian vertex-centered control volumes with fictitious values assigned at irregular stencils, determined to enforce physical jump conditions to second order. This hybrid maintains sharp enforcement of interface physics and exact global conservation, achieving $O(h^2)$ accuracy even for low-regularity solutions and intricate geometries (Cao et al., 2015).
Higher Order Compact Immersed Methods: Schemes such as HEJIIM utilize nine-point compact stencils and explicit Taylor-corrected jump series, achieving $O(h^4)$ convergence throughout the domain, including near non-smooth interfaces with singular sources (Singhal et al., 2021). Compact coupling techniques like CCIM combine Smereka's discrete delta regularization and coupling equations to recover second-order solution and gradient accuracy on arbitrarily complex domains (Zhang et al., 2021).

5. Machine Learning and Mesh-Free Methods

Recent advances leverage neural network methodologies to represent solution discontinuities and interface effects without explicit meshing.

Discontinuity Capturing Shallow Neural Networks: DCSNN represents the piecewise solution as a continuous function in $(d+1)$ -dimensional space, using an additional label for subdomain. A shallow (single-layer) mesh-free network is trained on the PDE, boundary, and interface jump residuals, achieving superior accuracy and efficiency, particularly for high-dimensional problems suffering from the curse of dimensionality in grid-based methods (Hu et al., 2021).
Cusp-Capturing Physics-Informed Neural Networks (PINNs): By augmenting the network input with a cusp-enforced feature (e.g., $|\phi(x)|$ for level-set $\phi$ ), the surrogate can sharply represent solution kinks at interfaces where flux jumps but the solution remains continuous. This approach delivers mesh-free, high-accuracy solutions for non-smooth interfaces in 1D–6D settings (Tseng et al., 2022).
Domain Decomposition Learning Algorithms: The Dirichlet–Neumann learning framework treats subdomain solutions via separate deep networks, enforces interface continuity via boundary penalties, and iterates using a mesh-free domain decomposition strategy, analyzed rigorously for penalty parameter scaling with coefficient contrast (Sun et al., 2023).

6. A Posteriori Error Control, Adaptivity, and Advanced Analysis

State-of-the-art methods address error estimation, adaptivity, and geometric/interface complexity.

Conservative Flux Recovery and Adaptive CutFEM: Using hybrid mixed formulations with local Lagrange multipliers and immersed Raviart-Thomas flux recovery, robust a posteriori error estimators are constructed that are reliable and locally efficient, driving adaptive mesh refinement toward interface and singularity regions (Capatina et al., 4 Jul 2025).
Boundary Correction for Curved Interfaces: For body-fitted meshes where geometric interface approximation mismatches the physical interface, Taylor-based boundary correction methods transfer Neumann (flux) and Dirichlet data from the curved interface to the mesh, restoring optimal convergence without requiring complex integration on curved elements (Hou et al., 5 Feb 2025).
Kernel-Free Boundary Integral Methods (KFBI): Interface problems on manifolds or embedded surfaces are treated by recasting singular boundary integrals as equivalent interface PDEs on Cartesian grids. Fast multigrid solvers and interface-corrected finite differences yield $O(h^2)$ global accuracy with $O(N^2)$ computational cost (Yin et al., 22 Aug 2025).

7. Practical Performance and Implementation Considerations

Uniform theoretical convergence rates (typically $O(h)$ in energy norm and $O(h^2)$ in $L^2$ ) are universally demonstrated across IFE, WG, DG, FV, and neural network methods, for a wide range of geometries, coefficient contrasts, and discontinuous data. Robustness to interface location and orientation is achieved by proper stabilization, local geometric constraints (e.g., maximum-angle selection), and penalty parameter tuning. Efficient solvers—including auxiliary space preconditioning, AMG for enriched spaces, and geometric multigrid for interface-corrected FD/FV schemes—ensure scalability and mesh/contrast independent convergence (Jo et al., 2021, Mu et al., 2018, Han et al., 2023, Capatina et al., 4 Jul 2025). Extensions and current research address adaptivity, moving and topologically complex interfaces, high-order geometry representation, and hybridization with neural operator surrogates.