A simple history dependent remeshing technique to increase finite element model stability in elastic surface deformations

by   Jessica R. Crawshaw, et al.

In this paper, we present and validate a simple adaptive surface remeshing technique to transfer history dependent variables from an old distorting mesh to a new mesh during finite element simulations of elastic surface deformation. This technique allows us to reduce the error arising from excessive mesh distortion whilst preserving information about the initial configuration of the mesh and the history dependent variables. The transfer technique presented here constructs the initial configuration of the new mesh by considering the distortion incurred by the elements of the old mesh and projecting backward in time. Using this new initial configuration, the stress and strain over the new mesh can be easily calculated. After presenting the necessary steps to reconstruct the initial configuration, we show that this relatively simple transfer technique adds stability to finite element simulations and reduces the spatial error and the strain error across the domain. The novel transfer technique presented in this paper is easy to implement and provides a simple strategy to add stability to simulations undergoing large deformations.



There are no comments yet.


page 13


Boundary-Conforming Finite Element Methods for Twin-Screw Extruders: Unsteady - Temperature-Dependent - Non-Newtonian Simulations

We present a boundary-conforming space-time finite element method to com...

An unfitted finite element method for two-phase Stokes problems with slip between phases

We present an isoparametric unfitted finite element approach of the CutF...

Inf-sup stability of the trace P2-P1 Taylor-Hood elements for surface PDEs

The paper studies a geometrically unfitted finite element method (FEM), ...

Inf-sup stability of the trace P_2-P_1 Taylor-Hood elements for surface PDEs

The paper studies a geometrically unfitted finite element method (FEM), ...

An isogeometric finite element formulation for geometrically exact Timoshenko beams with extensible directors

An isogeometric finite element formulation for geometrically and materia...

On the Sobolev and L^p-Stability of the L^2-projection

We show stability of the L^2-projection onto Lagrange finite element spa...

Algorithms for 2D Mesh Decomposition in Distributed Design Optimization

Optimization of thin-walled structures like an aircraft wing, aircraft f...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

When investigating large geometric deformations using finite element models it is critical to maintain mesh integrity to ensure model accuracy [luo2008dealing]. Accuracy depends on both the refinement and the regularity of the mesh used to discretise the deforming domain. This is particularly important during large deformations in which increasingly excessive mesh coarsening and element distortion degrades the quality of the mesh throughout the simulation [javani2014consistent]. As an example, Figure 1a shows a deforming capillary plexus, discretised using an unstructured triangular mesh (the plexus is derived from images in Ghaffari et al. [ghaffari2015simultaneous] and the deformation is based on the model presented in Krüger et al. [kruger2011efficient], using the algorithm from Osborne and Bernabeu [osborne2018fully]). As the pressurised plexus grows, the mesh becomes increasingly distorted and elements elongate in the tangential direction (Figure 1b). Excessive element distortion over a mesh, as shown in Figure 1c, introduces large, unacceptable error to the simulation, thereby jeopardising the accuracy of the model and possibly invalidating any further analysis. To overcome these complications and to continue the simulation, successive remeshing throughout the simulation is unavoidable [zeramdini2019numerical].


[width=0.9trim=1.6cm 3.7cm 0cm 5cm,clip, right]Fig1B.pdf Element aspect ratioover timetime (sec)

Figure 1: 1.4A finite element model of vascular deformation with a fixed mesh showing mesh degradation. A constant internal pressure deforms a capillary plexus from an initial configuration at () to a deformed state at time (. As the plexus grows, the mesh progressively coarsens, and the elements become increasingly distorted. The median element aspect ratio (solid black line) drops over the duration of the simulation, while the interquartile range (shaded region) expands.

Following the generation of a new mesh in the remeshing process, the simulation variables from the old degenerated mesh must be transferred to the new mesh. One can either completely recompute the simulation variables or transfer the variables from the old mesh to the new mesh [zeramdini2019numerical]. If the optimal mesh configuration changes continuously thoughout the simulation (as is the case in the example shown in Figure 1) the second approach is generally considered preferable. However, doing so is a delicate matter, as inadequate determination of the updated state variables can introduce remeshing error which will adversely affect the simulation [dureisseix2006information]. This process of transferring the history-dependent variables from the old to new mesh is known as history-dependent remeshing. History-dependent remeshing techniques have been developed over the past 30 years prompted by the need for increasing model accuracy [dureisseix2006information, srikanth2000shape, kumar2015parallel, hinton1974local, peric1996transfer, zeramdini2019numerical, zienkiewicz1992superconvergent, zienkiewicz1992superconvergent2, zienkiewicz1992superconvergent3]. These transferal techniques can be broadly categorised into two groups; (1) techniques to transfer discontinuous variables which are continuous within each element and stored at Gauss points (referred to here as techniques) such as stress and strain ( variables); and (2) techniques to transfer continuous variables stored at nodes (referred to here as techniques) such as displacement, velocity and temperate ( variables) [kumar2015parallel].

transfer techniques are more straight forward than

techniques, and are appropriate when describing continuous fields where information is stored at the nodes. These techniques typically utilise interpolation or extrapolation techniques to attain a continuous field over the original mesh from the nodal values

[dureisseix2006information]. This continuous field can then be used to compute the nodal values over the new mesh. Various papers have extended upon this simple concept to increase accuracy and robustness of the technique using the least squares method and shape functions [bussetta2010comparison, kumar2015parallel, pere2002mapping].

Owing to the discontinuous nature of variables such as stress and strain, techniques are inappropriate to transfer variables between meshes. Indeed, the transfer of variables requires considerably more sophisticated techniques. techniques typically (i) enforce the balance equation in a weak or strong sense, (ii) use first or second order interpolation techniques, or (iii) are based on nodal or element patches [kumar2015parallel]. Following the first approach (i), Javani et al. [javani2014consistent] transferred a minimum set of variables from the old to the new mesh, and then reconstructed the remaining variables via the governing constitutive and balance equations. This procedure ensured the consistency of the constitutive and equilibrium equations were preserved over the new mesh.

First and Second order interpolation techniques (ii) generally follow one of two approaches; (a) direct interpolation [kucharik2003efficient, liszka1984interpolation]; and (b) nodal recovery based techniques [peric1996transfer, lee1994error, dureisseix2006information]. Direct interpolation (a) approximates piecewise continuous fields over the new mesh by using a set of field values at Gauss points on the old mesh in the region of the considered Gauss point of the new mesh to build a local and continuous description of the variable [lobov2011dll4, liszka1984interpolation]. An example of the direct interpolation method was presented by Brancherie et al. [brancherie2008consistent], who proposed a field transfer operator based the Moving Least Squares and the diffuse approximation methods to reconstruct the stress field over the entire domain. Nodal recovery based techniques (b) transfer the variables by firstly projecting the values of the Gauss points to the nodal points of the old mesh and build a continuous interpolation of the variable. From this continuous interpolation, the value at the nodes of the new mesh can be calculated. The variables at the Gauss points of the new mesh are then obtained from the nodal points using shape functions. Nodal recovery based techniques may not always be appropriate; in 2007 Khoei et al. [khoei2007superconvergence] found that these techniques violate the equilibrium of the systems, even if the mesh is unchanged.

One of the most popular transfer techniques is the superconvergent patch recovery (SPR) technique (iii) proposed by Zienkiewicz et al. [zienkiewicz1992superconvergent, zienkiewicz1992superconvergent2] in 1992. Using the SPR technique, for each node, a continuous polynomial expansion is taken over a patch of elements sharing the given node to construct a global approximation of the variables across the mesh. At points with superconvergent properties, this polynomial expansion approaches the values provided by the finite element analysis. An example of this technique in practise was presented by Boussetta et al. [boussetta2006adaptive]

who used the SPR technique to produce a continuous description of the stress tensor over the domain, from which the required mechanical properties could be reconstructed across the new mesh. Over the years, many modified and alternative SPR techniques have been proposed to increase the accuracy and robustness of this technique. Boroomand

et al. [boroomand1997recovery] proposed the Recovery by Equilibrium in Patch technique. This technique is a modification of the classic SPR technique that avoids the need to specify superconvergent points. Gu et al. [gu2004modified] improved the performance of the SPR technique for nonlinear problems by using integration points as sampling points, rather than the traditional use of nodes as sampling points. Wiber et al. [wiberg1995improved, wiberg1997superconvergent] formulated the superconvergent patches by elements transfer technique by introducing element patches rather than nodal patches. This technique penalised violation of the equilibrium and included a least squares fit for the boundary conditions.

It is evident that techniques are far more simple and straight forward than techniques. Whilst robust, the currently available techniques are often complicated, computationally time consuming, or involve a significant implementation time [peric1996transfer, zienkiewicz1992superconvergent, hinton1974local, zienkiewicz1992superconvergent2, zienkiewicz1992superconvergent3]. Unfortunately stress and strain are defined at each element meaning techniques are ill-equipped to deal with the discontinuous nature of these variables, thus forcing the use of techniques. In this paper we present a simple alternative to the classic and techniques to transfer strain from a deforming mesh to the new mesh in finite element models. This is done by focusing our attention on the displacement of each node and the continuous distortion incurred by each element, variables, over the old mesh. From this we can construct a new initial configuration for the new mesh over the original geometry, , and then calculate the strain over the new mesh (a variable) without resorting to more complicated techniques. This novel history-depended remeshing technique is of particular use in finite element models of elastic deformation where the reverse deformation is not known. We have designed this novel technique to be relatively straight forward to implement and appropriate for iterative remeshing schemes. In the remainder of this paper we describe this novel technique in more detail, examine its accuracy, and demonstrate how it can be implemented in a model of vascular deformation.

2 Computational methodology

[width=0.6trim=5cm 6cm 1cm 4.5cm,clip]Fig2.pdf

Figure 2: 1.4 An example of the history dependent variable transfer technique following remeshing of a 2D deforming square. The initial geometry, , is discretised by (top left). The square is deformed by , resulting in a non-uniform mesh with distorted elements, (top right). In this example , however the deformation is usually unknown. Here is the position of Node in the initial configuration and is the position of Node in the deformed mesh. The deformed square, given by the boundary of , is remeshed with a new, independent mesh, (bottom right). History-dependent simulations require a knowledge of the initial configuration of the mesh. As such the initial positions for each node in are attained via a mapping, , of to the initial geometry producing (bottom left), where is the location of Node in the new mesh.

We begin with some initial geometry, , discretised by an initial mesh, . This mesh is comprised of a set of nodes, , with locations, , and elements, , (Figure 2, ). The mesh is subject to a deformation, , which is usually unknown, leading to a new configuration at a later time 111 Note that we don’t explicitly consider time at this point. Here we are considering only the initial state and the deformed state. Later we consider deformations over time., , with node locations (Figure 2, ). The elements in are distorted and no longer provide an appropriate approximation of the deformed geometry. This deformed geometry is then remeshed with a completely new set of nodes and elements to provide a more accurate representation of the geometry (Figure 2, ). The new mesh is denoted as , and is comprised of a new set of nodes and elements which are given by and , respectively where are the new node locations. In history-dependent models, such as models of elasticity, the state variables across the domain are dependent on a clearly defined initial configuration. As such, we must ensure that each node in the new mesh, , has an initial position defined on the undeformed geometry, , which we denote as where the initial node locations are given by . The set of nodes and elements define the undeformed new mesh (Figure 2, ).

[width=0.5]Fig3.pdf Original elementDeformed element

Figure 3: 1.4A schematic of the transfer process for a single node. Node (red) from the new mesh, located at , is shown in the deformed configuration within the nearest element, , from the old mesh (blue element, right). The relative position of Node within Element

is expressed using a change of basis with the edge vectors

and and the element normal (not shown, Equation (1)). By replacing these basis vectors with the equivalent vectors ( and ) from the undeformed state (black element, left), the position of Node , , over the original configuration can be determined.

2.1 Mapping the new mesh to the undeformed configuration

To determine , we map each node in to the undeformed geometry, . This is achieved by considering the relative position of each node within the nearest element, , in the old mesh, , and the deformation incurred by this element (Figure 3). The method to calculate the nearest element is discussed below (Section 2.2). It is assumed that the deformation within each element is linear [kruger2011efficient, newham2019finite, zauel2006comparison], however this assumption could be relaxed with further model development (see Discussion). The relative position, , of Node in

is expressed in the basis defined by three linearly independent vectors characterising Element



The three vectors, , and , uniquely characterise Element (Figure 3) and represent two edge vectors and the unit normal at time , respectively. Here the unit normal is defined as . The basis coefficients, and are determined from the set of simultaneous equations arising from a basis transformation from Cartesian coordinates to the new basis defined by Element . As we assume deformation across each element is linear, the deformation at each point within each element at any time can be described using the deformation of the basis vectors. Thus it follows that we can describe the position of Node in by interchanging the basis vectors in Equation (1) with the basis vectors characterising Element in the undeformed state,


Here , and are two edge vectors of Element in the undeformed state, and is the unit normal in this state. Using this relatively simple method, we are able to map each node from the new mesh, , to the initial geometry, , thus providing an initial configuration for each new element. This mapped initial configuration allows us to recalculate the history-dependent state variables, such as stress and strain, over the new mesh. In doing so, we are able to remesh the history-dependent domain during large elastic deformation and maintain a near optimal discretisation of the domain.

2.2 Nearest element search algorithm

[width=0.7]Fig4A.pdf Case

[width=1.1]Fig4B.pdf Case

[width=1.1]Fig4C.pdf Case

Figure 4: 1.4Examples for each case of the nearest element search (Algorithm 1). Case 1: Node (located at ) is contained within the triangular prism defined by the element with the closest centroid (). Case 2: Node resides in either prism defined by the elements (Elements and ) associated with the nearest edge, , but not necessarily the element with the nearest centroid (). Case 3: Node exists in the region between the two prisms, thus a plane, , is constructed to determine on which side of Edge Node exists. Here and define the unit normals for Elements and , respectively, and .

The new mesh, , is an improved approximation of the 3D surface defined by , and therefore by definition will not co-exist with . As such, we can not simply assume that all of the nodes in will exist within elements of . Here we describe a method to find the nearest element in for each node in (Algorithm 1). This process is divided into three successive cases. Case 1: For each node, , we begin by locating the nearest element centroid in (Element ). Element is accepted to be the nearest element if Node exists within the triangular prism defined by the element (Figure 4, Case 1). Case 2: If the test in Case 1 fails then we locate more potential nearest elements by identifying the nearest edge, Edge , in (Figure 4, Case 2). Edge is associated with two elements, and (external edges have only one associated element so we choose this element as the nearest element). If Node is contained, in the prism of either or , then we accept the relevant element. In Figure 4, Case 2, it is clear that the centroid (black star) of Element is the nearest centroid to Node , however Node resides in the prism defined by Element . If both Element and Element are rejected, then we need to consider a third case. Case 3: the node lies in the region between the two prisms defined by the elements. To establish whether Element or Element is nearest to Node , a plane, , dividing this region is constructed (Figure 4, Case 3), and we determine on which side of the plane Node exists, thus providing the nearest element. Plane is defined as the plane containing both Edge and the mean of the associated element normals, . Here and are the unit normals of elements and , respectively. The situation may arise where Node is contained in both prisms defined by elements and . In this case one should follow the process outlined for Case 3. The steps of this process are outlined in Algorithm 1, and the cases are illustrated in Figure 4. Over a mesh with sufficient curvature, many element prisms will overlap, and Node can be contained in multiple prisms. As such, it is important to ensure the selected element prism is also the nearest, as determined using the centroids and edges. The procedure to find the nearest element collapses down for 2D domains, where one simply needs to find the triangular element containing Node without concerns about the third dimension.

0.95 for Each Node in  do
       Find the nearest element centroid in , Element ;
       if Node is contained within the prism defined by Element  then
             Return Element ;
             Find the nearest edge in , Edge ;
             if Node is contained within the prism defined by either element associated to Edge , Element or  then
                   Return relevant element, or ;
             else if Node is contained within either both, or neither, prisms defined by either elements, or  then
                   Determine if node is above or below the plane, , containing Edge and , and orientated towards , where and are the unit normals associated with elements and , respectively.;
                   if Above  then
                         Return Element ;
                         Return Element ;
                   end if
             end if
       end if
end for
Algorithm 1 Finding the nearest element in for each node in .

3 Results

In this section we show that the history-dependent remeshing technique described in this paper provides consistent results for spatial deformation and strain energy across the domain after both 2D and 3D deformation. To do so, we examine the spatial error, , and strain error, , across that arises following remeshing. Additionally, we consider how this history-dependent remeshing technique improves simulation stability across the capillary plexus deformation shown in Figure 1.

3.1 Spatial information is preserved after history-dependent remeshing

To examine the spatial error across following remeshing, we consider a set of 2D meshes discretising ( square, Figure 2 ) with increasing mesh refinement of . Here we define refinement using the initial target edge length at mesh generation. Each mesh is deformed, remeshed once using a set of meshes with increasing refinement, and mapped back to the initial geometry, , as described above in the computational methodology and shown in Figure 2. The continuous deformation across the mesh is defined at the nodes, and is a variable. In this example, is deformed using , as shown in Figure 2, where is the position of Node , however in practice the deformation is not usually known. For each refinement, is deformed through and the spatial error

is considered between the deformed mapped mesh and the new mesh.

[width=trim=1.5cm 11cm 11cm 10cm,clip, left]Fig5.pdf Edge lengthSpatial error with increasing remeshedmesh refinement over 2D meshCoarse Fine

[width=trim=10.5cm 11cm 2cm 10cm,clip, right]Fig5.pdf Coarse Fine Edge lengthSpatial error with increasing initialmesh refinement over 2D mesh

Figure 5: 1.4 The median ( IQR) spatial error across the deformed square shown in Figure 2 following history-dependent remeshing. The spatial error is unchanged with increasing refinement of the remeshed from either a coarse initial mesh (yellow) or a fine initial mesh (red). The spatial error decreases with increasing refinement of the initial mesh with both a coarse (green) and a fine (blue), with only a slight difference between each. The green and blue dots below the -axis in denoted the refinement of the initial mesh, either coarse (, green) or fine (, blue). Similarly the orange and red dots below the -axis in denote the refinement of the remeshed mesh, either coarse (, orange) or fine (, red).

Figure 5 shows the median spatial error ( interquartile range, IQR) over for increasing refinement of either (a) or (b). In Figure 5a, the refinement of increases, while the refinement of held constant, being either coarse (yellow), or fine (red). The spatial error decreases by approximately two orders of magnitude when remeshing from a fine initial mesh, compared to remeshing from a coarse initial mesh. This is confirmed in Figure 5b where the refinement of the adapted mesh is held constant, either coarse (green), or fine (blue), while the refinement of the initial mesh is increased. The spatial error decreases with an accompanying narrowing of the IQR as the initial mesh refinement increases. There is little to no additional accuracy attained using a fine over a course one. In practice, this means that the refinement of the new mesh has little impact on the on the remeshing error, but increasing the refinement of the original mesh will reduce remeshing error and provide increased model accuracy. Note that in these (and all subsequent) plots the minimum and maximum error follow the same pattern as the IQR.

[width=0.8trim=-9cm 0cm 60cm 10cm,clip, left]Fig6A.png

[width=0.8trim=50cm 0cm 1cm 10cm,clip, right]Fig6A.png

[width=trim=1.5cm 11cm 11cm 10cm,clip, left]Fig6B.pdf Edge lengthCoarse Fine Spatial error with increasing remeshedmesh refinement over 3D mesh

[width=trim=10.5cm 11cm 2cm 10cm,clip, right]Fig6B.pdf Edge lengthCoarse Fine Spatial error with increasing initialmesh refinement over 3D mesh

Figure 6: 1.4 The median ( IQR) spatial error across a deformed cylinder following history-dependent remeshing. A cylindrical mesh is deformed by . The spatial error across is unchanged with increasing mesh refinement of , remeshed from either a coarse initial mesh (yellow) or a fine initial mesh (red). The spatial error across (both coarse, green, or fine, blue) reduces sharply with increasing refinement of the initial mesh. The green and blue dots below the -axis in denoted the coarse () and fine () refinement of the initial mesh, respectively. Similarly the orange and red dots below the -axis in denote the coarse () and fine () refinement of the remeshed mesh, respectively.

These results are confirmed in 3D, as shown in Figure 6, where a cylinder is deformed to a larger sinusoidal elliptical cylinder (Figure 6a), remeshed, and analysed as described above for the 2D example. Here the cylinder is deformed though , where and . Again, increasing refinement of has little effect on the spatial error (Figure 6b and 6c). For this 3D example, the spatial error declines with increasing refinement. Increasing the refinement of

is also accompanied by a narrowing of the standard deviation, as seen in 2D (Figure

5b). These results suggest that the spatial accuracy of the history-dependent remeshing technique is dependent only on the refinement of the initial mesh. This is unsurprising as additional spatial information can not be attained mid-way through a simulation by re-discretising the deformed geometry. Rather, the spatial accuracy of the original mesh limits the attainable information about the initial configuration.

3.2 The hyperelastic strain energy density across reproduces that of

We now examine the error introduced to the strain across the deformed, , and remeshed, , meshes. The strain, a variable, is defined on and is continuous over each element, but is discontinuous over the domain. Here we assess the error in the hyperelastic strain energy density subsequent to history-dependent remeshing to understand how simulations behave following this technique. We also consider error minimisation by optimising the initial, and adapted mesh refinement and the remeshing frequency. We model the deforming square domain, shown in Figure 2, as a hyperelastic membrane, where the stress-strain relationship can be derived from a strain energy density function, . Following Kruger et al. [kruger2011efficient] we adopt the Skalak strain energy density function [skalak1973],


where and are the area and strain dilation moduli respectively, and and are the 2D principal strain invariants, describing the strain, and the dilation state of the membrane, respectively. The Skalak strain energy density function was first developed by Skalak [skalak1973] in 1973 to describe the deformation of red blood cells, and has been since used in to model red blood cell movement [kruger2011efficient]. We assume each element is both mechanically isotropic and homogeneous, and that deformation is linear across each element. For further details on the numerical implementation of this model see Kruger et al. [kruger2011efficient].

[width=1trim=1cm 5cm 0cm 7cm,clip]Fig7A.pdf Strain across fixed mesh

[width=1trim=1cm 5cm 0cm 7cm,clip]Fig7B.pdf Strain across remeshed mesh

[width=trim=1cm 9cm 10.5cm 10cm,clip,left]Fig7C.pdf Edge lengthStrain error with increasingmesh refinement

[width=trim=10.5cm 9cm 1cm 10cm,clip]Fig7C.pdf Coarse Fine Edge lengthStrain error with increasingremeshed mesh refinement

Figure 7: 1.2 The median ( IQR) strain error across decreases with increasing initial mesh refinement. The Skalak strain energy density over the original mesh, , and the new mesh, , following deformation (, Figure 2) are shown in a) and b), respectively. c) The median ( IQR) strain error across (black) and the successive (purple) decreases with increasing mesh refinement. When remeshing, is ascribed the same edge length as its preceding , and as such the error of the associated meshes can be seen by comparing directly for each refinement level. d) The median ( IQR) strain error across , is reduced for all levels of refinement if the preceding is fine (red), rather than coarse (orange). For clarity, the preceding is identified with either a green circle (coarse, edge length 0.26) or a blue circle (fine, edge length 0.06), and the median strain error over the coarse and the fine are shown by the orange and red dashed lines, respectively.

Figure 7 shows the Skalak strain energy density across ((a), initial edge length ), and ((b), initial edge length ). Here the domain is remeshed only once following deformation, and we use and . The coarseness of the initial mesh is reflected in the Skalak strain energy density, which shows a pixelated, or granular, nature between elements (Figure 7a). This granularity is inherited by the adapted mesh, where the echo of the original strain energy density can be seen. This echo is particularly evident in the top right corner of the adapted mesh (Figure 7b) where element distortion is greatest and suggests that the refinement of the original mesh, and the remeshing frequency may have a large impact on the accuracy of the final mesh.

To examine the error in the Skalak strain energy density, we consider the difference between the exact, , and the numerical, (Equation (3)), strain energy density for each element (),

Owing to the simplicity of the prescribed deformation, we can derive the exact Skalak strain energy density at each element centroid from the Greens deformation tensor. For the deformation given by in Figure 2, the strain invariants in Equation (3) can be expressed as,

where is the displacement of the centroid of Element . For further details of the implementation and the derivation of the Skalak strain energy density, refer to Kruger et al. [kruger2011efficient] and Skalak [skalak1973], respectively. As we assume the deformation over each element is linear, the strain energy is constant within an element, and thus we can concisely compare the exact and numerical strain energy density at the centre of the element.

From Figure 7c, we see that the median strain error ( IQR) across (black) and (purple) reduces with increasing mesh refinement. Again, mesh refinement is determined by the edge length at the time of mesh generation. Both , and the successive were ascribed the same edge length at mesh generation, and as such the error of the associated meshes can be seen by comparing directly for each refinement level in Figure 7c. The median strain error over is of the same order of magnitude to that of its predecessor at each level of refinement, however is consistently larger than the strain error over the corresponding initial mesh. Adjusting the refinement of in isolation has little effect on the strain error, as seen in Figure 7d, where the initial mesh is fixed, either coarse (orange, edge length 0.26) or fine (red, edge length 0.06), with the refinement of increasing. For comparison, the strain error of the initial mesh for the coarse and the fine meshes are shown by the orange and the fine dashed lines, respectively. Again there is a notable decrease in the strain error between the coarse and fine initial mesh for each level of mesh refinement over the adapted mesh. The minimum and maximum strain errors follow the same pattern as the IQR.

3.3 Increasing remeshing frequency reduces error

The deformation examined above, shown in Figure 2, produced a heterogeneous three-fold dilation of the initial geometry, and caused significant element distortion, with of elements suffering an aspect ratio below (discussed further below in Section 3.4). Remeshing only once after large element distortion (defined here as an element aspect ratio less than ) is alone insufficient to increase strain accuracy. Though introducing stability to the simulation, this single remeshing event slightly increases the strain error (Figure 7). As such, we now examine the spatial error and the strain error over a deforming hyperelastic membrane while remeshing multiple times.

As described above in Section 3.2, we model the geometry shown in Figure 2 as a deforming hyperelastic membrane. We now remesh the domain at defined intervals (between remeshing events) and discretise the deformation, , over time, to produce the equivalent time dependent deformation function


Here is time in seconds, and is the duration of the simulation, and is the initial position of Node . As this simulation is run for seconds, the remeshing intervals are equivalent to remeshing frequencies (per second). Based on the results of previous sections, we use an initial edge length of for both the initial and remeshed meshes.

For each remeshing frequency, we examined the spatial and strain error, as shown in Figure 8. Note that a frequency of 1 corresponds to only remeshing at the end of the deformation whereas a frequency of corresponds to remeshing at the end and at equally spaced times throughout the deformation. For example a remeshing frequency of 3 corresponds to remeshing (with an initial edge length of ) at and seconds. The median ( IQR) spatial error peaks at a single remeshing event following excessive deformation (Figure 8a), reflecting the results shown in Figure 6. Following two remeshing events, the median spatial error falls by an order of magnitude and there is little further accuracy gained for larger remeshing frequencies. Figure 8b also shows a magnitude drop in the median ( IQR) strain error when the number of remeshing events increases past two. The strain error over the fixed mesh is shown at zero remeshing events, with a larger error and a wider IQR than all other number of remeshing events, with the exception of a single remeshing event. The strain error is increased following a single remeshing event, reflecting the results in Figure 7, and again occurring as a consequence of excessive deformation prior to remeshing. For both the spatial error, and the strain error the minimum and maximum error follows the same pattern as the IQR.

It is necessary to note that the required frequency of remeshing should be guided by the rate and spatial distribution of element distortion. Simulations with more rapid, or with slower element distortion than that shown in the the example here will require more or less frequent remeshing, respectively. We have examined remeshing frequency at constant intervals here for the sake of simplicity, and have shown that after some threshold excessive remeshing adds little further accuracy to the simulation, and that remeshing must occur prior to excessive element deformation.

[width=1trim=11cm 11.5cm 2cm 11.1cm,clip]Fig8.pdf Remeshing frequencySpatial error with increasingremeshing frequencya)

[width=1trim=2cm 11.5cm 11cm 11.1cm,clip]Fig8.pdf Remeshing frequencyStrain error with increasingremeshing frequencyb)

Figure 8: 1.4 The median ( IQR) spatial error (a) and strain error (b) for increasing remeshing frequencies across a square hyperelastic membrane deforming over time. Here remeshing frequency is measured as remeshing events per minute.

3.4 Element aspect ratios over the deformed mesh improves with remeshing at a small expense to the quality of the initial mesh

Mesh stability is dependent on the quality of the elements forming the mesh. The quality of a given element is often quantified using the element aspect ratio given by . Here and describe the radii of the inscribed and circumscribed circles, respectively [khan2020surface]. Note that there are several different definitions of the aspect ratio of a triangle, though all have a similar meaning. A higher aspect ratio is indicative of a high quality triangle, while a low aspect ratio is typical of a poor quality triangle for the purpose of a finite element model. An equilateral triangle has an aspect ratio of one, and as a triangle distorts, the aspect ratio falls, approaching zero for a triangle with no surface area. A mesh composed of elements with larger aspect ratios () will be more stable under deformation than a mesh deforming with distorted elements (elements with a low aspect ratio).

Remeshing produces a stable mesh with element aspect ratio closer to 1, as shown in Figure 9, thereby increasing the stability of the simulation. Figure 9 shows histograms depicting the distribution of the element aspect ratio across the initial undeformed mesh, (), the deformed mesh, (), the remeshed mesh, (), and the remeshed mesh mapped to the initial configuration (). As the original mesh deforms, the elements distort and the median aspect ratio falls from to , and the distribution is relatively uniform, indicating is a poor quality mesh (Figure 9b). To recover valid results from the simulation, remeshing is necessary to reinstate a high quality mesh, with a median element aspect ratio close to (Figure 9d). This high quality mesh discretising the deformed domain is introduced at a cost to the initial configuration, (Figure 9c). The aspect ratio across has a much larger distribution than that of

, with a left skew and a median aspect ratio of

. This left skewed distribution occurs due to the dependence of the mapping on the already distorted . Importantly, there are less distorted elements across compared to , and the distortion within these elements is less pronounced (Figure 9b and Figure 9c). The presence of highly distorted elements over contributes to the spatial error and strain error. Increasing remeshing frequencies, or remeshing prior to excessive distortion will reduce the variation in the element aspect ratio over , thereby reducing the spatial error and strain error, as seen in Figure 8.

[width=0.8trim=2.5cm 8.1cm 1.2cm 8.5cm,clip]Fig9.pdf Aspect ratioAspect ratioAspect ratioAspect ratio

Figure 9: 1.4 The element aspect ratio over each mesh. (top left), (top right), (bottom right) and (bottom left).


[trim = 2cm 3cm 50cm 3cm , width=0.55clip=true, left]Fig10A.png

Figure 10: 1.4A deforming capillary plexus modelled as a hyperelastic membrane with a fixed mesh (upper) or with periodic remeshing (lower). Periodic remeshing (remeshing every seconds) stabilises the mesh throughout the simulation, whilst the fixed mesh distorts. This is confirmed by the persistently high median ( IQR) element aspect ratio with period remeshing (red) and the fall of the element aspect ratio over the fixed mesh (black). Without remeshing (black), the total surface area diverges, however with periodic remeshing (red), the total surface area of the plexus grows to an equilibrium.

3.5 History dependent remeshing stabilises a hyperelastic model of capillary deformation

We now turn our attention to the capillary plexus shown in Figure 1, where instability arose from element distortion as the plexus deformed, corrupting the simulation. We model the plexus as a hyperelastic membrane, as described in Section 3.2, deforming with a constant internal pressure, and with either periodic remeshing (remeshing every seconds), or no remeshing (a fixed mesh). Figure 10a shows the plexuses at time seconds both with (lower) and without (upper) remeshing for comparison. It is again evident that without remeshing the mesh undergoes excessive distortion, with the median ( IQR) element aspect ratio falling below (Figure 10b). The plexus modelled with remeshing deforms smoothly, and maintains mesh stability (median element aspect ratio ). The distortion of the fixed mesh introduces error into the simulation, and drives divergence of the plexus surface area (Figure 10c, black), where convergence is expected, and seen during plexus deformation with remeshing (Figure 10c, red).

4 Discussion

In this paper we have presented a relatively simple history-dependent remeshing technique to stabilise finite element models of surface mesh deformation by removing the error arising with mesh degradation. Over a finite element mesh, we consider the domain as a patchwork of continuous linear elements, and utilise this concept to map any point back to its initial position. This mapping enables us to transfer the initial configuration to the new mesh. With this information, we can reconstruct the history dependent variables, such as stress and strain, instead of directly transferring state variables, as is tradition [peric1996transfer, zienkiewicz1992superconvergent, zienkiewicz1992superconvergent2, zienkiewicz1992superconvergent3]. These traditional variable transfer remeshing techniques currently available in the literature are highly technical and require significant code implementation [peric1996transfer, zienkiewicz1992superconvergent, zienkiewicz1992superconvergent2, zienkiewicz1992superconvergent3]. The remeshing technique described here provides a simple alternative to these more complex methods, whilst reducing the error accumulated during mesh distortion, and preventing the introduction of degradation related artefacts.

The mesh configuration is a key consideration when developing and implementing finite element models of deformation. If not handled appropriately, mesh distortion can introduce error into the simulation and sabotage the final results, as shown in Figure 1. This matter of mesh distortion is particularly problematic in large, or unidirectional deformations (Figure 1). In this paper we intentionally chose simple examples to analyse this remeshing technique because, while conceptually simple, they demonstrated that this remeshing technique handles 1) element distortion, 2) uniform and unidirectional extension and compression and 3) large deformation, with ease. The deformation shown in Figures 58 was known, indeed prescribed for the purposes of validation, however it is not normally known (Figure 10). As such, this technique recovers the reverse deformation for each element and node in simulations where the deformation is not necessarily known.

The history-dependent remeshing technique presented here removes mesh distortion related error, however, there is of course some error, although small, inherent to mapping the initial conditions from one mesh to another (Figures 68). This error decreases with increasing refinement of the original mesh, , suggesting that the remeshing error most likely arises from any spatial distance between any given new node in and the nearest element from . Moreover, this error associated with remeshing is sufficiently small compared to that resulting from mesh distortion (Figure 8a). The remeshing interval is also an important variable influencing the error associated with remeshing (Figure 8). That is to say, the degree of distortion incurred by the original mesh prior to remeshing impacts the accuracy of remeshing. Normally remeshing is triggered by the aspect ratio falling below a given threshold, for example . In most examples shown here, the mesh deformed to three times its original size, which we are considering to be a large deformation.

Whilst time consuming when applied to dense meshes, the history dependent remeshing technique described here is more efficient than using an excessive refinement over the initial mesh in an attempt to prevent mesh distortion. Moreover, increasing the mesh refinement alone will be insufficient to prevent element distortion in the case of unidirectional deformation. In addition, the more time consuming step in the remeshing process, that is the element search, could be expedited in future work by the introduction of spatial binning, in which only the elements from within the same predefined spatial bin, or region, as the node are considered.

In the development of this history dependent remeshing technique, we have relied on the assumption that all elements deform linearly. This assumption can be relaxed for models assuming higher order deformation across the element, by introducing higher order terms to the basis vectors () to encapsulate the non-linear nature of element deformation. It is expected that this extension will be relatively straight-forward to implement. Finally, this method has been developed and implemented for surfaces in 2D and 3D space, and further consideration should be given to 3D tetrahedral meshes, considering a more specific definition of and utilising the shape functions.

5 Conclusion

In this paper we have presented a novel history-dependent remeshing technique to prevent mesh distortion in deforming surface meshes. This technique relies on the assumption that the mesh can be considered as a patchwork of continuous elements, and uses this assumption to map the initial configuration from the old distorting mesh, to a new optimised mesh. In doing so, history dependent variables can be recalculated using the new initial conditions in a simple, clean manner. We have shown that the error associated with this technique is small and decreases with increasing refinement of the original mesh, and with an appropriate remeshing frequency.