# Isogeometric analysis based mesh adaptation for time dependent problems

Mustapha Bahari, Abderrahmane Habbal, Ahmed Ratnani

► To cite this version:

Mustapha Bahari, Abderrahmane Habbal, Ahmed Ratnani. Isogeometric analysis based mesh adaptation for time dependent problems. Engineering with Computers, 2024, 40 (6), pp.3851-3874. 10.1007/s00366-024-02009-8. hal-04892403

HAL Id: hal-04892403

https://inria.hal.science/hal-04892403v1

Submitted on 16 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Isogeometric Analysis based Mesh Adaptation for Time Dependent Problems

Mustapha Bahari<fnref n="1" />,<fnref n="2" />*, Abderrahmane Habbal<fnref n="2" />,<fnref n="1" /> and
Ahmed Ratnani<fnref n="1" />

*Corresponding author(s). E-mail(s): mustapha.bahari@um6p.ma;

Contributing authors: abderrahmane.habbal@univ-cotedazur.fr;

ahmed.ratnani@um6p.ma;

## Abstract

This article presents a new algorithm designed to create a dynamic r-adaptive mesh within the framework of isogeometric analysis. The approach is based on the simultaneous computation of adaptive meshes using a nonlinear parabolic Monge-Ampere equation simultaneously with a resolution of partial differential equations in multidimensional spaces. The technique ensures the absence of geometric boundary errors and is simple to implement, requiring the solution of only one Laplace scalar equation at each time step. It utilizes a fast diagonalization method that can be adapted to any dimension. Various numerical experiments were conducted to validate an original parabolic Monge-Ampere solver. The solver was respectively applied to Burgers, Allen-Cahn, and Cahn-Hilliard problems to demonstrate the efficiency of the new approach.

**Keywords:** Optimal Transport, Isogeometric Analysis, moving mesh, r-refinement, parabolic Monge-Ampere equation, Burger's equation, Allen-Cahn equation, Cahn-Hilliard equation.

<footnote n="1">¹The UM6P Vanguard Center, Mohammed VI Polytechnic University, Lot 660, Hay Moulay Rachid, Benguerir, 43150, Morocco.</footnote>

<footnote n="2">²Université Côte d'Azur, Inria, CNRS, LJAD UMR 7351, Parc Valrose, Nice, 06108, France.</footnote>
# 1 Introduction

Isogeometric Analysis (IGA) is a well established finite element methodology introduced by Hughes et al. [1]. A distinctive feature of IGA is its utilization of exact geometry, which eliminates errors arising from poor geometric approximations. IGA also achieves higher regularity and significantly reduces the number of degrees of freedom in the parameterization process [1, 2]. Mesh refinement plays a pivotal role in enhancing accuracy across all numerical methods. Well-known techniques such as knot insertion (h-refinement) and polynomial degree elevation are discussed within the context of IGA in [1], primarily operate on a global scale, posing challenges for local refinement.

From a physical perspective, solutions to many Partial Differential Equations (PDEs), such as the Cahn-Hilliard equation, develop sharp interfaces with small thicknesses [3]. Consequently, addressing the need for local moving mesh adaptation becomes crucial. Achieving dynamic adaptive mesh in response to the time-dependent solution poses significant challenges. Several factors must be carefully considered, including the physical properties inherent in the solution, the accuracy, and the computational cost associated with mesh refinement and coarsening. Therefore, the initial concept of dynamic (moving/remeshing) adaptive mesh in IGA uses the initial parameterization as the coarsest mesh and subsequently refines it through existing local h-refinement methods.

Therefore, techniques for local h-refinement have been widely developed. Such as hierarchical B-splines [4], T-splines [5], LR-splines [6], and truncated T-splines [7], these methods increase the degrees of freedom as it is for local degree elevation [8]. These methods leverage the characteristic of preserving geometry during refinement. Other approaches use the well-known r-refinement method, which attempts to respond to the need for local refinement by displacing a mesh into the area of interest while maintaining the same degrees of freedom [9–13].

Moving mesh method in IGA is an active research field, recent paper [14] introduces a new algorithm for moving local mesh based on analysis-suitable T-splines. The notion of analysis-suitable is discussed in [15], providing T-splines with some necessary properties. Otherwise, to compute a solution at time $t^{n+1}$, the solution at time $t^n$ must be projected onto the new finite element space (mesh) constructed at $t^{n+1}$. This projection process can result in the loss of certain physical properties of the approximate solution. To address this, P. Hennig et al. [16] introduced novel projection and transfer operators. Mathematical studies have investigated stability [17, 18] and error estimation [17, 19]. The literature on this topic is expanding rapidly, and for further details, we refer the reader to [20–23]. Another method worth mentioning is the Remeshing method, which involves generating/eliminating sets of elements (see [24]).

Capturing boundary layers and other concentration phenomena, which are highly sensitive to parameterization quality, represents a significant challenge, particularly given that the boundary is typically defined by Computer-Aided Design (CAD). Consequently, methods based on the displacement of a control point face limitations because of the fixed external control points required for exact geometry representation. The authors have previously applied this approach, building upon the optimization algorithm initially introduced by Xu et al. [9, 10]. For alternative concepts, such as the
configurational force-based r-adaptive method, please refer to [11, 12]. Additionally, Ye Ji et al. [13] introduce a different approach based on the absolute principal curvature.

Therefore, we introduce a new approach building upon our recent work for constructing r-adaptive meshes on more general geometries [25]. This approach offers several advantages, including the utilization of the fast diagonalization method [25–27], employing a single solver for all geometries, bypass geometrical complexity, and requires only $C^0$ regularity for the initial mapping.

In this paper, we introduce a novel solver based on the Parabolic Monge-Ampere Equation (PMAE) [28, 29]. The aim is to handle the necessary refinement in computational processes by concurrently computing adaptive meshes and the solutions to PDEs, while preserving critical properties such as scaling structure and achieving boundary mesh adaptation.

The expense of our method arises from composing two mappings: the initial mapping and the mappings that result from resolving PMAE in each time step, which we refer to as the optimal mappings. While it is feasible to approximate the composition, our paper simply assesses B-spline basis functions in the quadrature and its image using the optimal mapping.

The paper is structured as follows. In section 2, we introduce B-splines and the finite element spaces. In section 3.1, we describe the basic idea behind the automatic boundary preservation under r-refinement, and we present a new fast parabolic Monge-Ampere solver. In section 4, we develop a new scheme for time-dependent parabolic partial differential equations in an adaptive mesh. Finally, in the last section 5.3, we perform a numerical accuracy study to illustrate the performance gain of the new approach.

## 2 B-Spline basis functions and Finite element spaces

B-spline functions have advantageous mathematical properties [1, 30] and are easy to implement due to their recursive construction, De Boor [31]. The introduction of B-spline functions begins with a sequence of monotonically increasing real numbers {$\xi_i$}$_i$ with $\xi_i \le \xi_{i+1}$, and a positive integer $p$ that represents the polynomial spline degree.

### B-splines.

Given a positive integer $n$, we introduce the so called knot vector

$$ \Xi = \{\xi_1, \xi_2, \dots, \xi_{n+p+1}\}, \qquad (2.1) $$

where $n$ is the number of basis functions. p-Open knot vector means that the first and last knot points of $\Xi$ have multiplicity $p+1$ and multiplicity $\mu = \{\mu_2, \dots, \mu_{m-1}\}$ for the others:

$$ \Xi = \{\underbrace{\zeta_1, \dots, \zeta_1}_{p+1}, \dots, \underbrace{\zeta_j, \dots, \zeta_j}_{\mu_i}, \dots, \underbrace{\zeta_m, \dots, \zeta_m}_{p+1}\}. \qquad (2.2) $$

In this paper, we consider a uniform p-open knot vector, i.e., $\mu_i = 1$ for $i = 2, \dots, m-1$, where $\zeta_i$ are equally spaced. Therefore, the univariate B-spline basis functions $B_i^p$,
i = 1, ..., n are defined recursively from the vector Ξ and degree p by the well known Cox-de Boor recursion formula:

$$B_i^0(\xi) = \begin{cases} 1 & \text{if } \xi \in [\xi_i, \xi_{i+1}), \\ 0 & \text{otherwise,} \end{cases} \qquad (2.3)$$

and,

$$B_i^p(\xi) = \frac{\xi - \xi_i}{\xi_{i+p} - \xi_i} B_i^{p-1}(\xi) + \frac{\xi_{i+p+1} - \xi_i}{\xi_{i+p+1} - \xi_{i+1}} B_i^{p-1}(\xi), \qquad (2.4)$$

where in (2.4), we adopt the convention 0/0 = 0. These basis functions are piecewise polynomials of degree p on the subdivision {$\zeta_1, \dots, \zeta_m$}, where they have p-1 continuous derivative in $(\zeta_1, \zeta_m)$ and other properties, like positivity, minimal compact support, the partition of unity, local linear independence, etc.

An example of quadratic basis functions and their linear combination "quadratic B-spline curve" is shown in Fig.1. Note that the multiplicity p + 1 in the first (resp. last) knot leads to interpolate the first (resp. last) control point.

Respectively, a multi-dimensional B-spline can be constructed as a tensor product $B_{j_1, j_2, \dots, j_d}^{p_1, p_2, \dots, p_d} = \otimes_{i=1}^d B_{j_i}^{p_i}$ [1, 2].

<figure><img src="image_2.png" /><figcaption>**Fig. 1:** (left) A quadratic B-Spline curve and its control points using the knot vector Ξ = {$000 \quad \frac{1}{5} \quad \frac{2}{5} \quad \frac{3}{5} \quad \frac{4}{5} \quad 111$}, (right) the corresponding B-Splines.</figcaption></figure>

**Definition 1** (B-spline mapping). The B-spline mapping in $\mathbb{R}^2$ associated to Ξ = (Ξ₁, Ξ₂), p = (p₁, p₂), and a set of control points P, is described by:

$$\mathbf{F} = \sum_{j_1, j_2} \mathbf{P}_{j_1, j_2} B_{j_1, j_2}^{p_1, p_2} \qquad (2.5)$$
where $\mathbf{P}_{j_1,j_2} \in \mathbb{R}^2$.

## 2.1 Finite element spaces

From the vector knots $\Xi$, let us introduce a mesh $\mathcal{Q}_h$, a partition family on $\Omega_c$ by

$$ \mathcal{Q}_h = \{Q = \otimes_{i=1}^{d}(\zeta_{j_i,i}, \zeta_{j_i+1,i}), 1 \le j_i \le m_i - 1, d = 2, 3\}, \quad (2.6) $$

we denote by $h_Q$ the diameter of the element $Q \in \mathcal{Q}_h$. Each element in $\mathcal{Q}_h$ is mapped into another element of the physical domain, which introduces a mesh $\mathcal{K}_h$ in the physical domain $\Omega_p$, defined by:

$$ \mathcal{K}_h := \{K = \mathbf{F}(Q), Q \in \mathcal{Q}_h\}. \quad (2.7) $$

Analogously $\tilde{Q}$, the support extension of $Q$, is mapped into $\tilde{K} = \mathbf{F}(\tilde{Q})$. We introduce the B-spline space as follows:

$$ S_{\mu}^{p} = \operatorname{span}\{B_{\mu, i}^{p}, \text{ for } i = 1, \dots, n\} \quad (2.8) $$

where $B_{\mu, i}^p$ is a B-spline $B_i^p$ with regularity $\mu$. The induced discrete subspace of $H^s(\Omega_c)$ is defined by

$$ \mathcal{V}_h = \otimes_{i=1}^{d} S_{\mu_i}^{p_i}. \quad (2.9) $$

and the induced discrete subspace of $H^s(\Omega_p)$ is given by

$$ \mathcal{V}_h \circ \mathbf{F}^{-1} = \operatorname{span}\{\otimes_{j=1}^{d} B_{\mu_j, i_j}^{p_j} \circ \mathbf{F}^{-1}, \text{ for } i_j = 1, \dots, n, j = 1, \dots, d\}. \quad (2.10) $$

and we note that,

$$ \left\{ \frac{dv}{d\xi} : v \in S_{\mu}^{p} \right\} = S_{\mu-1}^{p-1}. \quad (2.11) $$

We introduce the discrete subspaces $V_h(\operatorname{div}, \Omega_c)$ and $V_h(L^2, \Omega_c)$ of $H(\operatorname{div}, \Omega_c)$ and $L^2(\Omega_c)$ respectively, as follows:

$$ \left\{ \begin{array}{l} V_h(\operatorname{div}, \Omega_c) = \left( S_{\mu_1}^{p_1} \otimes S_{\mu_2-1}^{p_2-1} \right), \\ V_h(L^2, \Omega_c) = S_{\mu_1-1}^{p_1-1} \otimes S_{\mu_2-1}^{p_2-1}, \end{array} \right. \quad (2.12) $$

and for $(\Psi^h, \Psi^h) \in \mathcal{V}_h \times V_h(\operatorname{div}, \Omega_c)$

$$ \left\{ \begin{array}{l} \nabla \times \Psi^h = \left( \frac{\partial \Psi^h}{\partial y} - \frac{\partial \Psi^h}{\partial x} \right), \\ \\ \nabla \cdot \Psi^h = \frac{\partial \Psi^h}{\partial x} + \frac{\partial \Psi^h}{\partial y}. \end{array} \right. $$
The choice of the discrete spaces $V_h(\mathrm{div}, \Omega_c)$ and $V_h(L^2, \Omega_c)$ is made to ensure that they form the DeRham sequence [32, 33]. In the 2D case, we have:

$$V_h(\mathrm{grad}, \Omega_c) := V_h \xrightarrow{\nabla \times} V_h(\mathrm{div}, \Omega_c) \xrightarrow{\nabla \cdot} V_h(L^2, \Omega_c). \quad (2.13)$$
### 3 Parabolic Monge-Ampere solver for boundary preservation under r-refinement

Consider a simplified physical domain, denoted $\Omega_p$, as shown in Fig. 2-(top). The first step of Isogeometric Analysis involves seeking a bijective, continuous, and one-to-one mapping, called initial mapping, denoted $F$, which establishes a connection between the computational domain $\Omega_c$ and the physical domain $\Omega_p$.

$$F : \Omega_c \rightarrow \Omega_p. \qquad (3.1)$$

Let $\rho_1$ represent a density function describing the desired adaptive mesh within the physical domain. Technique illustrated in Fig. 2, demonstrate that the adaptive mesh within the physical domain, corresponds to $\rho_1$, is equivalent to the image under initial mapping of an adaptive mesh equidistributed according to $\rho_1 \circ F$ in $\Omega_c$. For a comprehensive understanding, refer to subsection 3.2 [25].

<figure><img src="image_3.png" /><figcaption>**Fig. 2:** Illustration of the adaptive mesh technique and Illustration of uniform and adaptive meshes associated with each density function.</figcaption></figure>

The technique procedure can be summarized in two steps. Firstly, within $\Omega_c$, the equidistributed mesh is computed with respect to $\rho_1 \circ F$, yielding the adaptive one-to-one mapping $\Psi$. Secondly, the composition $F' \equiv F \circ \Psi$ is employed as a new mapping (in simulation) to generate the adaptive mesh within the physical domain. Therefore, for $\rho_0 = 1$ corresponds to a uniform mesh in $\Omega_c$, and

$$\rho = \rho_1 \circ F \qquad (3.2)$$
where $\sigma = |\Omega_c| / \int_{\Omega_c} \rho_1 \circ F(x) dx$. We are looking for the mapping $\Psi$ which is the gradient of a solution to the following Monge-Ampere system

$$ \left\{ \begin{array}{ll} \rho(\nabla\phi) \det(H(\phi)) = 1 & \text{in } \Omega_c, \\ \nabla\phi \cdot \mathbf{n} = \xi \cdot \mathbf{n} & \text{for all } \xi \in \partial\Omega_c, \\ \phi \text{ is convex.} & \end{array} \right. \qquad (3.3) $$

where $\mathbf{n}$ is the unit outward normal vector on $\Omega_c$. Therefore, the resolution remains consistent regardless of the geometry's complexity, as the initial mapping is primarily used in composition with the density function $(\rho_1 \circ F)$.

## 3.1 Parabolic Monge-Ampere solver for mesh generation

Several numerical techniques have been proposed to solve the Monge-Ampere equation (3.3). For a concise summary, please refer to [25, 34]. The parabolic Monge-Ampere equation was used in the literature for mesh equidistribution [28, 29, 35], where the solution of (3.3) represents the steady state of PMAE. The reason for choosing the parabolic Monge-Ampere equation is its natural integration with time-dependent PDEs. This choice is also motivated by the fact that we do not aim to solve the Monge-Ampere equation exactly; instead, our goal is to generate adaptive meshes closely following the behavior of such a monitor function. The PMAE is considered as relaxation technique [29], taking the form

$$ \left\{ \begin{array}{ll} \epsilon(I - \gamma\Delta)\partial_t\phi = (\rho(\nabla\phi)\det(H(\phi)))^{\frac{1}{d}} & \text{in } \Omega_c, \\ \nabla\phi \cdot \mathbf{n} = \xi \cdot \mathbf{n} & \text{for all } \xi \in \partial\Omega_c, \\ \phi(\xi, t=0) = \frac{\|\xi\|_2^2}{2} & \text{for all } \xi \in \Omega_c, \end{array} \right. \qquad (3.4) $$

The parameters $\epsilon$ and $\gamma$ play an important role in adjusting the response time of mesh adaptation in $\rho$. Specifically, $\gamma$ governs the degree of spatial smoothing, while $\epsilon$ controls the relaxation rate. Additional information about the relaxation and smoothing operator $(I - \gamma\Delta)$ can be found in [29, 36].

Therefore, the ability to preserve the boundary arises from the fact that $\Psi = \nabla\phi$ is not only injective but indeed a bijective mapping from the computational domain to itself. This is handled using the exact formulation at the boundary of $\Omega_c$, $\Psi \cdot \mathbf{n} = \xi \cdot \mathbf{n}$. The methodology of computing the adaptive mapping $\Psi$ with boundary preservation property is addressed in the next paragraph.

## Hdiv formulation of Parabolic Monge-Ampere solution

The Neumann boundary condition in the system (3.4) makes the mesh move along the boundary but not away from it. However, those boundary errors become bigger when a density function has a variation near the boundary. Therefore, a mixed variational formulation is employed to correct the gradient mapping $\Psi = \nabla\phi$ within $\Omega_c$, described as follows:

$$ \left\{ \begin{array}{ll} \Psi = \nabla\phi & \text{in } \Omega_c, \\ \Psi \cdot \mathbf{n} = \xi \cdot \mathbf{n} & \text{for all } \xi \in \partial\Omega_c. \end{array} \right. \qquad (3.5) $$
## 3.2 Variational formulation and PMAE Algorithms

Let us begin by introducing the time-explicit discretization of equation (3.4), as follows:

$$ \epsilon(I - \gamma\Delta)\phi^{n+1} = \epsilon(I - \gamma\Delta)\phi^n + dt(\rho(\nabla\phi^n)\det(H(\phi^n)))^{\frac{1}{d}}, \quad (3.6) $$

The rationale behind employing explicit discretization is to leverage the fast diagonalization method [25–27]. This choice is driven by the fact that, as outlined in section 3, the PMAE will always be solved within the computational domain.

The variational formulation of (3.4) is given by
$$ \text{Find } \phi^{n+1} \in H^2(\Omega_c) \text{ such that } \forall \varphi \in H^2(\Omega_c) $$

$$ \epsilon \int_{\Omega} \phi^{n+1} \varphi + \epsilon \gamma \nabla \phi^{n+1} \cdot \nabla \varphi = \int_{\Omega} \left( \epsilon(I - \gamma \Delta)\phi^n + dt(\rho(\nabla \phi^n) \det(H(\phi^n)))^{\frac{1}{d}} \right) \varphi \\ + \epsilon \gamma \int_{\partial \Omega} \boldsymbol{\xi} \cdot \mathbf{n} \varphi. \quad (3.7) $$

The variational formulation of (3.5) is given by introducing the decomposition $\Psi = \Psi_0 + \Psi_N$, where $\Psi_N$ is given such that $\Psi_N \cdot n = \boldsymbol{\xi} \cdot n$ on $\partial\Omega_c$ and the unknown $\Psi_0$ is such that $\Psi_0 \cdot n = 0$ on $\partial\Omega_c$. such that $\Psi_0$ is a solution of the following Galerkin approximation problem:

Find $\Psi_0 \in H_0(\text{div}, \Omega)$ such that $\forall v \in H_0(\text{div}, \Omega)$

$$ \int_{\Omega_c} \Psi_0 \cdot v = - \int_{\Omega_c} \partial \nabla \cdot v - \int_{\Omega_c} \Psi_N \cdot v. \quad (3.8) $$

In all computations, we have set the parameters $\epsilon$ and $\gamma$ to 1. To ensure time stability, when employing a density function with intensity lower than 10, we fixed the time steps at $dt = 0.5$. For details on stability analysis, refer to [29]. This chosen parameters have consistently proven to be efficient across various tests, demonstrating convergence to high-quality mesh adaptation.

The adaptive mapping is then given by the algorithm as follows:
**Algorithm 1:** Parabolic Monge-Ampère algorithm

**Require:** Given an initial value $\phi^0 = \frac{\|\xi\|_2^2}{2}$, and tolerance tol.

**Ensure:** Adaptive mapping $\Psi$.

1. Start time loop

2. While $|\phi^{n+1} - \phi^n|_{H^1(\Omega_c)} \ge tol$, do

1. Compute $\phi^{n+1}$ by solving the following variational formulation problem:
   Find $\phi^{n+1} \in H^2(\Omega_c)$ such that $\forall \varphi \in H^2(\Omega_c)$

$$ \int_{\Omega} \phi^{n+1} \varphi + \nabla \phi^{n+1} \cdot \nabla \varphi = \int_{\Omega} \left( (I - \Delta)\phi^n + 0.5(\rho(\nabla\phi^n)\det(H(\phi^n)))^{\frac{1}{d}} \right) \varphi \\ + \int_{\partial\Omega} \xi \cdot n\varphi. \qquad (3.9) $$

3. Compute $\Psi^{n+1}$ from (3.5) using the steady state solution $\phi^{n+1}$.

The in time adaptive mapping is given by the algorithm as follow:

**Algorithm 2:** In time Parabolic Monge-Ampere algorithm

**Require:** Given time dependent density $\rho_1(t, \cdot)$, time step $\delta t$ and the final time $t_{\max}$.

**Ensure:** Adaptive mapping $\Psi^n$ at $t_n$, $n \in \{0, \dots, t_{\max}\}$.

1. Compute the initial value $\phi^0$ at the initial time $t_0$ using Algorithm 1,

2. While $t \le t_{\max}$, do

1. $t_{n+1} = t_n + \delta t$

2. Compute $\phi^{n+1}$ from (3.9) using the last solution $\phi^n$.

3. Compute $\Psi^{n+1}$ from (3.5) using the last solution $\phi^{n+1}$.

**Remark 1.** 1. The mixed variational formulation (3.8) will be addressed through the fast diagonalization method [25].

2. Following the DeRham sequence (2.13) the Galerkin discrete solution $\Psi$ belongs to the discrete space $V_h(\mathrm{div}, \Omega_c)$. This means that $\phi \in V_h(L^2, \hat{\Omega})$. Specifically, if $\phi$ has a polynomial degree of $(p, p)$, a mapping $\Psi$ will have degrees $(p+1, p)$ and $(p, p+1)$.

3. Algorithm 2 is well-suited for time-dependent density function, which aligns with our utilization in the moving mesh method.

4. Our numerical code is written in Python, where intensive computations parts were accelerated using Pyccel [37], a static compiler for Python, based on Fortran and C. A laptop with a quad 2.6 GHz Intel Core i7, on ubuntu OS is used throughout the numerical tests.
### 3.3 Examples with an analytical density function

Next, we show the efficiency of mesh adaptation using analytical density functions. Let us define the global displacement of the grid points as follows:

$$ \mathcal{G}_{disp}^2 = \int_{\Omega_c} |\Psi - \xi|^2 \rho \, d\xi \qquad (3.10) $$

where $\sigma$ is defined in (3.2). We use a B-spline description for various models, employing quadratic B-splines. This implies that both $\phi$ and $\mathcal{F}$ have the same polynomial degree of $p = 2$ in each direction. The notation $\#\text{cells}$ stands for the number of cells in each direction.

**Example 1** (Quarter of a square with a circular hole). In this test, we show mesh adaptation using a density function defined as follows:

$$ \rho_1(x, y) = 2 + \sin(k^2 \pi \sqrt{(x - 0.6 - 0.25t)^2 + (y - 0.6)^2}) \qquad (3.11) $$

with a wave number of $k^2 = 4$. Fig. 3 shows the quality mesh adaptation at $t = 0$ using Algorithm 1. The mesh's behavior in the upper left-hand corner is caused by the geometrical singularity, which arises from a repeated control point in the initial mapping. Table 1 shows details on CPU time and the min-max values of the Jacobian function. Fig. 4 provides an example of how the mixed variation formulation ensures the exactness of the geometries. Fig. 5 illustrates the dynamic mesh refinement response using a time-dependent density function. Here, we use the Algorithm 2, with time step $\delta t = 10^{-2}$.

<table><caption>**Table 1:** Example 1: CPU-time needed to solve the Monge-Ampere equation and the global displacement $\mathcal{G}_{disp}$.</caption><thead><tr><th>#cells</th><th>G<sub>disp</sub></th><th>CPU-time (s)</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>1.407e-02</td><td>0.113</td><td>5.149e-01</td><td>2.198e+00</td></tr><tr><td>32</td><td>1.391e-02</td><td>0.351</td><td>6.325e-01</td><td>2.022e+00</td></tr><tr><td>64</td><td>1.383e-02</td><td>1.605</td><td>6.516e-01</td><td>1.984e+00</td></tr><tr><td>128</td><td>1.380e-02</td><td>8.784</td><td>6.535e-01</td><td>1.967e+00</td></tr></tbody></table>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3: Example 1: (left) Resulting adaptive mesh and (right) The monitor function at t = 0.</figcaption></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4: Example 1: A close-up on (left) resulting adaptive mesh using gradient mapping $\Psi = \nabla \phi$, and (right) using mixed variational formulation $\Psi \in V_h(\text{div}, \Omega_c)$.</figcaption></figure>

**Example 2 (Circular domain).** The example corresponds to a vortex located at (0, 0), defined by:

$$ \rho(\xi, \eta) = 1 + \kappa / \left( 1 + (10r(\xi, \eta) \cos(\theta(\xi, \eta) - 20r(\xi, \eta)))^2 \right), \quad (3.12) $$

where $\theta(\xi, \eta) = \arctan(\eta/\xi)$, and $r(\xi, \eta) = \sqrt{\xi^2 + \eta^2}$. Fig. 6 shows the high quality of mesh adaptation on the circle. Table 2 shows CPU-time and the positive Jacobian function with high accuracy. The objective in this example is to demonstrate that in PMAE, two factors can influence the mesh adaptation intensity: the intensity of a density function $\kappa$ and the time step $\text{dt}$. Figure 7 shows the intensity of mesh adaptation, measured by the difference between the maximum and minimum values of the Jacobian function, where a big value corresponds to high mesh adaptation intensity.
<figure><img src="image_8.png" /><figcaption>Fig. 5: Example 1: (top) Adaptive mesh and (bottom) corresponded monitor functions at t = 1, 3, 5.</figcaption></figure>

<table><caption>Table 2: Example 2: CPU-time needed to solve the Monge-Ampere equation and the global displacement $\mathcal{G}_{disp}$ using degree $p = 3$.</caption><thead><tr><th>#cells</th><th>G<sub>disp</sub></th><th>CPU-time (s)</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td>16</td><td>8.203e-02</td><td>0.237</td><td>1.075e-01</td><td>2.221e+00</td></tr><tr><td>32</td><td>8.394e-02</td><td>0.918</td><td>8.431e-02</td><td>3.044e+00</td></tr><tr><td>64</td><td>8.114e-02</td><td>3.828</td><td>1.022e-01</td><td>2.901e+00</td></tr><tr><td>128</td><td>7.727e-02</td><td>16.08</td><td>1.815e-01</td><td>2.905e+00</td></tr></tbody></table>

**Example 3** (three dimensional case). In this test, we demonstrate the effectiveness of mesh adaptation for three dimensional case using a density function defined as follows:

$$ \rho_1(x, y, z) = 7/(2 + \cos(8\pi\sqrt{(x-0.5)^2 + (y-0.5)^2 + (z-0.5)^2})) \quad (3.13) $$

Figs. 9-10 shows the high-quality mesh adaptation. Table 3 shows the global grid displacement, CPU-time, and Jacobian function values.
<figure><img src="image_9.png" /><figcaption>Fig. 6: Example 1: (left) Resulting adaptive mesh and (right) The monitor function at $\kappa = 9$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7: Example 1: (left) resulting adaptive mesh intensity with respect to $\kappa$ using $dt=0.5$, and (right) with respect to $dt$ using $\kappa=9$.</figcaption></figure>

<table><caption>Table 3: Example 3: CPU-time needed to solve the Monge-Ampere equation and the global displacement $\mathcal{G}_{disp}$.</caption><thead><tr><th>Geometry</th><th>#cells</th><th>G<sub>disp</sub></th><th>CPU-time (s)</th><th>min Jac(Ψ)</th><th>max Jac(Ψ)</th></tr></thead><tbody><tr><td rowspan="3">Cylinder</td><td>16</td><td>7.851e-01</td><td>9.997</td><td>4.774e-01</td><td>1.933e+00</td></tr><tr><td>32</td><td>7.808e-01</td><td>45.041</td><td>5.153e-01</td><td>1.791e+00</td></tr><tr><td>64</td><td>7.813e-01</td><td>358.806</td><td>5.702e-01</td><td>1.763e+00</td></tr><tr><td rowspan="3">3D-Puzzle piece</td><td>16</td><td>5.200e-01</td><td>9.92</td><td>4.666e-01</td><td>1.898e+00</td></tr><tr><td>32</td><td>5.336e-01</td><td>44.514</td><td>5.429e-01</td><td>1.808e+00</td></tr><tr><td>64</td><td>5.410e-01</td><td>325.298</td><td>5.653e-01</td><td>1.730e+00</td></tr><tr><td rowspan="3">Quarter of cube</td><td>16</td><td>2.659e+00</td><td>10.068</td><td>4.069e-01</td><td>2.185e+00</td></tr><tr><td>32</td><td>2.713e+00</td><td>45.806</td><td>4.542e-01</td><td>1.879e+00</td></tr><tr><td>64</td><td>2.741e+00</td><td>357.558</td><td>5.584e-01</td><td>1.757e+00</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Fig. 8: Example 3:** (left) A close-up of one extract part inside the cylinder (right) Adaptive meshes and the monitor function.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9: Example 3:** (left) A close-up of one extract part inside the 3D-Puzzle piece (right) Adaptive meshes and the monitor function.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10: Example 3:** (left) A close-up of one extract part inside the Quarter of cube (right) Adaptive meshes and the monitor function.</figcaption></figure>
# 4 Nonlinear parabolic time dependent problems

In this section, we introduce a finite element scheme within the framework of IGA for solving PDEs on moving meshes. We provide a detailed explanation of the Galerkin discrete variational formulation of the PDEs within the adaptive mesh. Finally, we conclude by introducing the algorithm of the moving mesh method.

Our focus is on r-refinement mesh techniques for parabolic equations, which can be expressed as follows:

$$ \frac{\partial u}{\partial t} = f(t, x, u, \nabla u, \Delta u), \quad \forall x \in \Omega_p, t > 0. \qquad (4.1) $$

To address this equation within the framework of moving mesh in isogeometric analysis, we compute $u$ in the physical domain $\Omega_p$ defined by $\Omega_p = \mathbf{F}(\mathbf{\Psi}(\Omega_c))$, where $\mathbf{\Psi}$ is a time-dependent mapping computed from Algorithm 2. This is achieved by applying the chain rule to $u \circ (\mathbf{F} \circ \mathbf{\Psi})$, yielding:

$$ \frac{\partial (u \circ (\mathbf{F} \circ \mathbf{\Psi}))}{\partial t} = \frac{\partial u}{\partial t} \circ (\mathbf{F} \circ \mathbf{\Psi}) + \nabla \mathbf{F} \circ \mathbf{\Psi} \frac{\partial \mathbf{\Psi}}{\partial t} \cdot \nabla u \circ (\mathbf{F} \circ \mathbf{\Psi}), $$

hence, we derive an equivalent system with an additional quasi-Lagrangian term by applying $(\mathbf{F} \circ \mathbf{\Psi})^{-1}$ and substituting $\frac{\partial u}{\partial t}$, as follows:

$$ \frac{\partial (u \circ (\mathbf{F} \circ \mathbf{\Psi}))}{\partial t} \circ (\mathbf{F} \circ \mathbf{\Psi})^{-1} = f(t, x, u, \nabla u, \Delta u) + \left( \nabla \mathbf{F} \circ \mathbf{\Psi} \frac{\partial \mathbf{\Psi}}{\partial t} \right) \circ (\mathbf{F} \circ \mathbf{\Psi})^{-1} \cdot \nabla u. \quad (4.2) $$

This formulation avoids the need for projection into the new space (mesh) during the subsequent adaptive mapping induced from parabolic Monge-Ampere equation. This approach is commonly adopted in various papers such as [28, 29, 38], all sharing the concept of avoiding interpolation. This ensures that the mesh inherits advantageous dynamic properties of the solution, such as scaling structures [29, 39].

## 4.1 Galerkin discretization of PDEs on moving meshes

Let us now outline the discrete Galerkin approach of (4.2). For simplicity, we close the equation (4.2) using homogeneous Dirichlet boundary condition. Let $W_h^n$ denote the discrete Galerkin space at time instant $t^n$, defined as

$$ W_h^n = \mathcal{V}_h \circ (\mathbf{F} \circ \mathbf{\Psi}^n)^{-1} \cap H_0^1(\Omega_p), \qquad (4.3) $$

where $\mathcal{V}_h \circ (\mathbf{F} \circ \mathbf{\Psi}^n)^{-1}$ is given by (2.10), and $\mathbf{\Psi}^n$ represents the optimal mapping at $t^n$. Let us denote the basis function of $\mathcal{V}_h$ as $\phi_j = B_{\mu_1, j_1}^{p_1} B_{\mu_2, j_2}^{p_2}$; thus, the basis function of $W_h^n$ is given by $\phi_j \circ (\mathbf{F} \circ \mathbf{\Psi}^n)^{-1}$.
To start, we employ time discretization using an implicit-scheme in the following manner

$$ \frac{u^{n+1} - u_{\Psi^{n+1}}^n}{\delta t} = f(t^{n+1}, x, u^{n+1}, \nabla u^{n+1}, \Delta u^{n+1}) + \left( \nabla \mathbf{F} \circ \mathbf{\Psi}^{n+1} \frac{\mathbf{\Psi}^{n+1} - \mathbf{\Psi}^n}{\delta t} \right) \circ (\mathbf{F} \circ \mathbf{\Psi}^{n+1})^{-1} \cdot \nabla u^{n+1} \quad (4.4) $$

where $u_{\Psi^{n+1}}^n = u^n \circ (\mathbf{F} \circ \mathbf{\Psi}^n) \circ (\mathbf{F} \circ \mathbf{\Psi}^{n+1})^{-1}$, and $u^n \in W_h^n$ denotes the solution at $t^n$ and $\delta t$ is a time step. Therefore, the Galerkin weak formulation of (4.4) is given by Find $u^{n+1} \in W_h^{n+1}$ such that $\forall v \in W_h^{n+1}$

$$ B(u^{n+1}, v) = 0 \quad (4.5) $$

where

$$ B(u^{n+1}, v) = \int_{\Omega_p} u^{n+1} v \, dx - \int_{\Omega_p} u_{\Psi^{n+1}}^n v \, dx - \delta t \int_{\Omega_p} f(t^{n+1}, x, u^{n+1}, \nabla u^{n+1}, \Delta u^{n+1}) v \, dx \\ - \delta t \underbrace{\int_{\Omega_p} \left( \nabla \mathbf{F} \circ \mathbf{\Psi}^{n+1} \frac{\mathbf{\Psi}^{n+1} - \mathbf{\Psi}^n}{\delta t} \right) \circ (\mathbf{F} \circ \mathbf{\Psi}^{n+1})^{-1} \cdot \nabla u^{n+1} v \, dx}_{L_t} \quad (4.6) $$

It is crucial to emphasize that subsequent to the integration of the quasi-Lagrangian term, the computed function $u_{\Psi^{n+1}}^n$ belongs to $W_h^{n+1}$ in (4.6). Upon a change of variables in (4.6), for $v = \phi_i \circ \mathbf{F}^{-1} \circ (\mathbf{\Psi}^{n+1})^{-1}$, we obtain the Lagrangian term as follows:

$$ \begin{aligned} L_t &= \sum_j u_j^{n+1} \int_{\Omega_p} \left( \nabla \mathbf{F} \circ \mathbf{\Psi}^{n+1} \frac{\mathbf{\Psi}^{n+1} - \mathbf{\Psi}^n}{\delta t} \right) \circ (\mathbf{F} \circ \mathbf{\Psi}^{n+1})^{-1} \cdot \nabla (\phi_j \circ \mathbf{F}^{-1} \circ (\mathbf{\Psi}^{n+1})^{-1}) \\ &= \sum_j u_j^{n+1} \int_{\Omega_c} \nabla \mathbf{F} \circ \mathbf{\Psi}^{n+1} \frac{\mathbf{\Psi}^{n+1} - \mathbf{\Psi}^n}{\delta t} \cdot (\mathbf{J}_{ad})^{-T} \nabla (\phi_j) \phi_i \det(\mathbf{J}_{ad}) \, dx \\ &= \sum_j u_j^{n+1} \int_{\Omega_c} \frac{\mathbf{\Psi}^{n+1} - \mathbf{\Psi}^n}{\delta t} \cdot (\nabla \mathbf{\Psi}^{n+1})^{-T} \nabla (\phi_j) \phi_i \det(\mathbf{J}_{ad}) \, dx, \end{aligned} \quad (4.7) $$

where

$$ \mathbf{J}_{ad} := \nabla (\mathbf{F} \circ \mathbf{\Psi}) = (\nabla \mathbf{\Psi})^T \nabla \mathbf{F} \circ \mathbf{\Psi}. $$

## 4.2 Moving mesh solver

Our aim in this subsection is to present a new moving mesh solver. Here, we solve simultaneously the Galerkin weak formulation (4.5) alongside the In time Parabolic Monge-Ampere Algorithm 2.
Therefore, let us define the residual vector as $\mathcal{R}$ where

$$ \mathcal{R}(u^{n+1}) = \{B(u^{n+1}, \phi_i)\}_i, \quad (4.8) $$

we denote the Jacobian function of $\mathcal{R}$ by $\mathcal{R}'$, we iteratively solve the nonlinear system through the Newton method. Hence, the moving mesh algorithm is introduced as follows:

**Algorithm 3:** moving mesh algorithm

1. Given the initial solution $u^0 = u(t_0)$ and uniform mesh $Q_h$.

2. Compute $\phi^0$ and $\Psi^0$ from Algorithm (1) using a density function based on $u^0$.

3. Compute $u^0 \in W_h^0$ in adapted mesh.

4. Start the time loop:

(a) $t_{n+1} = t_n + \delta t$.

(b) Compute the density function based on $u^n$.

(c) Compute $\phi^{n+1}$ from (3.9) using the last solution $\phi^n$.

(d) Compute $\Psi^{n+1}$ from (3.5) using the last solution $\phi^{n+1}$.

(e) Map the quadrature points through $\Psi^{n+1}$ and subsequently compute the B-spline function at these mapped points.

(f) Set $u^{(0),n+1} = u^n$.

(g) Multicorrector stage: do loop on m

(i) Solve $\mathcal{R}'(u^{(m),n+1})\delta u^{(m)} = -\mathcal{R}(u^{(m),n+1})$.

(ii) Update $u^{(m+1),n+1} = u^{(m),n+1} + \delta u^{(m)}$.

(h) $u^{n+1} = u_h^{(m+1),n+1}$.

**Remark 2.** -If the initial solution does not represent any boundary interface or there is no need to begin with an adapted mesh, $\Psi^0$ will remain as identity, and step 3 will be unnecessary.

-We use a distinct time step $\delta t$ for the PDE (4.2), which operates independently of the parabolic Monge-Ampere equation (3.9). This separation is necessary because aligning the time steps ($dt = \delta t$) would require adjusting parameters $\epsilon$ and $\gamma$ for each $\delta t$.

**Density function**

Our investigation extends to devising an appropriate monitor function. We introduce a modified arc-length monitor function (4.9), which is designed to generate high-quality adaptive mesh. By using an average rather than directly using the solution itself ($\rho_1 = \sqrt{1+c_1|\nabla u|+c_2|u|}$, as commonly employed in various papers [40–42]), mesh adaptation can extend around the solution and anticipate changes in its vicinity for future solutions. It is crucial to note that the following monitor function is well chosen from the motivation that the solution has a sharp interface, which means we need a density that adapts the mesh at the solution interface. Otherwise, the regularity of the density function is crucial for generating smooth adaptive meshes.
Let $\rho_Q$ be a piecewise constant monitor function, defined by:

$$ \rho_Q = \int_Q |(\nabla \mathbf{F})^{-T} \nabla(u_h \circ \mathbf{F})| \, d\xi, \quad \forall Q \in \mathcal{Q}_h. \qquad (4.9) $$

Therefore, in order to control the regularity of a density function and its variation, we introduce a regularized monitor function as follows:

$$ \rho_1(\mathbf{F}(x)) = c_0 \Pi_p(\rho_Q) + c_1, \quad \forall x \in Q, \forall Q \in \mathcal{Q}_h, \qquad (4.10) $$

where $\Pi_p$ is the $L^2$-projection (or quasi interpolation), so that, $\rho_1 \in C^{p-1}$. Moreover, $Im(\rho_1)$ is bounded during the computation thanks to the following choice of $c_0$ and $c_1$

$$ c_0 = \frac{10}{\max(\rho_Q) - \min(\rho_Q)}, \quad c_1 = \frac{\max(\rho_Q) - 11 \min(\rho_Q)}{\max(\rho_Q) - \min(\rho_Q)}. $$

**Remark 3.** - The selected intensity of our constructed monitor function $\rho_1$ aligns with the requirements of the analytic density function presented in section 3.3.

- The notation $\rho_1 \circ \mathbf{F}$ in (4.10) is due to our change of variable within the integral formula (4.9). However, it is crucial to note that the $u_h$ referenced in (4.9) is derived from projecting the solution in adapted mesh into a uniform mesh, intended to be included in Algorithm 3.

# 5 Numerical Results

In this section, to assess the efficiency of our approach, we select several examples of manufactured solutions. Our emphasis in these tests lies in employing $C^1$ regularity through quadratic B-splines. We choose to perform the computation using the time step $\delta t = \frac{1}{7(\#cells)^2}$, while setting a tolerance for the Newton iteration at $||\delta u^{(m)}|| \le 10^{-10}$. Finally, we solve the induced linear system using a preconditioned GMRES algorithm from the scipy library (see Jones et al. [43]).

Whenever you find **un** it stands for uniform mesh, and **ad** stands for adaptive mesh.

## 5.1 Burgers' equation

We evaluate the performance of our approach within the unit square and puzzle piece domains by demonstrating its accuracy and CPU-time performance. The governing equation involves updating $f(t, x, u, \nabla u, \Delta u)$ in (4.1), outlined as follows:

$$ \frac{\partial u}{\partial t} = \nu \Delta u - \frac{1}{2} \nabla \cdot u^2. \qquad (5.1) $$

where $\nu = 0.005$ and $t \in [0.25, 1.25]$. The equation is augmented with non-homogeneous Dirichlet boundary conditions and initial solution, which yield an exact
solution expressed as:

$$u(x, y, t) = 1 / \left( 1 + \exp \left( \frac{x + y - t}{2\nu} \right) \right). \qquad (5.2)$$

This example is extensively employed in the literature [29, 35].

**Example 4 (unit-square).** Table 4 shows a performance of adaptive mesh especially in the sense of $H^1$-error. Figure 11 illustrates the concurrent evolution of a linear shock and the adaptive mesh over time, maintaining a consistently high level of accuracy without oscillation within the adaptive mesh framework. Figure 12 shows the error convergence under h-refinement both with and without adaptive mesh, while Figure 13 demonstrates the CPU-time efficiency employing Algorithm 3. The Newton iterations yield results comparable to those obtained with a uniform mesh, all while preserving high accuracy in the approximation. The oscillations are effectively eliminated with just 32 elements in each direction, thanks to the adaptive mesh concentration around the solution interface.

<table><caption>**Table 4: Example 4:** Comparison of different results of computations error in a uniform and adaptive meshes using bi-quadratic B-spline at $T = 1.25$.</caption><thead><tr><th rowspan="2">#cells</th><th colspan="3">Uniform mesh</th><th colspan="3">Adaptive mesh using (4.10)</th></tr><tr><th>L<sup>2</sup>-error</th><th>H<sup>1</sup>-error</th><th>Newton-iter</th><th>L<sup>2</sup>-error</th><th>H<sup>1</sup>-error</th><th>Newton-iter</th></tr></thead><tbody><tr><td>16</td><td>4.34e-02</td><td>3.49e+00</td><td>5</td><td>9.15e-03</td><td>8.41e-01</td><td>5</td></tr><tr><td>32</td><td>9.57e-03</td><td>1.46e+00</td><td>4</td><td>1.07e-03</td><td>1.37e-01</td><td>4</td></tr><tr><td>64</td><td>1.77e-03</td><td>3.22e-01</td><td>3</td><td>3.22e-04</td><td>8.40e-02</td><td>3</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Fig. 11: Example 4:** (top) solution in adapted mesh at $t = 0.38, 0.51, 0.76, 1.25$ (bottom) projection of solution at $t = 1.25$ following the black line, and adaptive mesh at $t = 0.51, 0.7, 1.16$ (#cells=32).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12: Example 4:** (left) $L^2$-error evolution under h-refinement and (right) $H^1$-error evolution under h-refinement at $t = 0.35$. (un for uniform & ad for adaptive mesh)</figcaption></figure>

**Example 5 (puzzle piece).** Table 5 shows the performance of the adaptive mesh on the puzzle piece. Figure 14 shows the same results as in unit-square case, while preserving the exact geometry during r-refinement. Adaptive mesh is influenced by the quality of parameterization. As shown in the adaptive mesh in Figure 14, the intensity of the adaptive meshes varies across the entire domain due to the non-uniformity of the initial mesh.
<figure><img src="image_16.png" /><figcaption>Fig. 13: Example 4: (left) $L^2$-error under h-refinement, and (middle) $H^1$-error under h-refinement, and (right) mesh-size as a function of CPU-time at $t = 0.35$. (un for uniform & ad for adaptive mesh)</figcaption></figure>

<table><caption>Table 5: Example 5: Comparison of different results of computations error in a uniform and adaptive meshes using bi-quadratic B-spline at $T = 1.25$.</caption><thead><tr><th rowspan="2">#cells</th><th colspan="3">Uniform mesh</th><th colspan="3">Adaptive mesh using (4.10)</th></tr><tr><th>L<sup>2</sup>-error</th><th>H<sup>1</sup>-error</th><th>Newton-iter</th><th>L<sup>2</sup>-error</th><th>H<sup>1</sup>-error</th><th>Newton-iter</th></tr></thead><tbody><tr><td>16</td><td>5.34e-02</td><td>3.66e+00</td><td>8</td><td>2.20e-02</td><td>1.87e+00</td><td>8</td></tr><tr><td>32</td><td>1.98e-02</td><td>2.26e+00</td><td>6</td><td>1.56e-03</td><td>3.16e-01</td><td>5</td></tr><tr><td>64</td><td>1.19e-03</td><td>3.19e-01</td><td>3</td><td>3.248e-04</td><td>8.56e-02</td><td>3</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>Fig. 14: Example 5: (top) solution in adaptive mesh at $t = 0.38$, 0.64, 0.9, 1.25 (bottom) projection of solution at $t = 1.25$ following the dashed black line, and adaptive mesh at $t = 0.64$, 0.9, 1.25 (#cells=32).</figcaption></figure>
## 5.2 Allen-Cahn equation

In this subsection, we consider the Allen-Cahn equation [44, 45], introduced as follows:

$$ \frac{\partial u}{\partial t} = \Delta u - \frac{\mathcal{E}'(u)}{\varepsilon_m^2}, \qquad (5.3) $$

where $\mathcal{E}(u) = 0.25(u^2 - 1)^2$ is the potential energy function, $\varepsilon_m$ represents the thickness of the transition layer. By means of curvature flow, we choose $\varepsilon_m$ to be small so that a solution develops two bulk regions with $\pm 1$ min-max value separated by the so-called diffuse interface. We use for the analysis, the energy functional as:

$$ \text{Energy} = \int_{\Omega_p} \frac{\mathcal{E}(u)}{\varepsilon_m^2} + \frac{|\nabla u|^2}{2} dx. \qquad (5.4) $$

The equation (5.3) is formulated with homogeneous Neumann boundary conditions.

**Example 6.** We firstly investigate the parameter $\varepsilon_m$, formulated as:

$$ \varepsilon_m = \frac{m}{2\sqrt{2}\arctan(0.9)(64 + 1)}, \qquad (5.5) $$

where m is the number of cells representing the thickness, so the sharpness of the solution increases as m decreases [44]. To demonstrate the efficiency of the mesh adaptation for small values of m, we start by measuring the circle's radius that changes with time, as defined by:

$$ u(x, y, t) = \tanh \left( \frac{r(t) - \sqrt{(x - 0.5)^2 + (y - 0.5)^2}}{\sqrt{2}\varepsilon_m} \right), \qquad (5.6) $$

where $r(t) = \sqrt{0.25^2 - 2t}$, as calculated in [44]. We maintain #cells= 64 in each direction while varying the value of m.

Table 6 and Figures 15-17 show a good performance of mesh adaptation for a small thickness, displaying how well mesh adaptation works on a sharp moving interface. The energy in Table 6 decreases as m increases due to the fact that the interface diffuses for large values of m, which is clear in the projection of solution in Figures 15-17. Newton iterations remain stable and comparable to the uniform mesh. Otherwise, we are interested in smaller thicknesses m = 2. In Figure 16, at the earlier time T=5.5e-4, it is evident that increasing the degree has an effect on capturing the sharp interface results in more oscillations for a uniform mesh, whereas this effect is not observed in an adaptive mesh because the mesh is well adapted around the interface. Moreover, the L²-error demonstrates superior performance across varying degrees.
<table><caption>**Table 6: Example 6:** Comparison of different results of computations error in a uniform and adaptive meshes using bi-quadratic B-spline at $T = 0.015$.</caption><thead><tr><th rowspan="2">Degree</th><th rowspan="2">m</th><th colspan="3">Uniform mesh</th><th colspan="3">Adaptive mesh using (4.10)</th></tr><tr><th>L<sup>2</sup>-error</th><th>Energy</th><th>Newton-iter</th><th>L<sup>2</sup>-error</th><th>Energy</th><th>Newton-iter</th></tr></thead><tbody><tr><td rowspan="2">p=2</td><td>2</td><td>7.06e-02</td><td>1.50e+02</td><td>4</td><td>3.91e-03</td><td>1.44e+02</td><td>6</td></tr><tr><td>5</td><td>4.17e-03</td><td>5.76e+01</td><td>3</td><td>3.37e-03</td><td>5.76e+01</td><td>6</td></tr><tr><td rowspan="2">p=3</td><td>2</td><td>4.19e-02</td><td>1.48e+02</td><td>4</td><td>2.88e-03</td><td>1.44e+02</td><td>6</td></tr><tr><td>5</td><td>3.50e-03</td><td>5.76e+01</td><td>3</td><td>3.29e-03</td><td>5.76e+01</td><td>5</td></tr><tr><td rowspan="2">p=4</td><td>2</td><td>3.26e-02</td><td>1.47e+02</td><td>4</td><td>1.81e-03</td><td>1.44e+02</td><td>6</td></tr><tr><td>5</td><td>3.42e-03</td><td>5.76e+01</td><td>3</td><td>3.29e-03</td><td>5.76e+01</td><td>5</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Fig. 15: Example 6 m = 2:** (top) Adaptive solution at $t = 0$, 6.7e-3, 1e-2, 1.5e-2 (bottom) A projection of solution at $t = 0.015$ following the dashed black line, and adaptive meshes at $t = 6.7e-3$, 1e-2, 1.5e-2.</figcaption></figure>

**Example 7.** We evaluate our solver with $m = 2$, initialized from various configurations as computed in the paper [45]. Figures 18-19-20 sequentially depict the temporal evolution of a star, a torus, and a dumbbell, in parallel with their respective adaptive meshes.

This subsection is concluded with a random solution. The phase separation begins with an initial guess having an amplitude of $\pm 0.1$. Figures 21-22, using $m = 2$, illustrate the progression of phase separation, reaching an amplitude of $\pm 1$ over time.

Note that we have performed tests on a circular domain with a radius of 0.565 for $m = 2$ in order to maintain an equivalent physical surface size as in the unit square. Otherwise, the numerical solution for a circular domain with a radius of 1 with $m = 2$ falls outside the physical range of $[-1, 1]$.
<figure><img src="image_19.png" /><figcaption>Fig. 16: Example 6: $m=2$, $T = 0.00055$: (left) $L^2$-error evolution for various polynomial degrees $p$ along each direction (right) A projection of solution on $x = 0.5$ using number of elements 32×32.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 17: Example 6 m = 5: (top) Adaptive solution at $t = 0$, 6.7e-3, 1e-2, 1.5e-2 (bottom) A projection of solution at $t = 0.015$ following the dashed black line, and adaptive meshes at $t = 0$, 6.7e-3, 1e-2, 1.5e-2.</figcaption></figure>

The numerical results demonstrate a good assimilation of the adaptive moving mesh with respect to the behavior of the solution. As a result, even small oscillations at the sharp interface of the solution are captured accurately by the adaptive mesh.
<figure><img src="image_21.png" /><figcaption>Fig. 18: Example 7 m = 2: (top) Adaptive solution at $t = 0, 1.6e-3, 1.6e-2, 3.3e-2$ (bottom) A projection of solution at $t = 3.3e-2$ following the dashed black line, and adaptive meshes at $t = 1.6e-3, 1.6e-2, 3.3e-2$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19: Example 7 m = 2: (top) Adaptive solution at $t = 0, 2.7e-2, 3.5e-2, 4.5e-2$ (bottom) A projection of solution at $t = 4.5e-2$ on $x = 0.5$, and adaptive meshes at $t = 2.7e-2, 3.5e-2, 4.5e-2$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 20: Example 7 m = 2: (top) Adaptive solution at $t = 0$, 2e-3, 6.7e-3, 9.8e-3 and (bottom) A projection of solution at $t = 9.8e-3$ following the dashed black line, and adaptive meshes at $t = 2e-3$, 6.7e-3, 9.8e-3.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 21: Example 7 m = 2: (top) Adaptive solution at $t = 0$, 6.7e-4, 1.5e-3, 4.4e-3 (bottom) A projection of solution at $t = 4.4e-3$ following the dashed black line, and adaptive meshes at $t = 6.7e-4$, 1.5e-3, 4.4e-3.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 22: Example 7 m = 2: (top) Adaptive solution at t = 0, 6.7e-4, 1.6e-3, 4.05e-3 (bottom) A projection of solution at t = 4.05e-3 following the dashed black line, and adaptive meshes at t = 6.7e-4, 1.6e-3, 4.05e-3.</figcaption></figure>
## 5.3 Cahn-Hilliard equation

In this subsection, we study the Cahn-Hilliard model, which was studied in the context of IGA by H. Gomez et al. [3]. Let $c$ denote the concentration of one of a given binary mixture components, which we assume to be governed by the Cahn-Hilliard equation. Then the problem is given in strong form:

$$ \frac{\partial c}{\partial t} = \nabla \cdot (M_c \nabla(\mu_c - \lambda \Delta c)), \qquad (5.7) $$

where $M_c$ is the mobility, $\mu_c$ is the chemical potential of a regular solution in the absence of a phase interface, and $\lambda$ is a positive constant such that $\sqrt{\lambda}$ represents a length scale of the problem. We consider a non-linear relationship between the mobility and the concentration, given by

$$ M_c = Dc(1-c), \qquad (5.8) $$

where $D$ is a positive constant. Finally, $\mu_c$ is a highly non-linear function of the concentration, defined as

$$ \mu_c = \frac{1}{3} \log \frac{c}{1-c} + 1 - 2c. \qquad (5.9) $$

In relation to the equation (4.1), we reformulate the Cahn-Hilliard equation in its dimensionless coordinate by taking

$$ f(t^*, x^*, c, \nabla^* c, \Delta^* c) = \nabla^* \cdot (M \nabla^*(\mu - \lambda \Delta^* c)), \qquad (5.10) $$

where we make some changes of variable $t^* = t/(L_0^4/D\lambda)$, $x^* = x/L_0$, which results in

$$ M = c(1-c), $$

$$ \mu = \alpha\mu_c, \quad \alpha = \frac{L_0^2}{3\lambda}. $$

$L_0$ represent a length scale, while $\alpha$ is a dimensionless parameter that exclusively defines a solution.

**Remark 4.** In our numerical tests, as in [3], we consider periodic boundary conditions in all directions.

The initial solution is set as $c(t=0) = \bar{c}+r$, where $r$ represents a random variable with an amplitude of $\pm 0.05$. The time step $\delta t$ remains fixed at 5e-8. Typically, under the tolerance of 1e-4, the solver requires 4-10 GMRES iterations.

**Example 8.** In this example, we tested the Cahn-Hilliard solver with $\#cells = 64$ using $\bar{c} = 0.63$ and $\alpha = 6000$. We chose this example for comparison because, in this case, the thinner interface is more representative from a physical perspective. Figure 23 shows the high quality of mesh adaptation, the adaptive mesh focuses on the interfaces of Cahn-Hilliard solutions, identified through the utilization of the arc-length monitor function (4.10). By adjusting the coefficients in (4.10), the intensity of the adaptive mesh can be increased and a density can be manipulated to cover more cells. However, such adjustments necessitate adapting the time step to avoid stability issues. From the
numerical results shown in Figure 24, we observe that the adaptive mesh accurately
captures the interface of the Cahn-Hilliard solution. This has significant implications
for the behavior of the solution, as shown in Figures 25-26. This can be understood
from the study presented in [3], which required a fine grid with 128 elements in each
direction to compute a topologically correct solution associated with $\bar{c} = 0.63$.

<figure><img src="image_26.png" /><figcaption>Fig. 23: Example 8:(left) Approximate solution in uniform mesh (right) Approximate solution in adaptive mesh and corresponding adaptive mesh, and monitor function.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. 24: Example 8: Projection of solution following the dashed black lines from left to right at $t = 2.3 \times 10^{-4}$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 25: Example 8:(left) Approximate solution in uniform mesh (right) Approximate solution in adaptive mesh and corresponding adaptive mesh, and monitor function.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Fig. 26: Example 8: Projection of solution following the dashed black lines from left to right at $t = 7.5e-4$. (un for uniform & ad for adaptive mesh)</figcaption></figure>
# 6 Conclusion

In conclusion, we have introduced a novel moving mesh approach for isogeometric analysis. By utilizing a computed mapping derived from the parabolic Monge-Ampere equation to accommodate variations in a specified monitor, this method achieves high mesh equidistribution quality while preserving bijection across various cases. Its application in parabolic time-dependent problems shows promise in addressing singularities and capturing sharp interfaces, all while reducing the required degrees of freedom compared to a uniform mesh.

Through this adaptive framework, we have developed a new, fast PMAE solver, which is simple and easy to implement. The monitor function is well-designed to concentrate efforts where simulation is most needed. Moreover, the new scheme for PDEs in a moving mesh is presented and is proven through examples to be stable and efficient. Importantly, the scheme avoids any additional operations such as projection to compute the next solution.

To the best of our knowledge, the present paper is one of pioneering instances in employing dynamic r-refinement -moving mesh- for Isogeometric Analysis.

# 7 Declarations

The authors have no relevant financial or non-financial interests to disclose.

The authors have no competing interests to declare that are relevant to the content of this article.

All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript.

# References

[1] Hughes, T.J., Cottrell, J.A., Bazilevs, Y.: Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 194(39-41), 4135-4195 (2005)

[2] Cottrell, J.A., Hughes, T.J.R., Bazilevs, Y.: Isogeometric Analysis, Toward Integration of CAD and FEA. John Wiley & Sons, (2009)

[3] Gómez, H., Calo, V.M., Bazilevs, Y., Hughes, T.J.R.: Isogeometric analysis of the Cahn-Hilliard phase-field model. Computer Methods in Applied Mechanics and Engineering 197(49), 4333-4352 (2008) https://doi.org/10.1016/j.cma.2008.05.003

[4] Vuong, A.-V., Giannelli, C., Jüttler, B., Simeon, B.: A hierarchical approach to adaptive local refinement in isogeometric analysis. Computer Methods in Applied Mechanics and Engineering 200(49), 3554-3567 (2011) https://doi.org/10.1016/j.cma.2011.09.004
[5] Carraturo, M., Giannelli, C., Reali, A., Vázquez, R.: Suitably graded THB-spline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes. Computer Methods in Applied Mechanics and Engineering **348**, 660–679 (2019) https://doi.org/10.1016/j.cma.2019.01.044

[6] Johannessen, K.A., Kvamsdal, T., Dokken, T.: Isogeometric analysis using LR B-splines. Computer Methods in Applied Mechanics and Engineering **269**, 471–514 (2014) https://doi.org/10.1016/j.cma.2013.09.014

[7] Wei, X., Zhang, Y., Liu, L., Hughes, T.J.R.: Truncated T-splines: Fundamentals and methods. Computer Methods in Applied Mechanics and Engineering **316**, 349–372 (2017) https://doi.org/10.1016/j.cma.2016.07.020. Special Issue on Isogeometric Analysis: Progress and Challenges

[8] Liu, L., Casquero, H., Gomez, H., Zhang, Y.J.: Hybrid-degree weighted T-splines and their application in isogeometric analysis. Computers & Fluids **141**, 42–53 (2016) https://doi.org/10.1016/j.compfluid.2016.03.020. Advances in Fluid-Structure Interaction

[9] Xu, G., Mourrain, B., Duvigneau, R., Galligo, A.: Parameterization of computational domain in isogeometric analysis: Methods and comparison. Computer Methods in Applied Mechanics and Engineering **200**(23), 2021–2031 (2011) https://doi.org/10.1016/j.cma.2011.03.005

[10] Xu, G., Li, B., Shu, L., Chen, L., Xu, J., Khajah, T.: Efficient r-adaptive isogeometric analysis with Winslow’s mapping and monitor function approach. Journal of Computational and Applied Mathematics **351**, 186–197 (2019) https://doi.org/10.1016/j.cam.2018.11.003

[11] Hénap, G., Szabó, L.: On the configurational-force-based r-adaptive mesh refinement in isogeometric analysis. Finite Elements in Analysis and Design **124**, 1–6 (2017) https://doi.org/10.1016/j.finel.2016.10.002

[12] Umesh Basappa, A.R., Reddy, J.N.: Adaptive isogeometric analysis based on a combined r-h strategy. International Journal for Computational Methods in Engineering Science and Mechanics **17**(2), 73–92 (2016) https://doi.org/10.1080/15502287.2016.1153171 https://doi.org/10.1080/15502287.2016.1153171

[13] Ji, Y., Wang, M.-Y., Wang, Y., Zhu, C.-G.: Curvature-based R-Adaptive Planar NURBS Parameterization Method for Isogeometric Analysis Using Bi-Level Approach. Computer-Aided Design **150**, 103305 (2022) https://doi.org/10.1016/j.cad.2022.103305

[14] Habib, S.H., Kezrane, C., Hachi, B.E.: Moving local mesh based on analysis-suitable T-splines and bézier extraction for extended isogeometric finite element analysis - Application to two-dimensional crack propagation. Finite Elements in Analysis and Design **213**, 103854 (2023) https://doi.org/10.1016/j.finel.2022.
[15] Scott, M.A., Li, X., Sederberg, T.W., Hughes, T.J.R.: Local refinement of analysis-suitable T-splines. Computer Methods in Applied Mechanics and Engineering **213-216**(0), 206-222 (2012) https://doi.org/10.1016/j.cma.2011.11.022

[16] Hennig, P., Ambati, M., De Lorenzis, L., Kästner, M.: Projection and transfer operators in adaptive isogeometric analysis with hierarchical B-splines. Computer Methods in Applied Mechanics and Engineering **334**, 313-336 (2018) https://doi.org/10.1016/j.cma.2018.01.017

[17] Kumar, M., Kvamsdal, T., Johannessen, K.A.: Simple a posteriori error estimators in adaptive isogeometric analysis. Computers & Mathematics with Applications **70**(7), 1555-1582 (2015) https://doi.org/10.1016/j.camwa.2015.05.031. High-Order Finite Element and Isogeometric Methods

[18] Buffa, A., Giannelli, C.: Adaptive isogeometric methods with hierarchical splines: Error estimator and convergence. Mathematical Models and Methods in Applied Sciences **26**(01), 1-25 (2016) https://doi.org/10.1142/S0218202516500019 https://doi.org/10.1142/S0218202516500019

[19] Kumar, M., Kvamsdal, T., Johannessen, K.A.: Superconvergent patch recovery and a posteriori error estimation technique in adaptive isogeometric analysis. Computer Methods in Applied Mechanics and Engineering **316**, 1086-1156 (2017) https://doi.org/10.1016/j.cma.2016.11.014. Special Issue on Isogeometric Analysis: Progress and Challenges

[20] Ghorashi, S.S., Valizadeh, N., Mohammadi, S., Rabczuk, T.: T-spline based XIGA for fracture analysis of orthotropic media. Computers & Structures **147**, 138-146 (2015) https://doi.org/10.1016/j.compstruc.2014.09.017. CIVIL-COMP

[21] Bazilevs, Y., Takizawa, K., Tezduyar, T.E., Korobenko, A., Kuraishi, T., Otoguro, Y.: Computational aerodynamics with isogeometric analysis. Journal of Mechanics **39**, 24-39 (2023) https://doi.org/10.1093/jom/ufad002 https://academic.oup.com/jom/article-pdf/doi/10.1093/jom/ufad002/49478788/ufad002.pdf

[22] Takizawa, K., Tezduyar, T.E., Avsar, R.: A low-distortion mesh moving method based on fiber-reinforced hyperelasticity and optimized zero-stress state. Computational Mechanics **65**, 1567-1591 (2020) https://doi.org/10.1007/s00466-020-01835-z

[23] Xie, X., Yang, A., Wang, Y., Jiang, N., Wang, S.: Fully adaptive isogeometric topology optimization using MMC based on truncated hierarchical B-splines. Structural and Multidisciplinary Optimization **63**, 2869-2887 (2021) https://doi.org/10.1007/s00158-021-02850-1
[24] Maestre, J., Pallares, J., Cuesta, I., Scott, M.A.: A 3d isogeometric BE-FE analysis with dynamic remeshing for the simulation of a deformable particle in shear flows. Computer Methods in Applied Mechanics and Engineering **326**, 70–101 (2017) https://doi.org/10.1016/j.cma.2017.08.003

[25] Bahari, M., Habbal, A., Ratnani, A., Sonnendrücker, E.: Adaptive Isogeometric Analysis using optimal transport and their fast solvers. Computer Methods in Applied Mechanics and Engineering **418**, 116570 (2024) https://doi.org/10.1016/j.cma.2023.116570

[26] Lynch, R.E., Rice, J.R., Thomas, D.H.: Direct solution of partial difference equations by tensor product methods. Numerische Mathematik **6**(1), 185–199 (1964)

[27] Sangalli, G., Tani, M.: Isogeometric preconditioners based on fast solvers for the Sylvester equation. SIAM Journal on Scientific Computing **38**(6), 3644–3671 (2016)

[28] Budd, C.J., Williams, J.F.: Parabolic Monge–Ampère methods for blow-up problems in several spatial dimensions. Journal of Physics A: Mathematical and General **39**(19), 5425 (2006) https://doi.org/10.1088/0305-4470/39/19/S06

[29] Budd, C.J., Williams, J.F.: Moving Mesh Generation Using the Parabolic Monge–Ampère Equation. SIAM Journal on Scientific Computing **31**(5), 3438–3465 (2009) https://doi.org/10.1137/080716773 https://doi.org/10.1137/080716773

[30] N., G.R., T., L.: Knot Insertion and Deletion Algorithms for B-Spline Curves and Surfaces. SIAM, Philadelphia, USA (1993)

[31] R., J.: Review of A Practical Guide to Splines. Mathematics of Computation **34**(149), 325–326 (1980). Accessed 2024-05-15

[32] Buffa, A., Sangalli, G., Vázquez, R.: Isogeometric analysis in electromagnetics: B-splines approximation. Computer Methods in Applied Mechanics and Engineering **199**(17-20), 1143–1152 (2010) https://doi.org/10.1016/j.cma.2009.12.002

[33] Buffa, A., Rivas, J., Sangalli, G., Vázquez, R.: Isogeometric Discrete Differential Forms in Three Dimensions. SIAM J. Numerical Analysis **49**(2), 818–844 (2011)

[34] Sulman, M.M., Williams, J.F., Russell, R.D.: An efficient approach for the numerical solution of the monge–ampère equation. Applied Numerical Mathematics **61**(3), 298–307 (2011) https://doi.org/10.1016/j.apnum.2010.10.006

[35] Huang, W., Russell, R.D.: Adaptive Moving Mesh Methods. Applied mathematical sciences. Springer, New York, Heidelberg, London (2011)
[36] Huang, W., Russell, R.D.: Analysis of Moving Mesh Partial Differential Equations with Spatial Smoothing. SIAM Journal on Numerical Analysis 34(3), 1106–1126 (1997) https://doi.org/10.1137/S0036142993256441
https://doi.org/10.1137/S0036142993256441

[37] Bourne, E., Güçlü, Y., Hadjout, S., Ratnani, A.: Pyccel: a Python-to-X translator for scientific high-performance computing. Journal of Open Source Software 8(83), 4991 (2023) https://doi.org/10.21105/joss.04991

[38] Huang, W., Ren, Y., Russell, R.D.: Moving Mesh Partial Differential Equations (MMPDES) Based on the Equidistribution Principle. SIAM Journal on Numerical Analysis 31(3), 709–730 (1994). Accessed 2023-11-18

[39] Budd, C.J., Huang, W., Russell, R.D.: Moving Mesh Methods for Problems with Blow-Up. SIAM Journal on Scientific Computing 17(2), 305–327 (1996) https://doi.org/10.1137/S1064827594272025
https://doi.org/10.1137/S1064827594272025

[40] Cao, W., Huang, W., Russell, R.D.: Comparison of two-dimensional r-adaptive finite element methods using various error indicators. Mathematics and Computers in Simulation 56(2), 127–143 (2001) https://doi.org/10.1016/S0378-4754(01)00285-3. Method of lines

[41] Li, R., Tang, T., Zhang, P.: A moving mesh finite element algorithm for singular problems in two and three space dimensions. Journal of Computational Physics 177(2), 365–393 (2002)

[42] Tang, T.: Moving mesh methods for computational fluid dynamics. Contemporary Mathematics 383(8), 141–173 (2005)

[43] Jones, E., Oliphant, T., Peterson, P., et al.: Scipy: Open source scientific tools for Python (2001)

[44] Li, Y., Lee, H.G., Jeong, D., Kim, J.: An unconditionally stable hybrid numerical method for solving the Allen–Cahn equation. Computers & Mathematics with Applications 60(6), 1591–1606 (2010) https://doi.org/10.1016/j.camwa.2010.06.041

[45] Kim, Y., Ryu, G., Choi, Y.: Fast and accurate numerical solution of Allen–Cahn equation. Mathematical Problems in Engineering 2021, 1–12 (2021) https://doi.org/10.1155/2021/5263989