A Velocity-based Moving Mesh Virtual Element Method

We present a velocity-based moving mesh virtual element method for the numerical solution of PDEs involving moving boundaries. The virtual element method is used for computing both the mesh velocity and a conservative Arbitrary Lagrangian-Eulerian solution transfer on general polygonal meshes. The approach extends the linear finite element method to polygonal mesh structures, achieving the same degree of accuracy. In the context of moving meshes, a major advantage of the virtual element approach is the ease with which nodes can be inserted on mesh edges. Demonstrations of node insertion techniques are presented to show that moving polygonal meshes can be simply adapted for situations where a boundary encounters a solid object or another moving boundary, without reduction in degree of accuracy.


Introduction
We propose a new numerical method for the simulation of free boundary problems over general polygonal meshes based on combining a moving mesh algorithm with the Virtual Element Method (VEM). The new approach is presented focusing on the solution of the free boundary problem for the porous medium equation as a model problem.
Moving mesh methods form part of a large class of adaptive mesh refinement techniques alongside h and p refinement strategies. The primary advantage of moving mesh methods is the ability to optimize mesh structures without requiring any change in the mesh connectivity, which can cause computational challenges, particularly when parallel implementations are considered. They also provide a natural framework for tracking physical features of a time-dependent PDE, such as blow-up problems [21], phase change modelling [10], fluid-structure interaction problems [45], and more general time-dependent PDEs [5,41,31].
Here we consider a velocity-based moving mesh algorithm [6,8,9,35,38,7] for the numerical solution of free boundary problems. The method is closely related to the Geometric Conservation Law method (GCL) [25] and also forms part of a larger family of adaptive moving mesh methods [33]. It uses a Lagrangian formulation of a given PDE to solve directly for the mesh velocities which are integrated over time to evolve the mesh. The solution computed on a previous time-step is then transferred to a new time-step using an Arbitrary Lagrangian Eulerian (ALE) scheme [27] based on a weak distribution of a given monitor function.
We combine the velocity-based moving mesh algorithm with the virtual element method for spatial discretisation on general polygonal meshes [13].
The VEM was first introduced in [13] for the conforming discretisation of linear elliptic problems; cf. also [1,14,23,46] for extensions and implementation details. It provides a flexible discretization framework for the design of compatible general mesh methods that incorporate, at the discrete level, fundamental properties of the continuous problem at hand, such as topology, conservation, symmetry, and positivity. To achieve compatibility, discrete virtual element spaces are typically defined on individual polygonal elements implicitly through local boundary value problems. These are understood through a set of degrees of freedom instead of explicit basis functions. The use of local polynomial projections of the virtual test functions permits the definition of virtual element weak formulations which are computable at a cost on par with that of standard Finite Element Methods (FEMs). Here we consider the lowest-order VEM, which can be seen as a generalisation of the standard linear FEM to general polygonal/polyhedral meshes. Just like linear FEM, the lowest-order VEM can be set up using only vertex information and in particular without the use of computationally expensive quadrature for nodes on polygons. As such, our approach can be interpreted as a generalisation of the linear finite element moving mesh method [7]. Numerical results demonstrate that optimal orders of convergence are achieved when applying a virtual element discretization on a diverse range of polygonal mesh structures.
To simplify the presentation, the novel moving mesh VEM is developed for the two-dimensional porous medium equation (PME). The PME is a parabolic non-linear diffusion equation [51]. It is an ideal benchmark for a moving mesh method as it provides time-varying solutions with compact support and a moving boundary. However, the moving mesh VEM proposed here may be immediately generalised to other free boundary problems. Generalisations of the velocitybased method to other PDEs can be found in the literature [6,7] and the VEM has already been generalised to the solution of a wide range of PDEs [12,50,48], including flow in porous media [49,56,52].
In recent years, several works have been published which involve time-dependent and moving polygonal domains. These developments include: ALE maps between two meshes using virtual element and discontinuous Galerkin (dG) methods [37], a cut-based dG scheme for fluid-structure interaction problems [2], an ALE scheme on time-dependent reconstructed Voronoi meshes [30], adaptive methods using elliptic reconstruction techniques on moving meshes [22], polygonal mesh quality measures for an anisotropic MMPDE scheme [34], and a VEM for long term simulations of moving landforms [40]. The benefits of using polygonal discretizations for moving mesh methods include representing moving boundaries and interfaces with a minimal number of degrees of freedom and localized mesh refinement when a change in mesh connectivity is required, such as in contact problems. To show-case its potential in this respect, we test the moving mesh VEM for the numerical treatment of both obstacle and self-intersecting moving interface problems.
The layout of the remaining sections of this paper is as follows. In Section 2 we introduce the moving mesh algorithm, including the weak formulations necessary for the method. The components of the VEM discretisation are split into two sections. Section 3 reviews the fundamentals of the VEM and formulates a method for computing mesh velocities at a fixed point in time. Section 4 presents a framework for moving virtual elements along with a VEM for a conservative ALE update scheme. An overview of the algorithm and implementation details is given in Section 5. Node insertion/removal algorithms are discussed in Section 6. A set of numerical experiments is presented in Section 7 for the PME before drawing some conclusions in Section 8.

The Moving Mesh Framework
In this section we present the key components of the velocity-based moving mesh framework for a general evolutionary PDE of the form with t ∈ (0, T ], T ∈ R + , indicating time and L a generic spatial differential operator in R d . The choice of spatial operator will be made specific later on when we introduce the porous medium equation as the model problem. Further details on a finite element discretisation of this moving mesh method for general parabolic problems can be found in [7] and the references therein. The time-dependent support of the evolving solution to Equation (1) will be denoted by Ω t ⊆ R d , with boundary ∂Ω t which is allowed to move. The temporal superscript t will be swapped for n when considering discrete time levels t n or, when no confusion arises, the superscript notation will be dropped entirely for ease of reading. A moving coordinate system is defined by x = x(ξ, t) where x(ξ, 0) = ξ is the coordinate system on the initial domain Ω 0 , which is used as the reference domain for simplicity. It is assumed that, for all t ∈ [0, T ], the map x(·, t) is Lipschitz with Lipschitz inverse. The evolution of Ω t is determined by the velocity field v = ∂x ∂t . (2) Following [16], we shall also refer to the space-time domain

The Velocity Field
The movement of the mesh is derived by specifying the time evolution of the distribution of the spatial integral of some solution-dependent monitor function, M(ρ). Specifically, the velocity field v is determined by requiring that the initial distribution of the monitor function M(ρ) is conserved as time progresses. Namely, we look for v = ∂x ∂t such that, for all w ∈ L 2 (Q T ), the coordinates x(ξ, t) satisfy The derivation of the moving mesh method is based on application of the Reynolds transport theorem and on the assumption that the material derivatives of the weighting functions w(x, t) are zero with respect to the velocity field v (see Equation (2)). Namely, This is a common assumption made in finite element approaches to moving mesh algorithms [37,16,15] and is equivalent to assuming that w(x(ξ, t), t) = w(ξ, 0) in Equation (3). Equidistribution-based mesh movement algorithms typically attempt to reduce the global approximation error without changing the number of degrees of freedom by choosing a finite set of weighting functions w i (x, t), so that each one has a direct association with a mesh node or element. The monitor function M(ρ) is then selected to act as a local error indicator, and the mesh is adjusted in an attempt to equidistribute the values of the weighted monitor integrals and hence to equidistribute the local error across the mesh. Such an approach could be followed within our framework [38,7]. However, in this paper we adopt the following approach which is more akin to Lagrangian mesh movement algorithms, which attempt to move the mesh with the 'flow' velocity.
Choosing M(ρ) = ρ in equation (5) naturally leads to a weak approximation of the Lagrangian 'flow' velocity of the PDE when Lρ = ∇ · f for some flux f in Equation (1). This has two benefits: (a) it allows us to predict the movement of free boundaries; (b) it reduces interpolation error of the computed mesh and solution between time-steps because the mesh (and hence the solution) is transported with the velocity field inherent to the PDE.
In many PDEs (including the PME), ρ represents a density, i.e. its integral in space is a mass, so we will refer to M(ρ) = ρ as the mass monitor. For clarity of presentation, we will assume use of this mass monitor from now on. Hence, each weight function w(x, t) is assigned its own 'mass': Equation (3) then provides c(w) = µ 0 (w)/θ 0 from the initial conditions of the PDE, and we assume that these distribution coefficients remain constant in time. The evolution of ρ is governed by Equation (1), so Equation (6) provides us with a way to prescribe the evolution of the coordinate system in a way which retains the initial 'mass' distribution.
Remark 1. The original moving mesh finite element method [6,7] chose the weight functions to be the standard linear Lagrange finite element test functions on meshes of simplices. This associates a fixed proportion of the total mass of the system with each node of the mesh, so the values of µ t (w) depend not only on ρ, but also on the mesh node positions. In this situation, Equation (6) provides a way to compute mesh node velocities using standard finite element techniques, in a way which is consistent with local mass conservation when this is a property of the underlying PDE. Our results to follow demonstrate that the same principle can be applied on polygonal meshes within a virtual element framework.
Let t ∈ [0, T ] and consider all w ∈ H 1 (Q T ) such that D t w(t) = 0, cf. Equation (4). Differentiating the first equality in (6) with respect to time and applying the Reynolds transport theorem [53] leads tȯ Noting that (7) does not fully determine the velocity, we further require that each weighting function retains a fixed proportion of θ t as the mesh evolves. Hence we impose the constraintsμ Inserting (8) in (7), using (1), and integrating by parts results in with n the outward unit normal on ∂Ω t . Similarly, a direct application of the Reynolds transport theorem to the second equality in (6), yieldṡ In fact,θ t is typically known explicitly because ρv · n is known on the whole of ∂Ω t , and the boundary conditions provided with the PDE will enable evaluation of the integral of Lρ. Furthermore, for mass-conservative problems,θ t = 0. Note also that assuming that c(w) remains constant preserves the initial distribution of the mass (or a more general monitor integral) so a monitor which is initially equidistributed between a set of weighting functions should remain equidistributed as the coordinate system evolves with this velocity field.
Equations (9) and (10) are used to compute an instantaneous velocity v which is consistent with conserving the proportion of mass associated with each weighting function. This provides a form of local mass conservation whenθ t = 0. However, it still does not uniquely define v in multiple space dimensions. To overcome this issue, the velocity field is written in terms of its Helmholtz decomposition where φ is a scalar potential and q must be specified. This constraint is equivalent to imposing the curl of the velocity field because ∇ × q = ∇ × v. Moreover, for simplicity, we may further assume that q = 0. This is the natural choice for the porous medium equation, which is derived under the assumption of a curl-free flow velocity field. An example of the method applied with a rotational velocity field is presented in [7]. The problem for determining the velocity potential is therefore: find φ ∈ H 1 (Ω t ) such that where φ = 0 is specified at one point in Ω t to ensure uniqueness. The velocity field is finally obtained as the solution of equation (11) written in weak form with q = 0 and φ given by (12). That is, we find with v · n imposed on any part of the boundary where it is known. For instance, in case of contact with an obstacle, we would impose v · n = 0 on the contact boundary. Note that (13) is just the component-wise L 2 -projection.
Remark 2. The velocity field in the interior of Ω t could be discarded at this stage, and replaced by one derived using a different approach which is constrained by the boundary velocity derived above. For example, interior mesh nodes could be moved using our approach with a different monitor function or a Laplacian smoothing could be applied [38,45].

Solution Update
The velocity field v derived using Equations (12) and (13) can now be used in the update of the solution. Integration of Equation (7) by parts and substitution of Equation (1) results iṅ This is a standard conservative ALE update. Along with Equation (2), it gives a system of ODEs governing the evolution of the coordinate system and the distribution of the mass monitor which can be approximated using standard solvers such as explicit Runge-Kutta methods [7,35]. With these at hand, the solution can be recovered by solving the problem: find ρ ∈ H 1 (Ω t ) such that In the special case whereθ = 0 and c(w) is assumed constant in time, Equation (14) is redundant becauseμ t (w) ≈ 0 by design. In fact, without the velocity recovery from the potential through (12) and (13), we would haveμ t (w) ≡ 0. However, the practical steps of the recovery procedure may introduce small perturbations. Alternatively, one may use the knowledge thatμ t (w) ≡ 0 and hence resort directly to the initial values µ 0 (w), avoiding the need to calculate the ALE update (14) altogether: we refer to this as direct recovery. However, direct recovery can only be used with the specific choice of the mass monitor and with mass-conservative PDEs. In other situations the interior velocity field will not correspond toμ t (w) = 0 and the ALE update is essential.
Remark 3. An alternative, non-conservative, ALE formulation can be derived which would give an equation forρ instead ofμ (see [38]). This derivation does not require that the weight functions evolve with zero material derivative, but conservation of mass is lost as a result.

Porous Medium Equation
On an open, bounded domainΩ ⊂ R d , we consider the following initialboundary value problem for the Porous Medium Equation (PME): find ρ : where Φ = ρ m+1 /(m + 1), for some m > 0, and ρ 0 ≥ 0 having compact support inΩ. The PME belongs to the broader class of Generalized Porous Medium Equations (GPME), also known as filtration equations, obtained with Φ : R + → R + any increasing function. The mathematical analysis of the GPME is well developed; see the monograph [51] and the references therein. In particular, the notion of appropriate weak solutions is discussed in [51] where it is shown that, for non-negative ρ 0 ∈ L 1 (Ω) and for Ψ( there exists a unique non-negative weak solution to the GPME globally in time.
The PME models a number of physical processes such as fluid flow, heat transfer, and diffusion. It exhibits several interesting properties, including the existence of a family of radially symmetric similarity solutions [42], which are used to test the numerical method in Section 7. Other properties of the PME are discussed in [51,7,43]. It is a fundamental example of a degenerate parabolic equation, stemming from the condition that Ψ is non-negative, rather than simply positive. Solutions that exhibit an evolving compact support are ideal for the class of moving mesh methods considered herein because considering as unknown the support of the solution leads to a moving boundary problem that can be simulated over time without having to discretize the entire geometry ofΩ.
Introducing the time-dependent support of ρ as Ω t , we define the time-dependent coordinate system x(ξ, t) with a velocity fieldẋ = v(x, t) that corresponds, at ∂Ω t , with the movement of the free boundary of Ω t for all t ∈ [0, ∞].
Additionally, we consider the free boundary problem for the PME in which part of the boundary may be obstructed by a fixed object. To this end, the boundary is divided into a moving part ∂Ω t M and a fixed part ∂Ω t F , such that Thus, we arrive to the following classical free boundary problem for the PME, whose smooth solutions are weak solutions of (16)-(18), cf. [51].
Here, n is the outward pointing unit normal to the boundary ∂Ω t . Note that an additional (kinematic) boundary condition, which imposes zero flux through the moving boundary, is required to determine the boundary velocity v. On ∂Ω F , the boundary velocity is specified and ρv · n = 0.
Given a transformed coordinate system x and a distribution µ of the variable ρ, our method for solving the PME Problem 2.1 requires the solution of the following problems successively. (15)).

Problem 2.2 (Solution Reconstruction -Equation
where Problem 2.3 (Velocity Potential -Equation (12)). where This has been derived by substituting Lρ = ∇ · (ρ m ∇ρ) in (12) and applying the boundary conditions from Problem 2.1. To reconstruct the velocity field we introduce the follwoing space in which we have the following reconstruction problem. (13)).

Problem 2.4 (Velocity Reconstruction -Equation
where To compensate for the fact that the velocity is reconstructed from the potential (see the comments in Section 2.2), an a posteriorily computed ALE update may be used [27]. Specifically, substituting Lρ = ∇ · (ρ m ∇ρ) and the PME boundary conditions into (14) we arrive to the following. (14)).

Problem 2.5 (ALE Update -Equation
Once again we note that, whenμ t (w) = 0 is assumed, Equation (28) is redundant and Equation (21) is equivalent to a weak form of the PDE where ∇ × v = 0 has been assumed. In other words, the formulation recovers the Lagrangian flow field associated with the porous medium equation.
Equations (28) and (25) can now be combined to give a system of ODEs which determine the evolution of the coordinate system and the monitor distribution. These can be discretized in time using an appropriate solver for first-order ODEs. Specifically, the forward Euler method is used in Section 7.
We now proceed to describe how Problems 2.2-2.5 can be approximated in space on two-dimensional polygonal meshes using the virtual element method.

Virtual Element Method for the Velocity
In this section we introduce the VEM and apply it for the discretisation of Problems 2.3 and 2.4 to determine the domain velocity field. We follow the framework introduced in [13,1,24] for the discretisation of linear elliptic problems. For simplicity, we consider the lowest-order (linear) VEM in two dimensions. Hence, for the rest of the paper, it is assumed that d = 2. In the following we let 0 = t 0 < t 1 < . . . t Nt = T denote a given sequence of discrete time levels at which the PME is discretised.

A Moving Polygonal Mesh
The polygonal representation of the moving domain Ω t at the discrete time level t n is denoted by Ω n h and the corresponding polygonal mesh partitioning Ω n h is denoted by T n h . Polygonal elements and edges within the mesh are denoted by E and e, respectively. For a given polygonal element E, |E| denotes the area, h E the diameter, and (x c , y c ) the centroid.
The following conditions are sufficient to guarantee convergence of the VEM and are assumed to hold true at each time level t n , n = 0, 1, . . . , N t ; see, for example, [13]. 1 A polygon is simple if its boundary forms a closed graph with no intersecting edges other than at each vertex which intersects exactly two edges. Assumption 2 (Shape Regularity). Every E ∈ T n h is a star-shaped domain or a finite union of star shaped domains with respect to a ball of radius greater than γh E for some uniform γ > 0. Additionally, for all edges e ∈ ∂E, the length of e is greater than δh E for some uniform δ > 0.
Remark 4. The shape regularity assumption on the edges required by Assumption 2 can be entirely removed following, for example, [17]. Small edges are relevant, for instance, for the various types of contact problems considered below.

The Π ∇ and Π 0 Projections
Given a subset ω ⊂ R d , d = 1, 2, representing an element or an edge, with P k (ω) we denote the space of polynomials of degree k on ω. Furthermore, we use (·, ·) ω to denote the L 2 inner product over ω.
Critical to the construction of a VEM is the judicious choice of discrete spaces and degrees of freedom that permit the computation of local polynomial projections of virtual functions using only the degrees of freedom. Specifically, the definition of the VEM space to be introduced below is based on an H 1 -type projection operator Π ∇ onto P 1 (E) given, for any w ∈ H 1 (E), by Further, the VEM requires the availability of the L 2 -projection operator Π onto Additionally, the VEM requires the L 2 -projection onto constants of ∇w which we denote by Π 0 , thus The same operator applied component-wise to vector-valued functions will be denoted by Π.

Virtual Element Spaces
Virtual element discrete functions are defined implicitly as solutions of elementwise boundary value problems. These functions are only accessed through their degrees of freedom with which specific projections can be computed; in this way, no evaluation of the (implicitly known) discrete functions is required. Here we consider the linear virtual element space of [1] which is a modification of the original space proposed in [13] for which both the H 1 and L 2 projections introduced in the previous section are computable from the degrees of freedom.
Given a polygonal element E ∈ T n h , we first define the elemental boundary space as Then the local virtual element space of [13] is defined as Note that P 1 (E) ⊆ W (E). In fact, it is immediate to check that the space W (E) corresponds to the elemental linear finite element space whenever E is a triangle, that is W (E) ≡ P 1 (E). In this respect, the VEM can be seen as a generalisation of the linear FEM to polygonal meshes. It is clear from the definition that the only degrees of freedom in the choice of a function w h ∈ W (E) are related to its values on the boundary and, given that w h is linear on each edge, ultimately they depend only on the nodal values; a rigorous proof that these constitute a unisolvent set of degrees of freedom for W (E) is presented in [13]. Hence, in analogy to classical continuous linear finite element spaces, we identify the VEM functions by their nodal values. Additionally, the nodal values allow for the computation of the local Π ∇ projection of any w h ∈ W (E). Indeed, in order to compute the first equation in (30) we need to be able to evaluate, for any p ∈ P 1 (E), the right-hand side term: as ∆p = 0. This shows that such terms only depend on the value of w h on the boundary of E, which are known and are determined by the nodal values. The same is clearly true for the second equation in (30), hence the local Π ∇ projection is fully computable just by accessing the degrees of freedom. However, the L 2 projection is not computable for this space. Hence, following [1], we first consider the augmented spacẽ from which a new enhanced local virtual element space can be defined as Notice that the space V (E) is characterised by the following four properties: (i) we still have that P 1 (E) ⊆ V (E); (ii) the dimension of V (E) is equal to that of W (E) and we can still use the nodal values as degrees of freedom (see [1]); (iii) the local Π ∇ projection is still computable; (iv) if w h ∈ V (E) then Πw h = Π ∇ w h by construction, hence also the L 2 -projection Π is computable. The proof that the gradient projection into constants Π 0 ∇w h is also computable is similar. The global virtual element space is then defined for a given time t n as The dimension of this space is equal to the number of nodes in the mesh, which we shall denote by N dof .
Remark 5. When required, homogeneous Dirichlet boundary conditions can be embedded in the virtual element space by fixing the relevant boundary nodes. Hence, given a Γ ⊆ ∂Ω n h , we denote the constrained space by V n h,Γ = {w h ∈ V n h : w h | Γ = 0}.

Discretisation of the velocity problems
Having defined the virtual element space in Section 3.3, we are ready to present the VEM for the solution of the velocity problems 2.3 and 2.4. The VEM is based on the construction of approximate weak forms which are computable through the elemental projection operators. In particular, the VEM bilinear forms typically involve two terms. The polynomial consistency term acts on the projection of the discrete functions and is responsible for the accuracy of the method; the stabilization term is complementary to the consistency term and is required to ensure the coercivity of the VEM bilinear forms. See, for example, [13,1,23] for more details including proofs of stability and convergence.
Assume that a discrete solution at time level t n has been computed as ρ h ∈ V n h on the current mesh T n h . When n = 0, ρ h is defined as the virtual element interpolation of the initial condition of the PME Problem 2.1. Otherwise, ρ h is the current time discrete solution.
The VEM discretisation of the velocity potential Problem 2.3 reads: given with the approximate bilinear form A h and linear form d h built by summing element-wise contributions as typical of FEM, hence As anticipated above however, the preceding definition of the local VEM forms necessitates the use of projections as follows: where (Πρ h ) 0 is the constant component of Πρ h . Here, following [13], the stabilization form S E A (·, ·) is defined by with m E denoting the dimension of V (E), which for the space considered here is equal to the number of vertices of E, and with dof i (w) representing the i-th degree of freedom of the function w. Hence, dof l (v h ) = v h (x l ) with x l denoting the l-th vertex of E. The integration constant is fixed by constraining a single vertex value of φ h to zero. A variety of suitable stabilization choices are admissible [12,24] but we adopt the simplest choice in this paper. In particular, in the case when arbitrarily small edges appear in the mesh we refer to [18] for more appropriate stabilization terms. The velocity reconstruction Problem 2.4 is a global L 2 projection. Its VEM discretisation reads: given φ h ∈ V n h , the solution of (38), find v h ∈ [V n h ] 2 such that v h · n = 0 on the portion of Ω n h approximating ∂Ω n F and As before, the forms M h and b h are obtained summing the respective elemental forms with the stabilization term S E M (·, ·) given by [1] S

Moving the mesh
The mesh is transferred between discrete time levels by displacing the nodes of the mesh and maintaining the mesh connectivity between t = t n and t = t n+1 . For each mesh node x n the new position is obtained by the forward Euler method applied toẋ = v h (x), yielding x n+1 = x n + (t n+1 − t n )v h (x n ). Thus, consistent with the VEM philosophy, only the values of v h at the nodes, that is the degrees of freedom of v h , are required to compute the mesh movement.

Virtual Element Method for the Solution
Once the new mesh node positions have been computed, we consider the process of updating the solution. This is performed in two steps corresponding, respectively, to the conservative ALE update of the mass monitor Problem 2.5 and the actual solution update Problem 2.2. Before presenting the details on their VEM discretisation, a discussion on the hypothesis leading to such problems is in order. The original moving mesh method in [7] was based on the linear FEM for which the validity at the discrete level of the material derivatives assumption (4) leading to the ALE update (7) has been proven in [36]. In the VEM setting, instead, we exploit the fact that virtual element functions are only accessed through their nodal values: in close alignment with the work presented in [37], only the mesh skeleton velocity is known and used. Hence, in view of the solution update through the time step [t n , t n+1 ], we can assume that (4) is satisfied by the spacetime discrete basis which are then interpolated at the new time level, again, just by accessing the nodal values of the solution.

Discretisation of the solution problems
The initial condition ρ 0 h is approximated by interpolating the degrees of freedom of ρ 0 into the VEM space V 0 h . Then, the initial mass monitor distribution is computed via The next task is the update of the mass monitor over time levels. Once again, this is performed via the forward Euler method: the new monitor is thus given by µ n+1 h (w n+1 h ) = µ n h (w n h ) + (t n+1 − t n )μ n h (w n h ). This, in turn, requires the approximation of the ALE equation (28) which is performed still on the old time level (superscript omitted) bẏ Once the monitor is updated, we set the time level to the new time and update the solution using a VEM discretisation of Problem 2.2: find ρ n+1 Similarly to the velocity reconstruction, the discrete form m h (·, ·) is computed on the new time level by summing over element contributions where with the stabilization term S E m (·, ·) being given by

Partition of unity and conservation
A beneficial property of the linear virtual element method is that the basis functions form a partition of unity on Ω n h at any discrete time level, i.e.
where the set {ϕ i } N dof i=1 refers to the set of canonical VEM basis functions associated to the vertices of the mesh [13]; that is ϕ i (x j ) = δ ij for i = 1, ..., N dof where x j is the j-th node in the mesh.
For a given ρ h ∈ W n h the monitor integral θ n reads from which the polynomial consistency and partition of unity property of the VEM gives Therefore the global conservation of the mass monitor is only dependent on the polynomial component of the discrete solution and weighting functions. Further, the exact value of the monitor can also be recovered via Finally, considering the partition of unity property, the ALE update equation (14) and equation (57) for the PME, we geṫ which agrees with the conservation of mass principle for this particular PDE. In fact, virtual elements preserving relevant global conservation laws in different contexts can be constructed, an example of which is given in [50] for the heat equation.

Implementation Details
This section presents a complete overview of the moving mesh virtual element method. The construction of the required algebraic equations and imposition of boundary conditions are reviewed along with some practical remarks regarding the implementation of this method. The initial weak distribution of the monitor is stored in the vector µ 0 and can be computed using equation (46) whilst the mass matrix M n is computed by assembling the contributions from equation (50). In order to compute the discrete potential φ h ∈ V n h from equation (38), we solve the linear system A n φ = d n with A n and d n computed using equations (39) and (40), respectively. The solution of the resulting linear system determines φ h up to an additive constant which is inherited from the continuous formulation of the method; here we impose dof 1 (φ h ) = 0. Once φ h is recovered, the velocity field is reconstructed solving M n R v = b n , with M n R and b n given by equations (43) and (44), respectively.
The ALE update of vectorμ n is then obtained using equation (47) and, along with the mesh nodal velocitiesẋ n i = v h (ϕ i ), provides a system of ODEs which can be approximated using the forward Euler method. The main method is outlined in algorithm 1.
If required, we can strongly impose any Dirichlet boundary conditions on the solution over time whilst simultaneously conserving the monitor integral, the test space used in recovering the solution is augmented to preserve the partition of unity property (52). In the solution updates, Dirichlet boundary conditions can be enforced using the methodology presented in [35] which extends to a near identical virtual element approach when considering polygonal elements. This is not a necessary step to attain the accuracy results given in section 7 and so the implementation details are not discussed in this paper.

A Contact Algorithm
In this section we discuss two occurrences of contact and present corresponding basic node insertion algorithms that allow for localised and minimal changes to Algorithm 1: Moving mesh VEM input : The initial condition ρ 0 h ∈ W 0 h and mesh T 0 h , the final time T . Set n = 0; Compute µ 0 according to equation (46); while t n < T do Construct and solve A n φ = d n using equations (39) and (40) for φ h ∈ W n h ; Reconstruct the velocity via M n R v = b n ; Compute the ALE updateμ n from equation (47); Select ∆t and set t n+1 = t n + ∆t; Update the mesh node by x n+1 = x n + ∆tv; Update the monitor distribution by µ n+1 = µ n + ∆tμ n ; Figure 1: A demonstration of the mesh self-intersection problem. The initial condition has disconnected support inΩ (left). The disconnected ∂Ω t continues to move until some time t * where the boundary collides with itself (centre). A new connected boundary is formed over time and the topology of Ω t is now connected.
the mesh structure. In all mesh refinements considered, the change in mesh topology is only performed at the discrete time-levels. Hence, to ease notation, the superscript used to denote time-steps is omitted. Modified discrete functions, operators, and vectors are denoted using a hat symbol.

Contact scenarios
Here we present two contact scenarios that are numerically investigated in Section 7. The first scenario concerns the collision of the moving boundary with itself whereas the second situation involves collision with fixed geometric obstacles.
Self-intersection handles the situation where two parts of the moving boundary collide with each other. Typically, a remeshing is required in this instance. By using a VEM, the remeshing can be kept local and simple for colliding elements. A motivational case for a disconnected initial condition of the PME is given in Figure 1. Moving mesh finite element simulations of this type of problem are presented and discussed in [43].
Obstacle contact is encountered when the evolution of Ω t is obstructed by Figure 2: A demonstration of the obstacle contact problem. The initial condition is contained in a region of obstacles (left). After some time t * the moving boundary collides with the first obstacle (centre). A time-dependent interface now forms between the obstacles and ∂Ω t (right). external obstacles. An example of this is the presence of a solid phase in porous media. Figure 2 presents an example of collision with impermeable obstacles. By using a collision detection and node insertion algorithm, the moving mesh is capable of simulating the contact and the interaction between the moving mesh and a set of obstacles. As with the self-intersection problem, the VEM allows for this with minimal changes to the mesh topology. Additionally, the VEM is capable of performing local changes to polygonal elements such that the mesh boundary can move around the object boundaries without requiring additional mesh refinements.

Collision detection
For detecting mesh contact we use an adaptation of the classical node-tosegment collision detection algorithm [32,55,54]. We only consider boundary mesh and obstacle edges and nodes. We consider triplets of points ( 1 and x t 2 form a time-dependent edge e t and x t 3 is boundary node disconnected from e t . This triplet is referred to as a node-to-edge pair. Given a set of linear velocities (ẋ 1 ,ẋ 2 ,ẋ 3 ), by defining the vectors the contact time ∆t * between the hyperplane of e t and x t 3 is given as the minimum positive root of the quadratic equation where the coefficients are defined by In practise, we choose the contact time that makes physical sense (e.g. we discard negative roots as infeasible contact time) and only admit a singular contact time value for ∆t * . In the case of no feasible contact times we set ∆t * = ∞ and when two feasible contact times are given by equation (60) we choose the minimum of the two. By solving equation (60) for a set of node-to-edge pairing, a set of contact times can be computed. In the context of the moving mesh method, each node-to-edge pair on the boundary of the mesh and (when given) obstacle mesh is considered. If any cases indicate contact, the time step is scaled down to the minimum contact time and the corresponding node-to-edge pair is marked for contact. The detection method is outlined in Algorithm 2 for a single node-toedge pairing.
, a set of nodal velocities (ẋ t 1 ,ẋ t 2 ,ẋ t 3 ), the current time step ∆t Solve equation (60) and set ∆t * to the minimum positive root; for each value of ∆t * ∈ R do Compute x t+∆t * 3 and e t+∆t * ; if ∆t * ∈ [0, ∆t] and x t+∆t * 3 ∈ e t+∆t * then Mark the node-to-edge pair for contact; else Set the contact pair to no contact; end end output: the node-to-edge contact pair, the contact time step ∆t *

Node Insertion algorithm
Since the mesh allows for general polygonal element shapes, the insertion of a new node into a mesh edge can simply be performed by adding a vertex to the polygons sharing that edge. Then, a solution value associated with the new vertex must be introduced which requires an interpolation technique between the old and new global discrete spaces. In [22], an elliptic reconstruction operator is employed to preserve the quality of the discrete spatial derivative of the PDE. Instead, here we choose to preserve the polynomial component of the solution Πρ h between refinements through a redistribution of µ n . The reason for this choice is that, by preserving the polynomial component of ρ h , the global mass conservation of θ n h is maintained. This would not be the case if interpolation of the degrees of freedom was employed instead.
On a given element E ∈ T h with an inserted node on the boundary, we impose that the polynomial component of the solution is preserved so that the reconstructionρ h ∈V (E) satisfiesΠρ whereΠ denotes the projection operator constructed on the refined VEM spacê V (E). The arguments in Section 4.2 can be easily modified to show that this approach conserves both locally and globally the mass of the discrete solution. When introducing a new node onto a mesh edge of a given element E under the assumption of equation (62), the local contribution toμ is given by, The algorithm for node insertion is given in Algorithm 3.

Self-intersection algorithm
Due once more to the VEM flexibility in element geometries, the self-intersection problem does not require any introduction of additional degrees of freedom. Instead, the local connectivity of the disconnected mesh is updated to include the new node-to-edge pairing. Then, the node insertion algorithm is applied to recompute the solution and monitor distribution. As the boundary node velocities are not arbitrarily set, small edges are likely to appear during node insertion. This has not presented any stability issues within the numerical experiments of section 6 and we expect that the method remains robust in the presence of degenerate edges [18]. The subsequent mesh velocity problem is then solved based on the updated mesh and corresponding virtual element space. A simple demonstration is provided in Figure 3; the method at a given time step t n is presented in Algorithm 4.

Obstacle contact and pivot node algorithm
We denote a time-independent polygonal discretization of the set of obstacles by O h and consider the mesh node-to-edge pairings of boundary nodes from T n h and edges of O h and vice versa. In the case of obstacle-node to mesh-edge contact, the node insertion Algorithm 3 is applied to introduce an additional degree of freedom to the system; we refer to this new node, which is a fixed point on the obstacle geometry, as a "pivot node".
For contact between T n h and O h a no-penetration condition on the nodal velocities is strongly imposed on the formulation of the potential Problem 2.3 Algorithm 4: Self-intersection input : A mesh T n h , a solution ρ n h , a velocity field v, a time step size ∆t. Apply Algorithm 2 for boundary node-to-edge pairs in T n h ; Update ∆t; Compute T n+1 h , µ n+1 , ρ n+1 h ; if any contact pairs are marked then find the element E ∈ T n+1 h which contains the marked edge; Apply algorithm 3; end Update the boundary conditions of ρ n+1 h ; and velocity reconstruction Problem 2.4; namely, v · n = ∇φ · n = 0.
Hence, movement tangential to the obstacle's boundary is allowed. This is except when a pivot node is introduced, in which case we constrain its velocity to zero to preserve the geometry of the interface between the domain and the obstacle. The obstacle contact algorithm is outlined in Algorithm 5.
Given that the pivot node mesh velocity is constrained to zero, it is possible for the other boundary nodes laying on the obstacle (which have experienced meshnode to obstacle-edge contact) to pass through the pivot node. When this occurs, the connectivity of the mesh is updated to transfer the pivot node from one mesh edge to another as well as swapping the boundary node from one obstacle face to another.
Detection for pivot node collision is performed using Algorithm 2 for connected mesh boundary nodes moving from one obstacle edge to another.
A node is considered to be on ∂Ω n F only if both boundary edges sharing that node are in contact with the obstacle. We define a node to be "connected" if it lies on an edge of the obstacles. Degrees of freedom associated to connected nodes are constrained by equation (64) in the mesh velocity computations. If a node is connected by mesh edges to other connected boundary nodes we consider this to be an "interface" node and consequently change the boundary conditions from Dirichlet to Neumann defined in problem 2.1 (i.e. we change the degree of freedom from ∂Ω n M to ∂Ω n F ). If a connected node is not also an interface node, the homogeneous boundary condition and no-penetration condition are maintained. The structure of the boundary conditions are updated once every time step.
When a mesh node x i and a pivot node coincide (while the mesh node is moving along the obstacle boundary) the test function associated to a pivot node is chosen to satisfy ϕ pivot ≡ ϕ i . In other words, we duplicate the original VEM basis function and add it to the new discrete space. Remark 6. In practice, for instances of a pivot node and mesh node coalescing, the pivot node and corresponding degree of freedom can be temporarily removed from the VEM at a fixed time t n then reinserted before the mesh and ALE forward Euler updates. The monitor distribution and ALE update values are equally split between the pivot and node entries. which contains the marked edge; Apply Algorithm 3 to introduce a pivot node x pivot ; set v = 0 at x pivot ; end if Mesh node to obstacle edge is marked then Set v · n = 0 at mesh node; Update Neumann conditions for ρ n+1 h ; end

Numerical results
We report a series of numerical tests for the velocity-based moving mesh virtual element method proposed in the previous sections. Firstly, we present basic convergence test results using a known similarity solution of the PME Problem 2.1 for specific choices of velocity and solution recovery. Then, we investigate the effect on the solution of the node insertion algorithms described in Section 6. Finally, we present demonstrations of the contact algorithms of Section 6.

Sample meshes
We have tested four different mesh types used to subdivide the initial domain; representative examples of each are shown in Figure 4. The first mesh is the Voronoi Tessellation produced by randomly sampling mesh seeds in the domain [44,4]. The second mesh is a Centroidal Voronoi Tessellation (CVT) produced by the Lloyd algorithm which smooths a given Voronoi tessellation such that the generator points are the barycentric coordinates for each polygon [28]. The MATLAB package PolyMesher [47] was used to produce these two mesh types. The third mesh is constructed by overlaying the domain with a grid of uniform squares and cutting the mesh along the boundary. The last mesh type is a mixture of uniform Cartesian and polar tessellations. Note that the first three initial mesh types may present arbitrarily small edges and, moreover, arbitrarily small elements may appear near the boundary in the grid mesh type, as such potentially contradicting the mesh regularity assumptions stated in Section 3.1. In this respect, we note that the VEM is known to be quite robust, as we have also witnessed. We refer to the mesh size in each case as the largest element diameter in T 0 h .

The PME similarity solution
There exists a family of radially symmetric solutions on a given initial circular domain of radius r 0 for the Problem 2.1 defined in [42] and given by where d is the spatial dimension, r 0 is the initial radius, and Because of the nature of Equation (65), the solution is expected to have finite slope normal to the moving boundary for m ≤ 1 whilst, for m > 1, the solution presents an infinite slope normal to the boundary. Further properties of this analytical solution are discussed in [43,7].

Error computation
The numerical error is computed for both the solution and mesh by generalising the discrete approximations from [7]. An l 1 solution error is given by, while, for the mesh error, an l 1 norm is considered for the radial distance r(t) from the boundary of the mesh to the origin; thus Here, N B denotes the number of boundary nodes in the mesh and R n i denotes the radial distance from the origin of the i-th boundary node at time t n . A uniform ∆t that is small enough to ensure numerical stability is set for each initial sample mesh. Note that the meshes used in these numerical tests are not hierarchical. Furthermore, each Voronoi mesh is generated independently from randomly generated seeds. For each reduction of initial mesh size by a factor of 2, the time-step ∆t is reduced by a factor of 4 to ensure numerical stability. By reducing the time-step size by this factor we also expect that the temporal error to be O(h 2 ) when using the Forward Euler method for the refinement path of the four mesh types.

Convergence test
In this first convergence test the solution of the similarity solution for m = 1, d = 2, and r 0 = 0.5 is compared against the numerical solution for t = t 0 + T . The method is tested on each of the four mesh types for a circular domain. The time step sizes for the coarsest meshes are chosen according to where h mean is the average element diameter of the initial mesh T 0 h . In each mesh case we observe that the initial time-step size is approximately 1e − 4. From the initial time-step sizes we reduce ∆t by a factor of 4 each time the mesh is refined, which corresponds to the mesh size approximately halving with each refinement. In choosing the time step sizes we made conservative choices such that the numerical method was stable and presented the expected orders of convergence. A more robust approach would be to use adaptive time-stepping schemes but in this work we present convergence results for a uniform reduction in the time-step. As shown in Figure 5, second order accuracy is observed for the solution error for all mesh cases when T = 0.01. In the case of the Voronoi and Grid meshes, the empirical order of convergence (EOC) is less smooth compared to the CVT and mixed mesh types. This is most likely due to the weaker shape regularity of elements in these mesh types, but further studies are required. The mesh error EOC appears to have a long pre-asymptotic regime: it grows monotonically in all cases with the final computed values ranging between 1.55 (Voronoi) and 1.83 (mixed). This is consistent with finite element approximations of Darcy flow which observed the order of convergence of the velocity field to be lower than that of the pressure field [19,39]. Conservation of mass in the numerical solution is observed up to machine precision in all test cases.
Remark 7. Settingμ = 0 produces similar results to those reported in Figure 5. This is referred to as a "direct recovery" approach in the literature, see [7]. For brevity, the corresponding results are not presented here.

Node insertion convergence test
Our next numerical experiment considers the case of the one-dimensional PME extended in the x direction to a two-dimensional problem. This experiment has two interesting features. Firstly, the initial domain is geometrically exact (Ω 0 ≡ Ω 0 h ) unlike the circular meshes. Secondly, it allows us to test the obstacle contact and node insertion algorithm numerically against a known analytical solution derived from the one-dimensional case of equation (65). This is obtained by considering once again Equation (65) with the values m = 1, d = 1, r 0 = 0.5, and r = y on the initial domain given by Ω 0 = [−0.5, 0.5] 2 with initial condition ρ(x, 0) = 1 − 4y 2 .
The mesh is connected to two vertical planes at x = −0.5 and x = 0.5 with a no-penetration condition strongly imposed in the x direction; namely, x = 0 when |x| = 1/2. Solution snapshots at time T = 0.1 are shown in Figure 6 for the CVT mesh type. The mesh error is exclusively computed on the top and bottom faces of the rectangular domain.
In this test we only focus our attention on the CVT mesh type. We test the accuracy of the node insertion algorithm by including a discretization of the two planes into intervals in the y direction. In reference to the fixed domain PME 16 we defineΩ := [−1, 1] 2 and discretize the boundary ∂Ω into N uniformly spaced intervals to construct an N -gon. Vertices are then removed that intersect Ω 0 . This results in uniformly discretized intervals along both contact planes.
For the discretization of the contact planes N is set to an initial value of 32 and is doubled with each mesh refinement. Convergence results are reported in Figure 7. Here we observe second order accuracy in both the solution and mesh error. Also the mass is conserved, by design, up to machine precision throughout each test.

Contact demonstrations
We finally present two demonstrations of the node insertion algorithms for Problem 2.1 in challenging scenarios without known analytical solution. To ensure the quality of the mesh is preserved we introduce an alternative velocity field on the interior of the moving mesh based on the ALE approach. First, the method outlined in Section 5 is applied to approximate the Lagrangian boundary velocity v. The mesh velocity v for the internal node movement is then replaced byṽ, the harmonic extension of the Lagrangian boundary velocity: as such it is the solution to the problem In a VEM framework this alternative mesh velocity is computed by solving a standard Poisson's equation with Dirichlet boundary conditions [14]. Other choices for the interior mesh velocities include a modified monitor function M(ρ) [38], pseudo-elastic, and biharmoic formulations [45]. We remark that the quality of the mesh will still deteriorate over time. The purpose of these examples is to demonstrate the application of the node insertion algorithms. Optimising the choice of ALE velocity is left for future investigation.
In the first demonstration, we consider an initial condition of the PME that has a disconnected support such that self-intersection is expected to occur. The initial condition is given by An illustrative example of such initial condition is given in Figure 8 (top-left plot). The standard method is applied to simulate the PME for m = 1 with the contact detection Algorithm 2 applied at every time level to check for collision between elements. When contact occurs, Algorithm 4 is used to update the monitor distribution whilst the Dirichlet boundary degrees of freedom are flagged as interior degrees of freedom as the mesh connectivity evolves. Snapshots of the solution evolving over time are reported in Figure 8. The behaviour of the PME solution over time is in agreement with fixed mesh finite element approximations of this problem and similar benchmark tests performed for the PME in [43].
To demonstrate the obstacle contact algorithm we consider again the initial condition given by Equation (65). A set of obstacles are added to the computational domain in the shape of circles with random radii and centres. Each circle is approximated as a uniform polygon of comparable accuracy to the initial mesh. We tested the moving mesh VEM starting with a CVT type mesh made of 800 elements discretising the support of the initial solution. A few snapshots of the numerical solution are shown in Figure 9. Pivot nodes are inserted and removed along the contact. Mesh degeneration occurs after T = 0.2 and thus the simulation had to be terminated. In this case, and similar to finite element methods, a remeshing approach would rectify this issue.

Conclusions
In this paper we have combined, for the first time, a velocity-based moving mesh method with a virtual element method. This was achieved by extending the FEM formulation of the moving mesh method to a linear virtual element scheme on polygonal meshes. Numerical tests for the porous medium equation show that the proposed method obtains the same orders of accuracy as the original finite element approach. In fact, given that the linear VEM reduces to the linear FEM on triangular elements, our approach provides a natural extension of moving mesh FEM to polygonal mesh settings. Demonstrations of node insertion algorithms suggest that the virtual element method offers practical extensions to more complex problems. In particular, this work shows that it is straightforward to deal with situations where the moving boundary meets fixed obstacles or merges with other moving boundaries.
The VEM approach provides a flexible discretisation framework for adaptive moving mesh methods. For instance, VEM with curved elements are being developed with a level of generality out of reach for standard FEM, see eg. [11,37,3,29]. As such, the VEM is more suitable for the generalisation of moving mesh approaches in the higher-order setting, including for the solution of challenging problems such as fourth-order diffusion, phase-field, fluid-structure interaction and moving surface PDEs [7], including in three-dimensions [20,37,26]. We remark that these are largely open problems also for the more standard moving mesh finite element method. Extensions of the moving mesh VEM in these directions will be the subject of future works.