1 Introduction↩︎

The adaptive finite element method (AFEM) is one of the effective numerical approaches for solving partial differential equations, particularly showing strong adaptability for problems involving singularities or multiscale features. The fundamental idea behind adaptive finite element methods is to achieve an optimal asymptotic error decay rate by refining the finite element space. Generally, there are three primary approaches to refining a finite element space: the \(h\)-adaptive method, which improves the accuracy of finite element approximation by refining the mesh; the \(p\)-adaptive method, which uses a fixed mesh but increases the polynomial degree of the shape functions to improve its accuracy; and \(r\)-adaptivity (or moving mesh methods), which relocates mesh nodes to concentrate them in regions of interest without changing the total number of nodes. Naturally, these strategies can be synergistically combined to achieve greater efficiency and robustness than any single approach alone. Among them, the two most prominent combinations are the \(hp\)- and \(hr\)-adaptive methods, coupling \(h\)-adaptivity with \(p\)- and \(r\)-adaptivity, respectively. The \(hp\)-adaptive method is renowned for achieving an exponential convergence rate in the energy norm [1], [2]. However, its practical implementation, especially in high dimensions, is notoriously challenging due to the complexity of managing both mesh refinement and variable polynomial degrees. Given these challenges, this work focuses on the \(hr\)-adaptive framework.

Irrespective of the specific adaptive strategy, the key mechanism of an AFEM is to drive the mesh adaptation process (refinement or coarsen) based on an a posteriori error indicator. The a posteriori error estimator, which is a computable quantity depending on the numerical solution, provides information about the distribution of the error of the finite element approximation [3], [4]. In particular, for the \(hr\)-AFEM, the success of such an algorithm hinges on a reliable a posteriori error estimator to guide both mesh refinement and node distribution to achieve better accuracy with minimum degrees of freedom. A posteriori error analysis, along with convergence and optimality theories for adaptive finite element algorithms, are well established for elliptic equations [5]–[7]. Corresponding theories are available for parabolic equations [8]–[10]. However, translating this theoretical framework into a robust computational procedure is nontrivial, as applying AFEM to time-dependent problems presents unique challenges not encountered in stationary elliptic settings. The sequential nature of time-stepping and the evolution of solution features exacerbate the computational cost of mesh adaptation and the error propagation from frequent data transfer. Consequently, conventional AFEMs suffer from the following inherent fundamental bottlenecks when applied to such problems:

• Inefficient refinement: The refinement process in each time layer can require up to \(30-40\) iterations [10], [11].

• Interpolation challenge: The interpolation of the solution from the previous time step onto each newly refined mesh is algorithmically demanding and a significant computational bottleneck.

• Redundant initial mesh choice at each time step: At each new time step, the adaptive cycle begins with the final mesh from the previous step, resulting in a linear system with an excessively large number of degrees of freedom. Here, an additional and intricate mesh coarsening step (COARSEN) must be enforced.

In pursuit of an efficient and straightforward adaptive finite element method for parabolic problems, this work proposes a novel neural network-enhanced \(hr\)-adaptive finite element algorithm. Its key innovation lies in employing a data-driven neural network as our primary tool, which serves as a continuous function representation. The choice of this tool is natural, given the proven capacity of neural networks to act as universal function approximators and to learn efficient representations from data. Nowadays, neural networks have achieved remarkable success in artificial intelligence, enabling complex tasks in domains such as speech recognition, natural language processing, computer vision, and autonomous vehicles through data-trained models [12], [13]. In scientific computing, neural networks leverage their dual strengths—as universal approximators [14]–[18] and for mitigating high-dimensional complexity [19]—underpin their fundamental capability to handle complex, data- and physics-driven computational tasks. This capability is vividly evidenced by the growing number of studies that directly apply neural networks to solving partial differential equations [20]–[27].

Despite this rapid expansion, several key limitations have become apparent when deep learning is used as a stand-alone PDE solver, including high training costs, limited interpretability, and accuracy that often remains below that of well-established numerical methods. To overcome these limitations while leveraging the strengths of both artificial intelligence and classical numerical methods, a growing trend is to employ artificial intelligence as a tool for augmentation rather than replacement, shifting the role of neural networks from stand-alone PDE solvers toward specialized components embedded within mature numerical frameworks. In this vein, neural networks have been incorporated into finite element methods in various ways. Examples include error-informed adaptive FEM enhanced by neural networks [28] and \(h\)-adaptive finite element–interpolated neural networks designed to resolve localized features [29]. Further related developments and applications can be found in [30]–[33] and the references therein. These works illustrate how neural networks can be strategically integrated to complement, rather than replace, classical numerical methods.

Following the strategy of neural networks as specialized modules, our proposed method effectively addresses the three aforementioned limitations. This method is more efficient in terms of adaptive performance compared to the conventional refinement strategies discussed earlier. Specifically, the advantages of the proposed neural network-enhanced \(hr\)-AFEM in addressing the three limitations of the conventional AFEMs are as follows:

• Efficient refinement: By controlling the refinement of mesh nodes at each step, the number of iterations required for convergence is reduced to a maximum of \(7\) iterations.

• Interpolation free: By training a neural network to replace the prior finite element solution, the method bypasses the complex and costly steps of interpolation. With only a single fine-tuning required per time step, it then enables direct function evaluation at any spatial point.

• Cheap initial mesh choices per time step: By employing a combined error indicator that incorporates information from two adjacent time steps, the method ensures reliable re-initializing of the mesh to a prescribed, coarse uniform state at each time step.

Together, these mechanisms—drastic iteration reduction, elimination of interpolation costs, and avoidance of redundant initial computations—constitute a comprehensive efficiency gain. The rest of the paper is outlined as follows: Section 2 introduces the model problem and the fully discrete scheme for the parabolic equation. Then, we recap the conventional adaptive algorithm, offering a comprehensive analysis of its methodology while critically examining its key limitations and potential drawbacks. In response, Section 3 discusses the development of an efficient neural network-enhanced \(hr\)-adaptive finite element algorithm, specifically highlighting the main advantages over conventional methods. In Section 4, several numerical experiments are presented to demonstrate the robustness and effectiveness of the proposed method.

2 Standard adaptive finite element algorithm↩︎

2.1 Model problem↩︎

Let \(\Omega\) be a polygonal or polyhedral domain in \(\mathbb{R}^d\) (with \(d=2,3\)), and \(T>0\). Consider the following parabolic equation: \[\label{model} \left\{\begin{align} u_t - \nabla\cdot(a\nabla u) &= f, & \text{in } & \Omega \times (0, T], \\ u(\mathbf{x}, t) &= g, & \text{on } & \partial\Omega \times (0, T], \\ u(\mathbf{x}, 0) &= u_0, & \text{in } & \Omega, \end{align}\right.\tag{1}\] where \(f\in L^2(0, T; L^2(\Omega))\), \(u_0\in L^2(\Omega)\), and \(a\in L^{\infty}(\Omega)\) are uniformly bounded and positive in the sense that there exist positive constants \(\lambda_{\min}\) and \(\lambda_{\max}\) satisfying \[\lambda_{\min}\leq a(\mathbf{x}) \leq\lambda_{\max},\quad \forall \mathbf{x}\in\bar{\Omega}.\] This type of equation plays a critical role in applications such as heat conduction models [34], image denoising [35], material science [36], and financial mathematics [37]. The solution to the parabolic equation is generally smooth in most regions. However, discontinuities in initial conditions, nonuniform material properties, or interactions between multiple physical processes can give rise to singularities that emerge and gradually evolve over time. This requirement arises from the fact that the smoothness of solutions to parabolic equations may vary considerably across different regions, necessitating finer spatial discretization and more sophisticated numerical techniques to resolve the solution behavior accurately.

Throughout this work, we use the standard notation \(W^{m,q}(D)\) to refer to Sobolev spaces on an open set \(D\subset \mathbb{R}^d\), with the norm \(\Vert\cdot\Vert_{W^{m,q}(D)}\) and the semi-norm \(\vert\cdot\vert_{W^{m,q}(D)}\). The space \(W^{m,q}_0(D)\) is defined as \(\{w \in W^{m,q}(D): w|_{\partial D} = 0\}\). Furthermore, \(W^{m,2}(D)\) and \(W^{m,2}_0(\Omega)\) are denoted by \(H^m(D)\) and \(H^m_0(D)\), respectively. We also denote by \(L^{p}(0,T;H^k(D))\) the space of all \(L^p\) integrable functions from \((0,T)\) into \(H^k(D)\) with norm \(\|v\|_{L^{p}(0,T;H^k(D))}=\left(\int_0^T \|v\|_{k,D}^pdt\right)^{\frac{1}{p}}\) for \(p\in [1,\infty)\).

Given \(f\in L^2(0, T; L^2(\Omega))\) and \(u_0\in L^2(\Omega)\), the variational formulation of 1 reads: Find \(u \in L^2(0, T; H^1(\Omega))\bigcap H^{1}(0, T; H^{-1}(\Omega))\) such that \[\label{weak32form} \left\{\begin{align} (u_t,v)+(a\nabla u,\nabla v) & =(f,v), & \forall v \in H_0^1(\Omega),\quad t\in(0,T], \\ u(\mathbf{x},0) & =u_0. \end{align}\right.\tag{2}\]

2.2 Fully discrete FEM↩︎

Let the time step size be defined as \(\tau = T/S\), which partitions the interval \([0, T]\) uniformly into \(S\) subintervals and produces the discrete time levels \(t_n = n\tau\), for \(n = 0, \dots, S\). Let \(\mathcal{T}_h^n\) be a shape-regular triangulation of \(\Omega\) at the time level \(t_n\), whose associated linear finite element space is \[V_h^n=\{v\in H_0^1(\Omega):v|_K\in P_1(K),\forall K\in \mathcal{T}_h^n\},\] where \(P_1(K)\) denotes the space of linear polynomials on \(K\subseteq R^d(d=2,3)\). Let \(u_h^0=Pu_0\), where \(P: L^2(\Omega)\rightarrow V_h^0\) is the \(L^2\) projection operator onto the finite element space \(V_h^0\) defined on the initial mesh \(\mathcal{T}_h^0\). Applying the backward Euler method for time discretization and the finite element method for spatial discretization to the weak formulation 2 yields the following fully discrete scheme: Find a sequence of functions \(u_h^n \in V_h^n,\,n=1,2,\cdots, S\), such that \[\label{full-discrete} \left\{\begin{align} \left(\frac{u_h^n-\Pi_h^n u_h^{n-1}}{\tau},v\right)+(\nabla{u_h^n},\nabla v ) &= \left(f_h^{n},v\right), \quad & \forall v \in V_h^n, \\ u_h(\mathbf{x},0) & =u_h^0, \end{align}\right.\tag{3}\] where the operator \(\Pi_h^n: V_h^{n-1} \to V_h^n\) denotes the nodal interpolation operator from the previous time level’s finite element space onto the current one. Note that at each time level \(t_n\), the adaptive finite element algorithm generates a sequence of meshes \({\mathcal{T}^{n,0}_h, \mathcal{T}^{n,1}_h, \dots}\) along with their corresponding continuous piecewise linear finite element spaces \({V^{n,0}_h, V^{n,1}_h, \dots}\). Hereafter, we do not distinguish between the intermediate meshes generated within the \(n\)-th time step and denote them collectively by \(\mathcal{T}_h^n\), and similarly by \(V_h^n\) for the associated finite element spaces.

2.3 Recap of the standard \(h\)-AFEM↩︎

The standard \(h\)-adaptive finite element method for parabolic equations proceeds sequentially in time, with multiple mesh updates possible within a single time step. Given the numerical solution \(u_h^{n-1}\) from the previous time level \(t_{n-1}\), the iterative algorithm for computing the solution \(u_h^n\) at time \(t_n\) can be summarized as the following workflow:

Here, the INTERPOLATE step computes \(\Pi_h^n u_h^{n-1} \in V_h^n\) from \(u_h^{n-1}\), and is required whenever a mesh change occurs. The essential COARSEN step is applied once, at the end of the adaptive iteration loop within a time step. Its role is to coalesce elements with low estimated error. Typically, these elements lie in regions that were highly refined to capture features of the previous solution \(u_h^{n-1}\), but become redundant for subsequent computations once the new solution \(u_h^n\) is obtained. If not, these elements would accumulate, contributing to an unchecked growth in the total number of degrees of freedom and ultimately leading to computational blow-up. The process of refinement and coarsening will be further illustrated in 1. While mesh coarsening improves computational efficiency, it risks coarsening elements that were just refined in the previous iteration, which can degrade the accuracy of the subsequent solution.

We also recap the gradient recovery-based a posteriori error estimators, which use a certain norm of the difference between direct and post-processed gradient approximations as an indicator, have been widely used in the engineering and scientific computation community since the work of Zienkiewicz and Zhu [38]. The success of recovery-based error estimators stems from their computational efficiency, ease of implementation, and general asymptotic exactness [9], [39]–[42]. Let \(G: V_h^n\rightarrow V_h^n\times V_h^n\) denote the gradient recovery operator, and let \(\mathcal{N}_h^n\) denote the set of vertices of mesh \(\mathcal{T}_h^n\). For each node \(z\in\mathcal{N}_h^n\), let \(\omega_z\) be the element patch associated with \(z\), we define the recovered gradient \(G(\nabla u_h^n(z))\) as \[G(\nabla u_h^n(z)):=\sum_{K_i\in \omega_z}\frac{1/|K_i|}{\sum_{K_j\in \omega_z}1/|K_j|}\nabla u_h^n(c_i),\] where \(c_i\) denotes the center of the element \(K_i\). The recovered gradient \(G(\nabla u_h^n)\) over the whole domain is then obtained via interpolation: \[G(\nabla u_h^n)=\sum_{z\in\mathcal{N}_h^n}G(\nabla u_h^n(z))\phi_z^n,\] where \(\phi_z^n\) is the Lagrange basis of finite element space \(V_h^n\) associated with \(\mathcal{T}_h^n\). The corresponding local and global gradient recovery-based error estimators are defined separately as \[\label{estimators} \hat{\eta}_{h, K}^n:=\|G(\nabla u_h^n)-\nabla u_h^n\|_{0, K},\quad \hat{\eta}_{h}^n:=\|G(\nabla u_h^n)-\nabla u_h^n\|.\tag{4}\]

By integrating the ideas outlined above, we introduce 2.

2.4 Limitations of \(h\)-AFEM↩︎

It is important to highlight three primary limitations when applying 2 in practical computations. First, the bisection algorithm is a layer-by-layer refinement method, which can lead to excessive adaptive iterations. Coupled with the fact that each iteration executes the full INTERPOLATE-SOLVE–ESTIMATE–REFINE cycle, the overall computational expense becomes prohibitive. A canonical example is the Laplace equation \(-\Delta u = 0\) on the L‑shaped domain \([-1, 1]^2 \setminus (0,1) \times (-1,0)\) with Dirichlet boundary conditions, whose solution is \[u(x,y)=r^{2/3}\sin(2\theta/3),\quad r=\sqrt{x^{2}+y^{2}},\quad \theta=\tan^{-1}(y/x).\] Applying the standard adaptive finite element method with newest‑vertex bisection (see [11]) to this problem yields the meshes shown in 3. The minimal element size decreases from \(h_{K,\text{min}}=2^{-1}\sqrt{2}\) in the initial mesh to \(h_{K,\text{min}}=2^{-5}\sqrt{2}\) after 10 adaptive steps, yet achieving a quasi‑optimal convergence rate typically requires \(30‑40\) such refinements. This exemplifies the high iterative cost inherent to the classical approach. More numerical tests and analysis can be found in [11].

Second, the interpolation of the solution from the previous mesh onto the newly refined mesh poses a significant implementation challenge. This operation requires solving a point‑location problem by determining, for every newly created node, which element of the previous mesh contains it. A straightforward implementation based on exhaustive element search is prohibitively expensive. Although spatial search structures can accelerate this search, their reliable and robust integration remains algorithmically demanding. Crucially, this geometrically non‑trivial step must be executed after every mesh refinement, making it a recurrent computational bottleneck.

Third, the standard practice of beginning each new time step’s adaptive process directly from the final mesh of the previous step incurs substantial and persistent overhead. This forces all key computational components—matrix assembly, linear solver, error estimation, and data interpolation—to handle this inflated problem size in every iteration. Moreover, to avoid carrying forward this irreversible computational baggage, the standard method is forced to insert an additional and intricate mesh coarsening step. This step specifically targets and removes the fine mesh elements that were created to capture features of the old solution \(u_h^{n-1}\). We show the performance of 2 in 1.

3 Neural network-enhanced \(hr\)-AFEM↩︎

Having identified the causes of these limitations, this section introduces a solution framework composed of two essential tools and a central strategy. Specifically, we employ: (i) a mesh generation tool to perform local refinement; (ii) a data transfer tool to handle interpolation between non-matching meshes; and (iii) the mesh-reset and estimator-combine strategy to ensure a more economical number of degrees of freedom in each adaptive cycle. Together, these components give rise to our main proposal: a novel neural network-enhanced \(hr\)-AFEM.

3.1 Mesh size field↩︎

The first objective is to reduce the total number of adaptive iterations. In [11], an efficient adaptive finite element method was proposed for elliptic equations by employing a refinement strategy that can control the number of mesh vertices, ensuring that the adaptive process terminates within seven iterations. In this section, we extend that strategy to the adaptive finite element method for the parabolic equation 1 . To materialize this adaptive refinement, a robust and controllable mesh generation tool is required. Several established libraries such as Gmsh, MeshLab, TetGen, and NETGEN provide control, either directly or indirectly, to meet these requirements. Among them, Gmsh [43] is a widely adopted and well-documented option. Its key advantage is that refinement is driven by a user-supplied mesh size field, which determines the local mesh density. Given its widespread acceptance and precise controllability, we employ Gmsh as the primary mesh generation tool in this work.

The mesh size field is a scalar function defined over the computational domain that prescribes the target local mesh size. It serves as the direct link between the error estimator and the resulting mesh geometry: a smaller prescribed size in a region leads to a denser mesh there, and vice versa. In our framework, the size field is designed to strategically concentrate mesh density in regions with high error contribution by identifying nodes where the error estimator exceeds a threshold. Prior to introducing the algorithm, we present the necessary notation for the Hadamard product and its related operations as follows.

Then we provide 5 to generate such a mesh size field for Gmsh.

3.2 Neural networks↩︎

The meshes generated by Gmsh with the user-specified size field are unstructured and non-nested, making the implementation of the interpolation \(\Pi_h^nu_h^{n-1}\) both algorithmically complex and computationally expensive. As the second key tool, we employ a neural network-based mesh‑free surrogate to construct \(u_h^{n-1}\). Notably, unlike the mesh‑dependent interpolation \(\Pi_h^n u_h^{n-1}\) which must be recomputed after each mesh change, the neural network surrogate is mesh‑free and trained only once per time step.

3.2.1 General neural networks↩︎

Consider a general \(L\)-layer neural network \(u_{\theta}:\mathbb{R}^{d_0} \rightarrow \mathbb{R}^{d_L}\). Formally, the \(k\)-th layer is defined as a function \(\mathbf{N}^{(k)}(x)\) whose output is given by applying an element-wise activation function \(\sigma\) to an affine transformation of the input. Denoting the weight matrix and bias vector of this layer by \(\mathbf{W}_k \in \mathbb{R}^{d_k \times d_{k-1}}\) and \(\mathbf{b}_k \in \mathbb{R}^{d_k}\), respectively, the forward propagation is written as \[\mathbf{N}^{(k)}(\mathbf{x}^{(k-1)}):=\sigma(\mathbf{W}_k\mathbf{x}^{(k-1)}+\mathbf{b}_k).\] The activation function \(\sigma\) is crucial, as it provides the necessary nonlinearity that enables the network to learn complex patterns and approximate a wide range of functions. In the absence of such nonlinearity, the network would collapse into a sequence of linear transformations, severely restricting its ability to capture intricate patterns and generalize to unseen data. For regression, no activation is used in the output layer to allow for unconstrained output. Therefore, the general neural network is given by \[\label{neuralnetwork0} u_{\theta}(\mathbf{x})=\mathbf{W}_L\left(\mathbf{N}^{(L-1)}\circ \cdots \circ \mathbf{N}^{(2)}\circ \mathbf{N}^{(1)}(\mathbf{x})\right)+\mathbf{b}_L.\tag{5}\]

While numerous activation functions, such as ReLU, sigmoid, and tanh, have been proposed in the literature, we employ the hyperbolic tangent function in this work, \[\sigma(x)=\tanh{x} = \frac{e^x-e^{-x}}{e^x+e^{-x}}.\] This choice is motivated by its smoothness, bounded range, and symmetry about zero, which are well suited for approximating the smooth solution profiles commonly encountered in partial differential equations.

In general, the training of the neural network proceeds as follows. Let \(\theta=((\mathbf{W}_0, \mathbf{b}_0), \cdots, (\mathbf{W}_L, \mathbf{b}_L))\) represent all trainable parameters of the network. Given a function \(u\in C^0(\Omega)\), we train a neural network \(u_{\theta}(\mathbf{x})\) that approximates \(u\) using the dataset \(\{(\mathbf{x}_i, u(\mathbf{x}_i)\}_{i=1}^N\). The network parameters \(\theta\) are optimized by minimizing the empirical loss function \(\mathcal{L}(u_{\theta})\), typically the mean squared error (MSE): \[\theta^*=arg \min\limits_{\theta\in \Theta}\mathcal{L}(u_{\theta}),\quad \mathcal{L}(u_{\theta}):=\frac{1}{N}\sum\limits_{i=1}^N(u_{\theta}(\mathbf{x}_i)-u(\mathbf{x}_i))^2,\] where \(\Theta\) is the feasible set consisting of the parameters and hyperparameters of the neural network.

3.2.2 Mesh-free surrogate via neural networks↩︎

The objective here is to construct a neural network surrogate \(u_{\theta}^{n-1}\), based on the general neural network 5 , to approximate the finite element solution \(u_h^{n-1}\) obtained at \(t_{n-1}\). This surrogate provides a continuous, mesh-free representation of \(u_h^{n-1}\) for use in the subsequent adaptive procedure at time \(t_n\), thereby eliminating the need for the computationally expensive interpolation \(\Pi_h^n u_h^{n-1}\) in 3 . The surrogate \(u_{\theta}^{n-1}\) is trained using an input–output dataset constructed from the vertices of \(\mathcal{T}_h^{n-1}\) and the corresponding finite element approximations, namely \[\{ (\mathbf{x}, u_h^{n-1}(\mathbf{x})) \mid \mathbf{x} \in \mathcal{N}_h^{n-1} \},\] where \(\mathcal{N}_h^{n-1}\) denotes the set of vertices of the mesh \(\mathcal{T}_h^{n-1}\).

Motivated by [44], the Dirichlet boundary condition \(u|_{\partial \Omega} = g(\mathbf{x},t), \;\mathbf{x} \in \partial \Omega\) is imposed by constructing the neural network surrogate \(u_{\theta}^{n-1}\) as \[\label{neuralnetwork} u_{\theta}^{n-1}(\mathbf{x}) = d(\mathbf{x}) u_{\theta}(\mathbf{x}) + \tilde{g}(\mathbf{x},t_{n-1}),\tag{6}\] where \(u_{\theta}(\mathbf{x})\) denotes the general neural network defined in 5 . The function \(d(\mathbf{x})\) is assumed to be continuous and satisfies \[\label{dxcond} d(\mathbf{x})=0 \quad \text{for } \mathbf{x}\in\partial\Omega, \qquad d(\mathbf{x})\neq 0 \quad \text{for } \mathbf{x}\in\Omega,\tag{7}\] A natural candidate is to adopt the exact distance function, \[\label{dxformula} d(\mathbf{x}) =\min\limits_{\mathbf{y}\in\partial\Omega}|\mathbf{x}-\mathbf{y}|.\tag{8}\] Nevertheless, encoding the exact distance function 8 in neural networks entails considerable computational overhead [44]. In literature, an approximate distance function \(d(\mathbf{x})\) is usually employed. In our application, \(d(\mathbf{x})\) is not restricted to the approximate distance function. Any continuous function satisfying 7 is sufficient to enforce the Dirichlet boundary condition exactly through the network construction. The specific choice of \(d(\mathbf{x})\) utilized in the numerical examples is provided in 5. Furthermore, \(\tilde{g}(\mathbf{x},t_{n-1}) \in H^1(\Omega) \cap C^0(\Omega)\) is an extension of the boundary data \(g(\mathbf{x},t_{n-1})\) from \(\partial\Omega\) to the entire domain \(\Omega\), satisfying \(\tilde{g}|_{\partial \Omega} = g\).

For a homogeneous Dirichlet boundary condition \(u|_{\partial \Omega} = 0, \;\mathbf{x} \in \partial \Omega\), the neural network surrogate \(u_{\theta}^{n-1}\) in 6 reduces to \[\label{neuralnetwork43} u_{\theta}^{n-1}(\mathbf{x}) = d(\mathbf{x}) u_{\theta}(\mathbf{x}).\tag{9}\] As a simple illustration, consider a two-point boundary value problem on the interval \([a, b]\) with nonhomogeneous Dirichlet boundary conditions \(u(a, t_{n-1})=A\) and \(u(b, t_{n-1})=B\). In this case, one may choose \(\tilde{g}(x, t_{n-1}) = \frac{b-x}{b-a}A + \frac{x-a}{b-a}B\), which is simply the linear interpolant of the prescribed boundary data.

Then, the empirical loss function is defined as MSE: \[\label{lossFunction} \mathcal{L}\!\left(u_{\theta}^{n-1}; \mathcal{T}_h^{n-1}\right) = \frac{1}{\lvert \mathcal{T}_h^{n-1} \rvert} \sum_{\mathbf{x} \in \mathcal{N}_h^{n-1}} \left| u_{\theta}^{n-1}(\mathbf{x}) - u_h^{n-1}(\mathbf{x}) \right|^2,\tag{10}\] where \(\lvert \mathcal{T}_h^{n-1} \rvert\) denotes the number of vertices in \(\mathcal{T}_h^{n-1}\).

The optimal surrogate \(u_{\theta^*}^{n-1}(\mathbf{x})= d(\mathbf{x}) u_{\theta^{n-1}_*}(\mathbf{x}) + \tilde{g}(\mathbf{x},t_{n-1})\) is obtained by solving \[\theta_*^{n-1} = \arg\min_{\theta} \mathcal{L}\!\left(u_{\theta}^{n-1}; \mathcal{T}_h^{n-1}\right).\] The resulting minimizer \(u_{\theta^*}^{n-1}(\mathbf{x})\) provides an accurate approximation of \(u_h^{n-1}(\mathbf{x})\) and therefore serves as a mesh-free surrogate for \(u_h^{n-1}\) over the entire computational domain.

By integrating the optimized neural network surrogate \(u_{\theta^*}^{n-1}\) into the fully discrete finite element scheme 3 , we obtain the following reformulated variational problem.

The stability estimate ?? shows that the scheme is \(L^2\)-stable for any \(\tau>0\). Consequently, if the neural network surrogate contains only a small perturbation, then this perturbation does not amplify through the time-stepping procedure, and the resulting solution remains stable. Although the time step size \(\tau\) does not affect the unconditional stability of ?? , it has a significant impact on both the solution accuracy and the quality of the neural network initialization \(u_{\theta^*}^{n-1}\). A detailed numerical investigation is presented in 3.

3.3 \(hr\)-adaptive finite element algorithm↩︎

To effectively solve the neural network–enhanced finite element scheme ?? , we require an adaptive algorithm equipped with a reliable a posteriori error estimator to guide mesh refinement. Note that the accuracy of the term \((u_{\theta^*}^{n-1},v)\) in ?? may also affect the accuracy of the numerical solution \(u_h^n\). Consequently, the mesh \(\mathcal{T}_h^n\) at time \(t_n\) must be carefully generated to capture the key features of both the previous approximation \(u_{\theta^*}^{n-1}\) and the current solution \(u_h^n\), since the locations and intensities of singularities may vary from one time step to the next. Therefore, the error estimator should account for contributions from both \(u_h^n\) and \(u_{\theta^*}^{n-1}\), since both play a role in determining an appropriate mesh. In light of this, we employ a combined error indicator in 4 as \[\label{etarec} \eta_{h, K}^n=\max\{\hat{\eta}_{h, K}^n,\hat{\eta}_{h, K}^{n-1}\}=\frac{\hat{\eta}_{h, K}^n+\hat{\eta}_{h, K}^{n-1}+|\hat{\eta}_{h, K}^n-\hat{\eta}_{h, K}^{n-1}|}{2}, \qquad \eta_h^n=\left(\sum\limits_{K\in\mathcal{T}_h^n}\eta_{h, K}^2\right)^{1/2},\tag{11}\] where \(\hat{\eta}_{h, K}^n\) is defined in 4 . The estimator 11 incorporates information from both time levels and can guide the new mesh to capture the singularities at \(t_n\) while ensuring accurate computation of the term \((u_{\theta^*}^{n-1},v)\).

In contrast to the coarsening procedure commonly used in conventional AFEM, the final mesh at \(t_{n-1}\) is not inherited in our algorithm. Instead, the adaptation process employs a mesh‑reset mechanism, whereby the mesh is re‑initialized from a predefined coarse mesh \(\mathcal{T}_h^{n,0}\) at each time level \(t_n\). As a result, the number of degrees of freedom remains controlled throughout the simulation, leading to reduced computational cost.

At each time level, the adaptive procedure proceeds systematically. Starting from the initial mesh \(\mathcal{T}_h^{n,0}\), five adaptive iterations are performed. Based on the reliable data set \(\{(\eta_h^{n, i}, \lvert\mathcal{T}_h^{n, i}\rvert)\}_{i=2}^4\), a least-squares fitting is used to approximate the model between the error estimator and the NOV \[\label{fitModel} \eta_h^n \approx cN^{-p}.\tag{12}\] With the fitted parameters \((c, p)\), the desired NOV is predicted as \[N = \left\lceil\sqrt[p]{c/eTol}\right\rceil\] to ensure \(\eta_h^n \leq eTol\). The predicted \(N\), together with the local error estimators, is then used to compute the mesh size field via 5. Subsequently, a new mesh \(\mathcal{T}_h^{n,5}\) is generated using Gmsh according to the computed size field, and equation ?? is solved on this mesh. Finally, a seventh adaptive iteration is performed on \(\mathcal{T}_h^{n,6}\) as a safeguard to verify the prescribed tolerance.

Based on the discussion above, given the approximation \(u_h^{n-1}\in V_h^{n-1}\) at the time step \(t_{n-1}\) and the selected error estimators, the workflow of the new \(hr\)-adaptive algorithm for the parabolic equation 1 at the \(n\)-th time step can be summarized as follows:

From the perspective of the adaptive workflow, the proposed method introduces three main modifications compared to the conventional one. Here, LEARN represents the training of the neural network surrogate to transfer the solution between successive meshes; SOLVE procedure solves the neural network–enhanced finite element scheme ?? ; SIZE refers to the computation of mesh size field described in 3.1; and GENERATE stands for creating a new adapted mesh based on that field. Following this workflow, we present a neural network-enhanced \(hr\)-adaptive algorithm for the parabolic equation 1 in 7.

Compared with 2, the main changes in 7 are as follows. The generation of the initial mesh \(\mathcal{T}_h^0\) is now performed within the time-stepping loop. The adaptive process is limited to a maximum of seven iterations per time step. The SIZE and GENERATE steps replace the REFINE and INTERPOLATE steps from 2, although their operations are fundamentally different. Similarly, the LEARN step takes the place of COARSEN, but functions in a distinct manner.

4 Numerical Examples↩︎

This section presents numerical experiments to demonstrate the efficiency of the neural network-enhanced \(hr\)-adaptive 7 in solving the parabolic equation in both 2D and 3D. The main purpose of the tests centers on three primary test cases which feature distinct moving singularities—rotation, diffusion, and splitting—to validate the algorithm’s performance in tracking dynamic sharp gradients. These 2D cases are then extended to 3D to assess scalability in resolving such singularities and show performance beyond planar geometries. Finally, two numerical examples, namely the Allen–Cahn equation and a parabolic problem defined on a domain with reentrant corners, are provided to further illustrate the robustness and applicability of the proposed adaptive algorithm to nonlinear parabolic equations and domains with geometric singularities. All examples are implemented using the fealpy library (version 1.2.0) [46].

In the following examples, unless stated otherwise, the conforming linear finite element is used in all numerical experiments, with the time step size \(\tau=0.01\). The iteration tolerance \(eTol\) is employed to control adaptive termination, with specific values provided in each respective example. It primarily affects the computation of \(iteRO\) in the sixth adaptation and determines whether a seventh adaptation is necessary. For neural network training, weights are initialized via the Kaiming method with biases initialized to zero. The \(\tanh\) function serves as the activation function for the hidden layers, while the output layer uses a linear activation with no bias term (\(b_L=0\)). The L-BFGS optimizer is employed to minimize the loss function 10 to ensure convergence (stagnation). The neural network hyperparameters, including the number of layers, the number of neurons per layer, and the learning rate, are selected through empirical tuning rather than a systematic grid search. The distance function \(d(\mathbf{x})\) and the continuous extension \(\tilde{g}(\mathbf{x}, t_{n-1})\) of the boundary data used in the neural network 6 for all numerical examples are provided in 5.

4.1 2D Example↩︎

The neural network 6 used in this section employs an architecture \([2, 40, 40, 40, 1]\), comprising two input neurons, three hidden layers of 40 neurons each, and one output neuron. It is adopted as a demonstration case and has undergone no optimization or parameter tuning.

10 illustrates the evolution of the neural network training epochs over time for different choices of the time-step size \(\tau\), where a vertical scale is adopted to improve the visibility of variations in training epochs. For every tested value of \(\tau\), the first time level is the most computationally demanding, requiring more than \(2000\) training epochs. Once this initial network has been trained, the training effort is substantially reduced at subsequent time levels. When \(\tau\) is small, the epoch counts quickly stabilize and remain nearly uniform over time. As \(\tau\) increases, however, the epoch counts exhibit larger fluctuations, suggesting that larger temporal increments produce greater variations in the complexity of the learning task between successive time levels. The total number of training epochs (excluding the initial step) and the average number of training epochs per time step (excluding the initial step) are reported in 1. From these results, we observe that smaller \(\tau\) reduces the discrepancy between consecutive time levels, thereby decreasing the number of training epochs required at each step. However, since a smaller \(\tau\) results in more time levels overall, the total computational cost in training may increase. In summary, as \(\tau\) increases, the total number of training epochs initially decreases and then increases. Moreover, as the time step size \(\tau\) becomes larger, the advantage of reusing the neural network parameters from the previous step gradually diminishes. Nevertheless, parameter reuse remains beneficial, as it still reduces the number of training epochs compared with retraining from scratch at every time level.

To further examine the optimization behavior under the adopted setting, 11 reports the training loss for \(\tau=0.01\) at different time levels. The loss curves show that the MSE at each time level decreases to a value between \(10^{-6}\) and \(10^{-5}\). At the first time level, approximately \(2000\) training epochs are required to reach this accuracy, whereas fewer than \(100\) epochs are sufficient at all subsequent time levels.

Next, we revisit the experiment reported in 1 with the time-step size fixed at \(\tau=0.01\) and present the corresponding average training epochs under different network architectures and initialization strategies. Specifically, 2 summarizes the average number of training epochs (excluding the initial time step), computed over three independent runs, required for the L-BFGS optimizer to reach its default stopping criteria. The results indicate that both the network architecture and the initialization strategy have some influence on the training cost. Nevertheless, all tested configurations converge within a reasonable number of epochs, demonstrating the robustness of the proposed approach. Among the network architectures considered, deeper networks generally require fewer training epochs than shallower ones, particularly for the wider architectures. This can be attributed to the greater representational capacity of deeper networks, together with the warm-start initialization from the previous time step, which facilitates rapid convergence through minor parameter updates. In contrast, shallower networks may require larger parameter updates to achieve comparable accuracy, resulting in higher training cost. Regarding initialization strategies, no single method consistently achieves the lowest training cost across all architectures. Nevertheless, Uniform, Xavier, and Kaiming initializations often lead to slightly faster convergence, whereas Normal initialization is less frequently the best-performing choice. Overall, the differences among initialization methods remain moderate, especially for deeper networks, suggesting that the convergence behavior of the proposed method is relatively insensitive to the choice of initialization.

Additionally, 12 compares our neural network–enhanced \(hr\)-adaptive Algorithm 7 with the conventional \(h\)-adaptive method in [10]. The results show that our algorithm achieves a significant reduction in the degrees of freedom throughout the computation. Regarding computational cost, although the initial step is more time-consuming due to the extensive neural network training required, the subsequent time steps are considerably faster. As a result, the proposed method attains a substantial overall time advantage, as quantified in 3.

10pt

In 13, the top-left picture shows the initial uniform mesh \(\mathcal{T}_h^{0,0}\) (\(98\) nodes), which also serves as the initial mesh for adaptive iterations at each time step. The next five pictures show the meshes over five successive refinements, with the NOV being \(180\), \(340\), \(681\), \(1388\), and \(12818\), respectively. It can be observed that the NOV approximately doubles at each step, and the mesh generated by 7 at the fifth step is significantly larger than that of the previous mesh. The process terminates before the sixth refinement as the global error estimator is smaller than the tolerance \(eTol=0.01\), demonstrating that our fitted model 12 indeed plays its intended role.

14 illustrates the evolution of the numerical solutions and the corresponding mesh at different time steps. The mesh dynamically refines around sharp gradients or singularities, efficiently allocating resources where needed, and reveals a peak rotating around the coordinate origin. Notably, the grids shown here differ significantly from those in 8 due to the smaller time step used, which causes the positions of two distinct singularities at adjacent time steps to nearly overlap. The quasi-optimal convergence rate of the adaptive algorithm for the gradient error is demonstrated in the left panel of 9.

In this example, the settings are the same as 3, except for the tolerance \(eTol=0.05\). Training loss histories at different time levels are presented in 15. Similar to 3, approximately 2000 training epochs are required at the initial time level to achieve the MSE below \(10^{-5}\). In contrast, fewer than 100 epochs are sufficient at all subsequent time levels to attain comparable accuracy. 16 illustrates the adaptive refinement process of the initial mesh over six iterations. The NOV at each iteration are \(556\), \(1088\), \(2182\), \(4436\), \(9069\), \(18678\), and \(38659\), respectively. Unlike the 3, the NOV of the sixth mesh is approximately twice that of the fifth mesh. Moreover, the need for a seventh mesh indicates that the error estimator obtained on the sixth mesh is still larger than the tolerance. This suggests that the mesh generated by the previous refinement did not meet the desired tolerance, even though the tolerance was used when calculating the NOV for the sixth mesh. This implies the necessity of one extra step after the least square fitting, namely the case \(k=5\) in 7. Four different time steps are shown in 18, where the mesh and numerical solutions illustrate a scenario in which a ring-shaped crater gradually decreases its radius over time. The quasi-optimal convergence rate of the adaptive algorithm for the gradient error is demonstrated in the middle panel of 9. 17 compares the exact and finite element solutions at \(t=1\). The numerical solution accurately captures the localized ring-shaped profile. Moreover, the \(L^\infty\)-error is about \(6 \times 10^{-3}\), which implies that the finite element approximation achieves satisfactory accuracy.

Consistent with 3, the computational settings remain unchanged with the tolerance \(eTol=0.01\). The evolution of the training loss across different time levels is shown in 19. To reach the desired MSE, the initial time level requires more than 500 training epochs. In contrast, fewer than 20 epochs are sufficient at all subsequent time levels to achieve a comparable level of accuracy. 20 shows the adaptive process of the initial mesh over six iterations, with the NOV at each iteration being \(118\), \(206\), \(401\), \(813\), \(1677\), and \(3463\), respectively. 21 shows the evolution of the numerical solutions and the corresponding mesh at four different time steps. The refinement process dynamically captures the splitting of the central peak into two distinct moving peaks, maintaining a coarser resolution in smoother regions while ensuring high resolution in areas with sharp gradients. The right panel of 9 shows that the adaptive algorithm reaches the quasi-optimal convergence rate in terms of the gradient error.

4.2 3D Example↩︎

Extending to three dimensions, the neural network is updated to the architecture \([3, 32, 32, 32, 32, 1]\) with the meaning as given earlier. Owing to the lack of fully developed data structures for three-dimensional adaptive refinement in the fealpy framework, no direct three-dimensional comparisons are included in this study. Nevertheless, this circumstance serves to highlight a salient advantage of the proposed algorithm: its adaptive mechanism is inherently dimension-agnostic. Consequently, the entire implementation requires no modification when extending from 2D to 3D, underscoring the flexibility and robustness of our approach.

Apply 7 with \(eTol = 0.1\). The evolution of the training loss across different time levels is shown in 23. To reach the desired MSE, the initial time level requires more than 1000 training epochs. In contrast, fewer than 20 epochs are sufficient at all subsequent time levels to achieve a comparable level of accuracy. 24 illustrates the mesh refinement process at initial time \(t=0.0\), showing the transition from an initial uniform mesh to its evolution over five refinement steps. The NOV at each iteration are \(944\), \(1685\), \(2819\), \(4827\), \(8733\), and \(15908\), respectively. 25 presents the numerical results and corresponding meshes at four different time instances, demonstrating the adaptive method’s capability to efficiently handle time-dependent 3D problems with high precision. The left panel of 22 provides numerical evidence that the adaptive algorithm achieves the quasi-optimal convergence rate for the gradient error.

We set a relatively large mesh refinement tolerance, \(eTol = 2.5\), to reduce the computational cost. This choice leads to an earlier termination of the adaptive procedure, prior to the least-squares fitting step, while still allowing us to assess the effectiveness of both the adaptive algorithm and the neural network training strategy. The evolution of the training loss across time levels is presented in 26. Unlike the previous examples, the loss in this case decreases only to the order of \(10^{-4}\) for most time levels beyond the initial step. This behavior is primarily due to the inherent noise in the training data resulting from the relatively large adaptive tolerance, which limits further convergence. Nevertheless, the training loss remains well below the prescribed tolerance, indicating that the achieved accuracy is sufficient to reliably advance the computation. 27 illustrates the adaptive process of the initial mesh over four refinements, listing the NOV at each iteration as \(8080\), \(14502\), \(27327\), \(56739\), and \(128527\), respectively. 28 presents the evolution of the numerical solutions and corresponding meshes at four distinct time points, with refinement specifically focusing on the contracting spherical shell. As illustrated in the middle panel of 22, the adaptive algorithm attains the quasi-optimal convergence rate for the gradient error.

With the mesh refinement tolerance \(eTol = 0.1\), 30 shows the adaptive process at \(t=0\), terminating after seven steps and NOV at each iteration being \(390\), \(615\), \(1062\), \(2085\), \(4606\), \(10860\), and \(27113\), respectively. The mesh refinement process highlights that the adaptive algorithm achieves a highly efficient convergence procedure compared to conventional adaptive methods. 29 presents the training loss at various time levels during the learning process. To reach the desired MSE, the initial time level requires more than 200 training epochs. In contrast, fewer than 20 epochs are sufficient at all subsequent time levels to achieve a comparable level of accuracy. 31 illustrates the evolution of the numerical solutions and corresponding meshes at four distinct time points, capturing the symmetric splitting of the initial peak into two separate regions moving in opposite directions over time. The quasi-optimal convergence rate of the adaptive algorithm for the gradient error is demonstrated in the right panel of 22.

4.3 Additional benchmark examples↩︎

As a further demonstration beyond the linear model equations on rectangular domains considered above, we apply the proposed method to benchmark tests involving nonlinear dynamics and non-rectangular geometry.

Conclusion↩︎

In this paper, we proposed a neural network-enhanced \(hr\)-adaptive finite element method for parabolic equations that overcomes key limitations of conventional adaptive approaches, including excessive refinement iterations, costly interpolation between non-nested meshes, and the accumulation of unnecessary degrees of freedom across time steps. By employing a mesh-free neural network surrogate to replace the traditional interpolation operator, together with a controllable mesh size field strategy and a mesh-reset mechanism, the algorithm significantly reduces computational overhead while maintaining high accuracy. The combined error estimator ensures effective tracking of evolving singularities, and the adaptive process terminates within \(7\) iterations per time step. Numerical experiments in both \(2D\) and \(3D\) demonstrate the efficiency, robustness, and quasi-optimal convergence behavior of the proposed method, highlighting the potential of integrating neural networks as specialized modules within classical finite element frameworks.

In this work, we assume and numerically verify that the chosen network capacity is sufficient for the model problems under consideration, which are already challenging for conventional AFEM. This provides a reasonable starting point for the development of a complete and practical algorithmic framework. Nevertheless, several important questions remain open for future investigation, particularly the identification of minimal or near-optimal network architectures that can achieve the desired approximation accuracy with reduced computational complexity.

A potential limitation of the current approach arises when the complexity of the solution increases substantially during the simulation, or when highly localized and multiscale features emerge. In such situations, a fixed network architecture may no longer provide sufficient approximation power. This challenge is closely related to the well-known spectral bias of neural networks, which can make the representation of increasingly complex high-frequency components more difficult. Although such behavior was not observed in the numerical experiments presented in this work, it may become relevant for more demanding applications. One possible remedy is adaptive network enrichment. For example, the network capacity could be increased whenever a suitable error indicator exceeds a prescribed tolerance over a sustained period. The parameters associated with the existing network could be inherited from the previously trained model, while newly introduced parameters could be initialized independently. Another possible strategy is to employ alternative neural network architectures. In our follow-up work on neural network-based interpolation, we compared several neural network architectures for representing finite element solutions, including feedforward neural networks and extreme learning machines [47]. Our results indicate that certain architectures can achieve better approximation performance while maintaining low computational complexity. Consequently, selecting a more suitable network architecture may further improve the efficiency and robustness of the proposed framework. The design of reliable enrichment criteria, together with a systematic investigation of adaptive network-capacity strategies and alternative neural network architectures for problems with strongly localized or multiscale features, lies beyond the scope of the present work and will be explored in future research.

Acknowledgments↩︎

The authors express sincere thanks to the referees for their useful comments and suggestions, which led to improvements of the presentation. The authors thank Chunyu Chen for discussions about the implementation of the adaptive algorithm.

Huang and Yi were supported by the National Key R & D Program of China (2024YFA1012600) and NSFC Project (12431014). Yin was supported by the University of Texas at El Paso Startup Award.

5 Functions in neural network 6↩︎

In this section, we provide the distance functions \(d(\mathbf{x})\) and the continuous extension \(\tilde{g}(\mathbf{x},t_{n-1})\) of the boundary data used in the neural network 6 for the examples in 4.

2D squared domain. For the two-dimensional examples on \(\Omega=[-1,1]^2\), we take \[d(\mathbf{x}) = (x^2-1)(y^2-1).\] The continuous extension \(\tilde{g}(\mathbf{x},t_{n-1})\) in 6 is defined by \[\begin{align} \tilde{g}(\mathbf{x},t_{n-1}) = &\;\frac{1-x}{2}u(-1,y,t_{n-1}) + \frac{1+x}{2}u(1,y,t_{n-1}) + \frac{1-y}{2}u(x,-1,t_{n-1}) + \frac{1+y}{2}u(x,1,t_{n-1}) \nonumber \\ &- \frac{(1-x)(1-y)}{4}u(-1,-1,t_{n-1}) - \frac{(1-x)(1+y)}{4}u(-1,1,t_{n-1}) \nonumber \\ &- \frac{(1+x)(1-y)}{4}u(1,-1,t_{n-1}) - \frac{(1+x)(1+y)}{4}u(1,1,t_{n-1}). \end{align}\]

2D L-shaped domain. For the two-dimensional L-shaped domain \(\Omega = [-1,1]^2 \setminus \big([0,1]\times[-1,0]\big)\), we choose \[\begin{align} d_1(x,y) &= (1-x)(1+x)y(1-y),\\ d_2(x,y) &= (1+x)(-x)(1+y)(1-y), \end{align}\] and define \[\begin{align} d(x,y) = d_1(x,y)+d_2(x,y) +\sqrt{d_1^2(x,y)+d_2^2(x,y)}. \end{align}\] Equivalently, \[\begin{align} d(x,y) &= (1-x)(1+x)y(1-y) +(1+x)(-x)(1+y)(1-y) \\ &\quad + \sqrt{ \left[(1-x)(1+x)y(1-y)\right]^2 +\left[(1+x)(-x)(1+y)(1-y)\right]^2 }. \end{align}\] The continuous extension \(\tilde{g}(x,y,t_{n-1})\) is given by \[\begin{align} \tilde{g}(x,y,t_{n-1}) &= \frac{x(x-1)}{2}\,u(-1,y,t_{n-1}) + (1-x^2)\,u(0,y,t_{n-1}) + \frac{x(x+1)}{2}\,u(1,y,t_{n-1}) \\ &\quad + \frac{y(y-1)}{2} \Bigg[ u(x,-1,t_{n-1}) - \frac{x(x-1)}{2}\,u(-1,-1,t_{n-1}) \\ &\qquad\qquad - (1-x^2)\,u(0,-1,t_{n-1}) - \frac{x(x+1)}{2}\,u(1,-1,t_{n-1}) \Bigg] \\ &\quad + (1-y^2) \Bigg[ u(x,0,t_{n-1}) - \frac{x(x-1)}{2}\,u(-1,0,t_{n-1}) \\ &\qquad\qquad - (1-x^2)\,u(0,0,t_{n-1}) - \frac{x(x+1)}{2}\,u(1,0,t_{n-1}) \Bigg]. \end{align}\]

3D domain. For the three-dimensional examples on \(\Omega=[-1,1]^3\), we choose \[\begin{align} d(\mathbf{x}) = (1-x^2)(1-y^2)(1-z^2). \end{align}\] The continuous extension \(\tilde{g}(\mathbf{x},t_{n-1})\) is given by \[\begin{align} \tilde{g}(\mathbf{x},t_{n-1}) &= \frac{1-x}{2}\,u(-1,y,z,t_{n-1}) + \frac{1+x}{2}\,u(1,y,z,t_{n-1}) \\ &\quad + \frac{1-y}{2}\,u(x,-1,z,t_{n-1}) + \frac{1+y}{2}\,u(x,1,z,t_{n-1}) \\ &\quad + \frac{1-z}{2}\,u(x,y,-1,t_{n-1}) + \frac{1+z}{2}\,u(x,y,1,t_{n-1}) \\ &\quad - \frac{(1-x)(1-y)}{4}\,u(-1,-1,z,t_{n-1}) - \frac{(1-x)(1+y)}{4}\,u(-1,1,z,t_{n-1}) \\ &\quad - \frac{(1+x)(1-y)}{4}\,u(1,-1,z,t_{n-1}) - \frac{(1+x)(1+y)}{4}\,u(1,1,z,t_{n-1}) \\ &\quad - \frac{(1-x)(1-z)}{4}\,u(-1,y,-1,t_{n-1}) - \frac{(1-x)(1+z)}{4}\,u(-1,y,1,t_{n-1}) \\ &\quad - \frac{(1+x)(1-z)}{4}\,u(1,y,-1,t_{n-1}) - \frac{(1+x)(1+z)}{4}\,u(1,y,1,t_{n-1}) \\ &\quad - \frac{(1-y)(1-z)}{4}\,u(x,-1,-1,t_{n-1}) - \frac{(1-y)(1+z)}{4}\,u(x,-1,1,t_{n-1}) \\ &\quad - \frac{(1+y)(1-z)}{4}\,u(x,1,-1,t_{n-1}) - \frac{(1+y)(1+z)}{4}\,u(x,1,1,t_{n-1}) \\ &\quad + \frac{(1-x)(1-y)(1-z)}{8}\,u(-1,-1,-1,t_{n-1}) + \frac{(1-x)(1-y)(1+z)}{8}\,u(-1,-1,1,t_{n-1}) \\ &\quad + \frac{(1-x)(1+y)(1-z)}{8}\,u(-1,1,-1,t_{n-1}) + \frac{(1-x)(1+y)(1+z)}{8}\,u(-1,1,1,t_{n-1}) \\ &\quad + \frac{(1+x)(1-y)(1-z)}{8}\,u(1,-1,-1,t_{n-1}) + \frac{(1+x)(1-y)(1+z)}{8}\,u(1,-1,1,t_{n-1}) \\ &\quad + \frac{(1+x)(1+y)(1-z)}{8}\,u(1,1,-1,t_{n-1}) + \frac{(1+x)(1+y)(1+z)}{8}\,u(1,1,1,t_{n-1}). \end{align}\]

References↩︎