Siddhant Kumar

Kostas Danas
email: konstantinos.danas@polytechnique.edu

Dennis M Kochmann

Enhanced local maximum-entropy approximation for stable meshfree simulations

Keywords: meshfree/particle methods, updated Lagrangian, maximum entropy, tensile instability, finite deformation

We introduce an improved meshfree approximation scheme which is based on the local maximum-entropy strategy as a compromise between shape function locality and entropy in an information-theoretical sense. The improved version is specifically designed for severe, finite deformation and offers significantly enhanced stability as opposed to the original formulation. This is achieved by (i) formulating the quasistatic mechanical boundary value problem in a suitable updated-Lagrangian setting, (ii) introducing anisotropy in the shape function support to accommodate directional variations in nodal spacing with increasing deformation and eliminate tensile instability, (iii) spatially bounding and evolving shape function support to restrict the domain of influence and increase efficiency, (iv) truncating shape functions at interfaces in order to stably represent multi-component systems like composites or polycrystals. The new scheme is applied to benchmark problems of severe elastic and elastoplastic deformation that demonstrate its performance both in terms of accuracy (as compared to exact solutions and, where applicable, finite element simulations) and efficiency. Importantly, the presented formulation overcomes the classical tensile instability found in most meshfree interpolation schemes, as shown for stable simulations of, e.g., the inhomogeneous extension of a hyperelastic block up to 100% or the torsion of a hyperelastic cube by 200 • -both in an updated Lagrangian setting and without the need for remeshing.

Introduction

The finite element method (FEM) has been widely successful in simulating the mechanics and physics of solids. Unfortunately, conventional FEM is often severely limited when solving problems involving complex geometries and, especially, large deformations that may lead to severe mesh distortions. In this latter context, problems associated with classical FEM are two-fold. First, formulating the governing equations in the initial or reference configuration (referred to as total-Lagrangian description) often becomes inapplicable in case of large deformations, where the initial mesh loses relevance. Second, when formulated in the updated-Lagrangian or Eulerian settings, severe mesh distortion may lead to entangled or ill-shaped elements. Solutions to the latter scenario have included adaptive remeshing and mesh refinement as well as the arbitrary Lagrangian-Eulerian (ALE) method [START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF], in which the mesh is updated independently of the geometry. Such solutions can be complicated and computationally expensive, particularly for adaptive remeshing and refinement, which is prone to errors when mapping stresses and strains from the old mesh to the new mesh.

By contrast, meshfree methods are more suited for simulations involving large deformations. Detailed reviews of meshfree methods can be found in Belytschko et al. (1996b) and [START_REF] Chen | Meshfree methods: Progress made after 20 years[END_REF]. Notable examples of meshfree methods include smoothed-particle hydrodynamics (SPH) [START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF][START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF], the element-free Galerkin method [START_REF] Belytschko | Element-free Galerkin methods[END_REF], reproducing kernel particle methods [START_REF] Liu | Reproducing kernel particle methods[END_REF], and material point methods [START_REF] Sulsky | A particle method for history-dependent materials[END_REF][START_REF] Sulsky | Application of a particle-in-cell method to solid mechanics[END_REF][START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF], all of which circumvent mesh-related problems by treating nodes as interacting particles. However, meshfree methods have problems that are different from those of FEM and involve, in particular, the treatment of discontinuities and boundary conditions as well as numerical instabilities. [START_REF] Belytschko | Smoothing and accelerated computations in the element free Galerkin method[END_REF] showed that all kernel-based meshfree methods experience a rank-deficiency instability that causes spurious modes. This instability can be eliminated by sampling away from the nodes for numerical integration, using, e.g., stress-point integration [START_REF] Dyka | An approach for tension instability in smoothed particle hydrodynamics (SPH)[END_REF][START_REF] Dyka | Stress points for tension instability in SPH[END_REF][START_REF] Randles | Smoothed particle hydrodynamics: Some recent improvements and applications[END_REF] or material-point integration [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF][START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF]). Additionally, [START_REF] Belytschko | Smoothing and accelerated computations in the element free Galerkin method[END_REF] showed that, when the discretization is formulated in an updated-Lagrangian setting, a purely numerical instability called tensile instability arises from the changing nodal spacing and the associated localization of the shape function support. Addressing the tensile instability usually involves adaptive evolution of the characteristic lengths of the approximants, which is not trivial for most approximation schemes such as, e.g., moving least-squares (MLS) approximants. Meshfree schemes also face challenges in accurately capturing material interfaces due to the smearing of discontinuous gradients induced by the non-local support [START_REF] Cordes | Treatment of material discontinuity in the element-free Galerkin method[END_REF][START_REF] Krongauz | EFG approximation with discontinuous derivatives[END_REF].

More recently, maximum-entropy approximation schemes have attracted interest among the class of meshfree particle methods. [START_REF] Sukumar | Construction of polygonal interpolants: a maximum entropy approach[END_REF] used the maximization of information entropy [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] to formulate meshfree interpolants on polygonal elements, while [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] used a Pareto compromise between locality of approximation and maximization of information entropy to create the local maximum-entropy (local max-ent) shape functions. The local max-ent scheme is particularly attractive because it provides a seamless transition between nonlocal approximation and simplicial interpolation on a Delaunay triangulation like in FEM. It also possesses the weak Kronecker property, which allows the direct application of boundary conditions. [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF] used the local maxent shape functions to develop the optimal transportation meshfree (OTM) method for simulating problems such as fluid-structure interactions and fragmentation [START_REF] Li | Large scale optimal transportation meshfree (otm) simulations of hypervelocity impact[END_REF]. Though these methods are promising for dynamic problems, they are less well-suited for quasistatic simulations like the ones reported here. Additionally, [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF] proposed the use of an isotropic adaptive evolution of shape functions to overcome the problems associated with tensile instability [START_REF] Belytschko | Smoothing and accelerated computations in the element free Galerkin method[END_REF]. Despite this adaptive evolution of shape functions, that scheme is still prone to tensile instability in case of anisotropic deformations (as will be demonstrated in Section 2.4). Further advances in the area of max-ent approximations include the convergence analysis of [START_REF] Bompadre | Convergence analysis of meshfree approximation schemes[END_REF], the variational formulation of the optimal support size of max-ent shape functions [START_REF] Rosolen | On the optimum support size in meshfree methods: A variational adaptivity approach with maximumentropy approximants[END_REF], max-ent schemes with arbitrary order of consistency [START_REF] González | A higher order method based on local maximum entropy approximation[END_REF], as well as tools to evaluate derivatives of max-ent shape functions near the boundary [START_REF] Greco | Derivatives of maximum-entropy basis functions on the boundary: Theory and computations[END_REF].

Here, we present an enhanced local max-ent scheme for stable, quasistatic meshfree simulations. In Section 2.1, we present modified local max-ent shape functions that are based on an anisotropic Pareto compromise between maximizing information entropy and minimizing shape function width, which will play a crucial role in eliminating the tensile instability under large deformations. In Section 2.2, we begin with a quasistatic total-Lagrangian formulation of the mechanical boundary value problem and transform it into a spatially discretized, quasistatic updated-Lagrangian formulation. In Section 2.3, we examine the tensile instability through an instructive one-dimensional toy problem and, motivated by the same, in Section 2.4 we present a new, anisotropic adaption scheme for the evolution of the shape function support that is able to cope with tensile instability. Next, in Section 2.5, we discuss the application of essential boundary conditions and propose a simple, computationally inexpensive scheme that leverages the weak Kronecker delta property of local max-ent shape functions to accurately capture discontinuous derivatives across material interfaces. Section 3 summarizes a suite of numerical benchmark simulations of representative boundary value problems. We demonstrate that the enhanced local max-ent scheme introduced here provides better convergence and better handling of severe distortions than FEM. It accurately captures discontinuous derivatives across material interfaces and, most remarkably, avoids the common tensile instability associated with anisotropically increasing nodal spacings (for better understanding of the numerical implementation, a pseudo-code for solving simple hyperelasticity problems is included in B). Finally, Section 4 concludes our investigation.

Enhanced local maximum-entropy interpolation in updated-Lagrangian setting

First introduced by [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF], the local maximum-entropy (or local max-ent) scheme belongs to the class of convex approximation schemes and provides a seamless transition between finite elements (FE) and meshfree interpolations. The approximation scheme is based on a compromise between minimizing the width of the shape function support and maximizing the information entropy of the approximation. Building upon the idea of anisotropic shape functions [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF], we here introduce an enhanced version of the original local max-ent scheme, which uses an anisotropic support to deal with tensile instability in an updated Lagrangian formulation.

Anisotropic local maximum-entropy approximation

Consider a finite set of n n distinct nodes in d dimensions, X = {x a , a = 1, . . . , n n } with nodal positions x a ∈ R d . The domain bounded by convex hull of the node set is denoted by Ω ⊂ R d . Let u : Ω → R be a function that we aim to approximate based on the known values of the function at nodes, {u a = u(x a ), a = 1, . . . , n n }. We wish to construct an approximation of the type

u h (x) = n n a=1 u a N a (x), x ∈ Ω, (1)
where N a : Ω → R denotes the shape function corresponding to node a. Shape functions are subject to the constraints

N a (x) ≥ 0 ∀ x ∈ Ω, (2a)
n n a=1 N a (x) = 1 ∀ x ∈ Ω, (2b)
n n a=1 x a N a (x) = x ∀ x ∈ Ω. (2c)
The first constraint (2a) ensures the non-negativity of shape functions1 . The second constraint (2b) enforces zerothorder consistency (i.e., constant functions are exactly approximated), and whereas third constraint (2c) enforces firstorder consistency and guarantees the exact interpolation of affine functions. These together ensure that the scheme is consistent under h-refinement. Local max-ent schemes that satisfy higher-order consistency are also possible [START_REF] González | A higher order method based on local maximum entropy approximation[END_REF] but the formulation becomes increasingly complex. Hence, we limit ourselves to first-order consistency for the scope of this contribution. We define the width of a shape function N a as 2

U a [N a , β] = Ω N a (x) x -x a 2 β dV = Ω N a (x) (x -x a ) • β (x -x a ) dV, (3)
where β ∈ R d×d is a constant, positive-definite metric tensor defining the Euclidean distance in Ω. Further, β is a tensorial analogue of the scalar locality parameter in the original formulation of [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF]. A measure of the total width for the set of shape functions N = {N a , a = 1, . . . , n n } follows as the sum of individual shape function widths, viz.,

U[N, β] = n n a=1 U a [N a , β] = Ω n n a=1 N a (x) x -x a 2 β dV. (4)
Minimum shape function support or maximum locality requires minimizing the functional U[N] with respect to all shape functions subject to constraints (2). The resulting scheme is equivalent to using linear interpolants on a Delaunay triangulation of the node set X [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF]. By interpreting the shape functions as probability distributions, the information entropy encoded in the shape functions evaluated at a point x ∈ Ω is defined as [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF])

H[N](x) = - n n a=1 N a (x) ln N a (x).
(5)

The total information entropy is obtained by integration over the entire domain Ω, i.e.,

H[N] = - Ω n n a=1 N a (x) ln N a (x) dV. (6)
Intuitively, the approximation scheme should be based on a minimum inference bias by the nodal positions, which is equivalent to the maximum information entropy [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. This requires maximizing the information entropy functional H[N] with respect to all shape functions subject to constraints (2). Consequently, the resulting shape functions have a global support that is too large for a reasonable approximation. [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] proposed a Pareto optimality between the competing optimizations (minimum locality vs. maximum entropy) to find the shape functions. In adopting their approach, we define the functional

F[N, β] = U[N, β] -H[N] = Ω n n a=1 N a (x) x -x a 2 β + N a (x) ln N a (x) dV, (7)
where the metric tensor β acts as a parameter controlling the Pareto optimality and is referred to the locality parameter in subsequent discussions. For given β, the local max-ent shape functions are thus computed as

N β = arg min F[N, β] s.t. (2) . (8)
For example, consider the simple isotropic case of β = βI, where β is a scalar parameter [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF].

In the limiting case of β → ∞, the problem reduces to minimizing locality, which results in affine interpolation on a Delaunay triangulation. In contrast, when β → 0, the problem reduces to maximizing information entropy with no regard to locality of the approximation. A detailed discussion for general anisotropic locality parameter is presented later in this section.

The structure of the minimization problem in (8) admits a pointwise optimization, so the shape functions at a point x ∈ Ω are given by

N 1 (x), . . . , N n n (x) β = arg min        n n a=1 N a (x) x -x a 2 β + N a (x) ln N a (x)        subject to              N a (x) ≥ 0, a = 1, . . . , n n n n a=1 N a (x) = 1 n n a=1 N a (x)x a = x. (9)
Analogous to the isotropic formulation of [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF], there exists a unique set of minimizers given by

N a (x) = 1 Z(x, λ * (x)) exp -x -x a 2 β + λ * (x) • (x -x a ) , a = 1, . . . , n n , (10)
where we define the partition function

Z : R d × R d → R as Z(x, λ) = n n a=1 exp -x -x a 2 β + λ • (x -x a ) , (11)
and a unique minimizer λ * as λ * (x) = arg min

λ∈R d ln Z(x, λ). ( 12
)
The derivation of ( 10)-( 12) with the newly introduced tensorial locality parameter β is summarized in A.

Again following [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF], the spatial derivatives of the shape functions are given by where

∇N a (x) = -N a (x) J (x, λ * (x)) -1 (x -x a ) (13)
J(x, λ) = ∂ 2 ln Z(x, λ) ∂λ ∂λ = n n a=1 N a (x, λ)(x -x a ) ⊗ (x -x a ) -r(x, λ) ⊗ r(x, λ) (14)
and

r(x, λ) = ∂ ln Z(x, λ) ∂λ = n n a=1 N a (x, λ)(x -x a ). ( 15
)
For a detailed derivation, the reader is referred to the analogous isotropic case presented by [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF]. Note that β is assumed to be constant when constructing the approximation and when computing the above gradients of the shape functions. The justification of this assumption will be discussed in Section 2.3, where the spatial discretization scheme is presented.

The transition from a scalar to a tensorial locality parameter is essential and allows for anisotropic adaptivity of the shape functions, which forms the basis for our approach to dealing with tensile instability, as will be demonstrated in Section 2.4. Note that, for variational problems, the locality parameter β can be included in the global functional to yield shape function supports that are optimal with respect to the total potential energy [START_REF] Rosolen | On the optimum support size in meshfree methods: A variational adaptivity approach with maximumentropy approximants[END_REF]. However, this introduces additional complexity, both theoretical and computational, and will not be pursued here. Instead, the locality parameter β is defined by the user and chosen based on the initial nodal spacing and the physics of the problem.

Figure 1 illustrates the transition of the shape function support from global to local and from isotropic to anisotropic, as the locality parameter β changes. We emphasize that, instead of plotting a particular shape function evaluated at 1, note that we do not plot a given shape function N a as a function of position x but we compute shape functions N a at a given point x.

different points in space as commonly shown in the literature, Figure 1 illustrates the shape function values of each node evaluated at a fixed point (thus demonstrating the influence of each nodal contribution at that location in space). This is required since shape functions are no longer associated with nodes only -viewed from different material points with different β-tensors, a shape function evaluated at a node may have distinct values for each material point under consideration. In addition, shape functions are only evaluated at material points and there is no unique definition of those at nodes. In the limit of any eigenvalue of β tending to +∞, the functional F in ( 7) is dominated by the shape function width. As a result, the shape functions reduce to simplicial interpolation on a Delaunay triangulation of the node set in the corresponding eigen-direction (see Figures 1(c) and 1(e)). In particular, Figure 1(e) illustrates global support in one direction and convergent linear interpolation in the perpendicular direction. Unlike the FE method, the local max-ent scheme does not satisfy the strong Kronecker property at the nodes, i.e., it is an approximation and not an interpolation. However, it does satisfy the weak Kronecker property at the boundary nodes [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF]. Figure 2 shows the shape function support when evaluated at points on the convex hull of the node set. When the evaluation point lies on an edge (see Figure 2(a)), shape functions of interior nodes vanish at the point and the approximation only depends on boundary nodes. Further, if the evaluation point coincides with a corner node, shape functions evaluate to zero for all nodes but the corner node itself. More generally for convex schemes like local maxent, the shape function support of interior nodes terminates at the convex hull and, consequently, the approximation on the boundary is independent of interior nodes. This fact is advantageous for straightforward application of essential boundary conditions. In contrast, non-convex approximation schemes such as smoothed particle hydrodynamics (SPH), element-free Galerkin method (EFG), reproducing kernel particle method (RKPM), etc. do not satisfy this property and require additional methods for imposing essential boundary conditions [START_REF] Bonet | Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations[END_REF][START_REF] Belytschko | Element-free Galerkin methods[END_REF][START_REF] Liu | Reproducing kernel particle methods[END_REF][START_REF] Fernández-Méndez | Imposing essential boundary conditions in mesh-free methods[END_REF].

The minimization problem in (12) can be solved in a few Newton-Raphson iterations. The derivative and Hessian matrix of the objective function are given by ( 15) and ( 14) respectively, and the Newton step

λ ← λ -J(x, λ) -1 r(x, λ) (16)
is iterated until convergence. In the limit of the approximation converging to linear interpolation along any direction (some eigenvalue of β → +∞), the Hessian matrix becomes singular when the guess for λ is far away from the unique minimizer λ * . [START_REF] Li | The optimal transportation method in solid mechanics[END_REF] suggested the use of a steepest descent technique in this case. Alternatively, we found that the regularized Newton method proposed by [START_REF] Polyak | Regularized Newton method for unconstrained convex optimization[END_REF] and applied to local max-ent by [START_REF] Foca | On a Local Maximum Entropy interpolation approach for simulation of coupled thermo-mechanical problems[END_REF] gives a faster convergence in the aforementioned limit of large eigenvalues of β. We consider the modified objective function

F(x, λ, ζ) = ln Z(x, ζ) + 1 2 ∂ ln Z(x, λ) ∂λ ζ -λ 2 . ( 17
) Since ln Z(x, ζ) is a convex function and ∂ ln Z(x,λ) ∂λ ζ -λ 2 ≥ 0, the regularized function F(x, λ, ζ) is strongly convex in ζ.
The minimizer of ( 12) is given by

λ * (x) = arg min ζ∈R d F(x, λ, ζ)| ζ=λ . (18)
The derivative and Hessian matrix of the objective function in (18) are obtained as

∂F(x, λ, ζ) ∂ζ ζ=λ = r(x, λ), ∂F(x, λ, ζ) ∂ζ∂ζ ζ=λ = J(x, λ) + r(x, λ) I. (19)
Clearly, the regularized Hessian is non-singular even when J(x, λ) is singular. Therefore, the regularized Newton step for any λ λ * , viz.

λ ← λ -J(x, λ) + r(x, λ) I -1 r(x, λ), (20)
is expected to provide more robust convergence for the aforementioned limiting case.

Quasi-static updated-Lagrangian formulation

In a classical total-Lagrangian description, the initial configuration is always the reference at any time or load steps. However, in simulations involving large deformations and distortions, the initial reference configuration often loses its meaning. An updated-Lagrangian formulation overcomes this limitation by using the deformed configuration of the previous step as the reference in an incremental fashion. For instructive purposes, we first describe a quasi-static variational problem in the total-Lagrangian setting and then reformulate it in an updated-Lagrangian framework. The proposed scheme is general to all approximation schemes and is not limited to local max-ent shape functions. Also, even though we limit our derivation to finite elasticity, the method can easily be extended to a more general class of variational material models (e.g., by using variational constitutive updates for inelasticity [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF]).

We assume an incremental setting in which Ω n refers to the n th configuration of a body Ω (likewise, all fields are denoted by subscripts (.) n ), and in particular n = 0 corresponds to the initial configuration with domain Ω 0 . Consider the total potential energy at step n with respect to the initial domain Ω 0 :

I[ϕ n ] = Ω 0 W(F n ) dV - Ω 0 ρ 0 B n • ϕ n dV - ∂Ω N 0 Tn • ϕ n dS , (21)
where W : GL + (d) → R denotes the strain energy density, F = Grad ϕ is the deformation gradient, ϕ : Ω → R d is the deformation mapping, B represents body forces and T constant surface tractions applied over the Neumann boundary

∂Ω N 0 ⊂ ∂Ω 0 .
All quantities are defined at a step n ≥ 0. We seek solutions as infimizers of the total potential energy3 , specifically,

ϕ n = arg inf I[ϕ] : ϕ ∈ U n , (22)
where

U n = ϕ ∈ H 1 (Ω 0 ) : ϕ = φn on ∂Ω D 0 (23)
defines the set of admissible solutions with φn denoting essential boundary conditions at step n defined on the Dirichlet boundary ∂Ω D 0 ⊂ Ω 0 . Without loss of generality, we make two assumptions: the initial configuration is undeformed (i.e., F 0 = I), and the initial and deformed coordinate systems coincide. The latter allows us to define the total displacement field at step n as u n (x) = ϕ n (x 0 ) -x 0 , and we seek solutions in terms of the displacement field. We thus write the deformation gradient as

F n = ∂ϕ n ∂x 0 = I + ∇ 0 u n , (24)
where ∇ m = ∂ ∂x m is the spatial derivative with respect to the m th configuration.

Let us introduce a discretized approximation based on the above modified local max-ent scheme, which gives

u n (x 0 ) ≈ u h n (x 0 ) = n n a=1 u a n N a (x 0 ) ⇒ F n ≈ F h n = I + n n a=1 u a n ⊗ ∇ 0 N a (x 0 ), ( 25
)
where the shape functions N a (x 0 ) are computed in the initial configuration Ω 0 . The set of all nodal displacements are denoted by U n = {u 1 n , . . . , u n n n , }. Computing the first variation of ( 21) and inserting the approximate fields at step n yields the discretized stationary condition in terms of nodal force balance:

f a n (U n ) = f a int,n (U n ) -f a ext,n = 0, a = 1, . . . , n n . ( 26
)
The internal force vector acting on node a is given by

f a int,n (U n ) = Ω 0 P(F n ) • ∇ 0 N a (x 0 ) dV (27)
with the first Piola-Kirchoff stress tensor P = ∂W/∂F. Analogously, the external force on node a is

f a ext,n = Ω 0 ρ 0 B n N a (x 0 ) dV + ∂Ω N 0 Tn N a (x 0 ) dS . ( 28
)
Solving the quasi-static system of equations in ( 26) calls for a robust and efficient iterative solver, which for optimal convergence rates uses the consistent tangent such as the family of Newton-Raphson methods.

To this end, the submatrices of the tangent matrix T n and its components are computed as

T ab n = ∂ f a n ∂u b n (U n ), T ab n ik = Ω 0 C i jkl (F n ) ∇ 0 N a j (x 0 ) ∇ 0 N b l (x 0 ) dV - ∂ f a ext,n i ∂ u b n k , (29)
where we used indicial notation with the classical summation convention and the incremental stiffness tensor C i jkl = ∂P i j /∂F kl . Once the nodal displacements are found, the total strain energy of the domain at step n is given by the first term in (21), i.e., E n (U n ) = Ω 0 W(F n ) dV. This completes the total-Lagrangian description at the n th step with respect to the initial configuration as reference.

To transform the formulation into an updated-Lagrangian framework, we introduce an incremental deformation gradient F n→n+1 which maps the configuration from step n to n + 1 as

F n→n+1 = ∂ϕ n+1 ∂x n = I + ∇ n (u n+1 -u n ). ( 30
)
Under the assumption that the step size is sufficiently small, the incremental deformation gradient can be related to total deformation gradient by the composition

F n+1 ≈ F n→n+1 F n . (31)
Note that the mapping of gradients between the reference and current configurations is

∇ 0 (•) = ∂(•) ∂x 0 = ∂x n ∂x 0 ∂(•) ∂x n = F n ∇ n (•), (32)
while the density and volume integral are transformed, respectively, according to

ρ 0 = J n ρ n , Ω 0 (•) dV = Ω n 1 J n (•) dV with J n = det F n . ( 33
)
To arrive at a consistent updated-Lagrangian description, we start with the total-Lagrangian description at step n + 1, analogous to ( 27), but now computing all shape functions with respect to the Ω n instead of Ω 0 . By exploiting the transformations (32) and ( 33), the internal nodal forces at step n + 1 are thus given by

f a int,n+1 (U n+1 ) = Ω 0 P(F n+1 ) • ∇ 0 N a (x n ) dV = Ω n 1 J n P(F n→n+1 F n )F T n ∇ n N a (x n ) dV = f a int,n+1 (∆U n ) ( 34
)
where ∆U n = U n+1 -U n is the vector of nodal displacement increments from Ω n to Ω n+1 . Note that the term 34) can be interpreted as the pull-back of the Cauchy stress tensor from Ω n+1 to Ω n or, alternatively, as the push-forward of the first Piola-Kirchoff stress tensor from Ω 0 to Ω n . This formulation is a compromise between classical incremental/rate formulations and variational modeling which typically uses a total-Lagrangian setting.

J -1 n P(F n→n+1 F n )F T n in (
The transformation of the external nodal forces (28) occurs in a similar fashion (using ( 32), ( 33)), leading to

f a ext,n+1 = Ω n ρ n B n+1 N a (x n ) dV + ∂Ω N n Tn+1 N a (x n ) dS , (35)
where B n+1 and Tn+1 are now defined in the n th configuration. Further, the tangent matrix components of ( 29) and the total strain energy transform into, respectively,

T ab n+1 ik = Ω n C i jkl (F n→n+1 F n ) F q j ∇ n N a q (x n ) F rl ∇ n N b r (x n )J -1 n dV - ∂ f a ext,n+1 i ∂ ∆u b n+1 k (36) and E n+1 (∆U n ) = Ω n 1 J n W(F n→n+1 F n ) dV. ( 37
)
To construct a numerical approximation of the spatial integrals in the above expressions for energy, forces and tagent matrix, we introduce a second set of n p so-called material points which discretize the mass distribution in the domain according to [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF]

ρ n (x n ) = n p p=1 ρ p n V p n δ(x n -x p n ) = n p p=1 m p δ(x n -x p n ), (38)
where ρ p n , V p n , and m p are, respectively, the mass density, volume, and mass associated with the p th material point located at x p n ∈ Ω n . In order to distinguish nodes from material points, in the following superscripts {a, b} and p are reserved for nodes and material points, respectively.

Material points are meshfree analogues of quadrature points in conventional FEM, where sampling is performed for the purpose of numerical integration. At the initial step, the p th material point is assigned a density ρ p 0 = ρ(x 0 ) based on the mass density field ρ in the initial configuration. The initial material point volume V p 0 is computed from an ad-hoc triangulation of the initial configuration (which is discarded afterwards, as discussed below), akin to the volume of an element in the FEM context. Here, we enforce the local conservation of mass [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF], so that each material point carries a constant mass (m p = ρ p n V p n = const. ∀ n ≥ 0), while the density and volume are updated accordingly at every step, resulting in

ρ p n+1 = ρ p n det F p n→n+1 = ρ p 0 J p n+1 , V p n+1 = V p n det F p n→n+1 = V p 0 J p n+1 , (39)
where F p n→n+1 = F n→n+1 (x p ) is the incremental deformation gradient at the material point p. Nodes do not carry mass nor volume and solely provide the kinematic interpolation required at the material points, while the constitutive laws are evaluated at the material points only. Figure 3 illustrates the evolution of material points and nodes in the totaland updated-Lagrangian settings.

Integration over the domain Ω n is now approximated by a sum over the material points weighted by the their respective volumes:

Ω n φ(x n ) dV ≈ n p p=1 φ(x p n )V p n , (40) Initial configuration
. . .

(n + 1) th configuration F 0→1 F 1→2 F 2→3 F n→n+1 F 1 F 2 F n+1
Figure 3: Illustration of the total-Lagrangian description (dashed arrows) where the current configuration is referenced to the initial configuration via the total deformation gradient F n+1 , and of the updated-Lagrangian description (solid arrows) where the current configuration is referenced to the previous configuration via the incremental deformation gradient F n→n+1 . The red and black points denote material points and nodes, respectively.

where φ : Ω n → R is an arbitrary function. When applied to the meshfree governing equations, the update-Lagrangian description with material point sampling leads to the following set of equations:

F p n→n+1 (∆U n ) = I + n n a=1 ∆u a n ⊗ ∇ n N a (x p n ) and F p n+1 = F p n→n+1 F p n , (41a)
f a int,n+1 (∆U n ) = n p p=1 V p n J p n P(F p n→n+1 F p n )F p n T ∇ n N a (x p n ), (41b)
f a ext,n+1 = n p p=1 ρ p n V p n B n+1 (x p n )N a (x p n ) + ∂Ω N n Tn+1 N a (x n ) dS , (41c)
T n+1 ab ik (∆U n ) = n p p=1 C i jkl (F p n→n+1 F p n ) (F p n ) q j ∇ n N a q (x p n ) (F p n ) rl ∇ n N b r (x p n ) V p n J p n - ∂ f a ext,n+1 i ∂ ∆u b n k , ( 41d
) E n+1 (∆U n ) = n p p=1 V p n J p n W(F p n→n+1 F p n ). ( 41e
)
2.3. Tensile instability [START_REF] Belytschko | Smoothing and accelerated computations in the element free Galerkin method[END_REF] provided a unified analysis of stability for meshfree particle methods. In particular, they identified two distinct instabilities -(i) an instability due to presence of spurious modes introduced by the rank-deficiency of the discretization operator, and (ii) tensile instability. The first instability is particular to nodal integration, and stress-point integration has been shown to eliminate this type of instability; see, e.g., [START_REF] Dyka | An approach for tension instability in smoothed particle hydrodynamics (SPH)[END_REF][START_REF] Dyka | Stress points for tension instability in SPH[END_REF][START_REF] Randles | Smoothed particle hydrodynamics: Some recent improvements and applications[END_REF][START_REF] Belytschko | Smoothing and accelerated computations in the element free Galerkin method[END_REF]. Similar to stress-points, material point integration also alleviates the instability due to rank-deficiency by performing quadrature away from the nodes. By contrast, tensile instability is more nefarious and appears when updated-Lagrangian/Eulerian discretization kernels are used. Tensile instability can be eliminated by using (total) Lagrangian kernels [START_REF] Belytschko | Smoothing and accelerated computations in the element free Galerkin method[END_REF], but this defeats the purpose of using particle-based meshfree methods. Here, we study the tensile instability with local maxent shape functions in a simplified 1D setting, which will later motivate our solution for dealing with this type of instability in general 3D problems.

As an instructive example, consider a 1D chain of n n = 50 equally-spaced nodes whose positions in the initial configuration Ω 0 = [0, 1] are x a 0 = a-1 n n -1 , a = 1, . . . , n n . Assume there is initially only a single material point (n p = 1), located at x p 0 = 0.5 and carrying an initial volume V 1 0 = 1 equal to the unit length of the chain, and a constant unit 10 stiffness (linear elastic constitutive law). The chain is now stretched by a ratio λ > 1, while keeping the locality parameter β 1 constant. Because the shape functions are dependent on the nodal positions, increasing the nodal spacings under extension results in a localization of the shape function support at the material point, see Figure 4(a); i.e., fewer and fewer nodes contribute to the approximation of the fields of interest at the material point.

Recall that the internal nodal forces (41b) are given by the product of the first Piola-Kirchoff stress tensor with a discretization kernel

F p n T ∇ n N a (x p n )V p n /J p n = F p n T ∇ n N a (x p n )V p 0 . Figure 4(b)
shows the changing nodal kernel values with increasing nodal spacing (i.e., increasing applied stretch λ). From the point of view of any particular node, increasing the stretch λ leads to changes in its kernel value (evaluated at the material point p), which first increases strongly and then decays to zero as the shape function support increasingly localizes, as illustrated in Figure 4(c). Internal nodal forces can be computed using the kernel values and the stress (from linear elastic constitute law). The non-monotonic behavior of the kernel (with respect to stretch ratio λ) results in a non-monotonic nodal force that starts to decrease after a certain strain, as shown in Figure 4(d). This phenomenon of non-physical stiffness is called tensile instability, and is a purely numerical artifact arising due to the changing nature of the discretization kernel; see [START_REF] Swegle | Smoothed particle hydrodynamics stability analysis[END_REF] for a related discussion in the context of SPH. In summary, increasing nodal spacings causes a localization of the shape function support, which in turn leads to tensile instability. . As suggested by [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] and [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF], the tensile instability can be overcome by evolving the locality parameter of the local max-ent shape functions at every step. In their isotropic max-ent framework, those authors introduced a non-dimensional constant γ such that β = γ/h 2 , where h is some appropriate measure of average nodal spacing at a given time or load step4 . In the simple problem of a 1D node set as in Figure 4, evolving β as per the aforementioned scheme, i.e.,

β p n+1 = γ h p n+1 2 = β p n       h p n h p n+1       2 , (42)
produces exactly coincident shape function support at every step and hence correctly reproduces constant nodal force kernels and linear nodal force responses as shown in Figure 5. This adaptivity scheme can then be extended to higher dimensions as

β p n+1 = 1 (J p n→n+1 ) 2/3 β p n ( 43
)
where β p n is the locality parameter at n th step, and where it is assumed that (J p n→n+1 ) 1/3 is a reasonable estimate of the relative average nodal spacing around material point p. Although well suited for 1D, the isotropic nature of this scheme does not properly account for direction-dependent changes in nodal spacing and hence is severely limited in the amount of anisotropic deformation that can be simulated before the onset of tensile instability. As a simple example, consider a uniaxial extension test in 2D, which leads to strongly increasing nodal spacings in the extension direction, whereas the Poisson effect results in nodal spacings shrinking in the transverse directions at the same time.

Local updates of the locality parameter of type ( 43) have an essential consequence, which is key to the proposed methodology and a departure from classical discretization schemes including the original maximum-entropy approximants of [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF]. Here, β p is defined and updated only at each material point p, where it is needed to evaluate shape functions N a (x p n ) and their gradients ∇N a (x p n ) for the numerical integration of energy, forces and stiffness tensors according to (40). Therefore, neither a shape function N a nor a node a is associated with a specific β but, depending on the material point p at x p where shape functions are to be evaluated, all shape functions N a (a = 1, . . . , n) are computed using the very same β p of that material point. As a practical consequence, shape functions cannot be evaluated at nodes, only at material points. It also means that there is no such thing as unique shape functions N a (x), rather we have N a = N a (x, β) and β varies based on where the shape function is being evaluated -here, N a (x p ) = N a (x p , β p ). Hence we also cannot plot the shape function N a for any node a (this is why we chose the discrete visualization in Fig. 1). This choice is necessitated by the specific adaptivity scheme used here (see ( 43) and its anisotropic generalization in Section 2.4 below).

While a node might inherit different β when computing shape functions at different material points, it will be equal for all nodes in the context of the local approximation at a given material point. We also note that shape functions and shape function derivatives at step (n + 1) are calculated based on the previous values β p n . Therefore, we may drop the spatial dependence on the shape functions through β when calculating their gradients in (13), as is also done for the nodal locations in the updated-Lagrangian setting.

Anisotropic adaptivity and bounded shape function support

To eliminate tensile instability in higher dimensions, we need to anisotropically update the support of the max-ent shape functions; i.e., the Pareto optimality between shape function width and information entropy should be updated differently for each dimension based on the local deformation, which is exactly why a tensor-valued locality parameter β was introduced in (3). Assuming the step size between configurations n and n + 1 is sufficiently small, the invariance of the shape function support in (10) implies

N a n (x p n ) = N a n+1 (x p n+1 ) =⇒ (x p n -x a n ) • β p n (x p n -x a ) = (x p n+1 -x a n+1 ) • β p n+1 (x p n+1 -x a n+1 ). ( 44
)
If the effective support size of N a is reasonably small to approximate the local deformation with an affine map, we may exploit the incremental deformation gradient to obtain the approximation

β p n+1 = (F p n→n+1 ) -T β p n (F p n→n+1 ) -1 . (45)
Hence, even if we start with an isotropic shape function support, anisotropy is introduced incrementally through the deformation gradient at a given material point p. We will refer to this scheme as anisotropic adaptivity of the shape function support. Figure 6 compares the evolution of a representative shape function support based on a 2D node set that is stretched uniaxially. It is evident from Figure 6(d) that the effective support size along the direction of stretch is maintained for anisotropic adaptivity of the locality parameter, whereas it becomes localized for constant and isotropic adaptive locality parameters in Figure 6(b) and Figure 6(c), respectively. When local max-ent shape functions span over all nodes in the domain, their global support makes the implementation computationally impractical. However, shape functions decay exponentially with distance from their respective node according to (10), implying that only a few nodes relatively close to a material point are significant to the approximation accuracy at that material point. Hence, it is practical to construct the approximation based only on nodes that lie inside a cut-off region centered around each material point. An analogous remapping of nodes to material points was discussed for the OTM method [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF]. Since the shape functions decay as exp(-β xx a 2 ), the original formulation by [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] suggests the use of a cut-off sphere whose radius is given by R cut = (-log )/β, where 0 < ≤ 1 is a small tolerance, so the shape function support is truncated beyond the cut-off sphere.

In an updated-Lagrangian setting, the cut-off limits must evolve according to the local deformation; otherwise, a loss in connectivity will effectively result in the localization of shape functions. For example, analogous to the isotropic adaptivity scheme (43), the radius of a cut-off sphere centered at a material point p can be updated by

R p cut,n+1 = (J p n→n+1 ) 1/3 R p cut,n . (46)
However, it is not surprising that updating the cut-off radius isotropically as in (46) will result in a severe loss of connectivity during anisotropic deformation and ultimately give rise to tensile instability. Instead, we propose the use of cut-off ellipsoids to improve the connectivity updates. Using ellipsoidal cut-off regions also provides a better approximation when the spatial distribution of nodes is highly anisotropic. We define that, at the n th step, the connectivity of each material point x p n includes all nodes x a n which lie inside the ellipsoid described by (x a n -

x p n ) • M p n (x a n -x p n ) ≤ 1, (47)
where M p n is a d × d matrix with d positive eigenvalues. At the initial configuration, the cut-off ellipsoid can be chosen appropriately based on the tolerance and the average nodal spacing in each direction around the material point. To avoid any significant loss in connectivity, the ellipsoid is updated based on the deformation map at every material point p according to (x a n -

x p n ) • M p n (x a n -x p n ) = (x a n+1 -x p n+1 ) • M p n+1 (x a n+1 -x p n+1 ), (48)
which, under the assumption that the step size and changes in ellipsoid size are sufficiently small, gives the scheme for updating the anisotropic connectivity:

M p n+1 = (F p n→n+1 ) -T M p n (F p n→n+1 ) -1 . ( 49
)
This is analogous to the update of the locality parameter β, see (45). Figure 7 illustrates how an anisotropic ellipsoidal cut-off region better maintains connectivity compared to the isotropic scheme for the case of significant shear deformation.

Treatment of essential boundary conditions and material interfaces

Local max-ent is not interpolatory, i.e., it does not satisfy the strong Kronecker property. However, it does satisfy the weak Kronecker property on the boundary. Precisely, if x lies on the convex hull Ω of the domain, then N a (x) is zero if the node located at x a lies in the interior of the Ω. This implies that interpolation on the boundary depends only on boundary nodes and is independent of the interior nodes. Therefore, essential boundary conditions can be applied without modifications of the shape function support. Note that the max-ent shape functions only satisfy consistency up to first-order and hence can only impose at most affine boundary conditions exactly; however, this error vanishes with h-refinement as the essential boundary conditions are recovered as a converging sequence of piecewise-linear functions. Also, for problems involving non-convex domains, the weak Kronecker property is not satisfied on the non-convex part of the boundary. Several techniques like the visibility criterion [START_REF] Belytschko | Element-free Galerkin methods[END_REF][START_REF] Krysl | Element-free Galerkin method: Convergence of the continuous and discontinuous shape functions[END_REF], diffraction method (Belytschko et al., 1996a;[START_REF] Organ | Continuous meshless approximations for nonconvex bodies by diffraction and transparency[END_REF] and the transparency method [START_REF] Organ | Continuous meshless approximations for nonconvex bodies by diffraction and transparency[END_REF] have been proposed to deal with non-convex domains in element-free Galerkin methods and can be extended to max-ent approximants. However, all of those have limitations for interpolation on non-convex boundaries; see, e.g., Belytschko et al. (1996b) for a review. Alternatively, the max-ent shape function support can be seamlessly made highly local (by increasing the eigenvalues of the locality parameter β) in the proximity of a non-convex boundary and hence minimizing the extraneous non-convex region included in the convex hull of shape function support [START_REF] González | A higher order method based on local maximum entropy approximation[END_REF]; but requires manual control for problems involving complicated geometries.

Discontinuous derivatives like material discontinuities also pose a challenge in meshfree methods. In conventional FE methods, nodes can conveniently be chosen such that a discontinuity like a material interface does not pass through any element. As a result, capturing solutions with discontinuous derivatives while ensuring displacement continuity across the material interface is straightforward. In contrast, without proper treatment, the support of nonlocal meshfree shape functions like max-ent spreads across any discontinuity, which leads to incorrectly smearing out of the interface. The problem is fundamentally the same as the inability to apply essential boundary conditions. Several complicated techniques have been introduced such as additional special shape functions with discontinuous derivatives [START_REF] Krongauz | EFG approximation with discontinuous derivatives[END_REF], including Lagrange multipliers [START_REF] Cordes | Treatment of material discontinuity in the element-free Galerkin method[END_REF], and the transition to finite elements near interfaces [START_REF] Belytschko | A coupled finite element-element-free Galerkin method[END_REF][START_REF] Krongauz | Enforcement of essential boundary conditions in meshless approximations using finite elements[END_REF][START_REF] Huerta | Enrichment and coupling of the finite element and meshless methods[END_REF].

The local max-ent scheme allows for a straightforward solution for capturing discontinuous derivatives. Our approach is inspired by the technique developed by [START_REF] Cordes | Treatment of material discontinuity in the element-free Galerkin method[END_REF] for material interfaces in the elementfree Galerkin method. Consider the example problem of a 2D domain containing two materials separated by an interface, as illustrated in Figure 8. The spatial discretization is chosen such that the interface contains only nodes and no material points. Each material point and node is assigned uniquely to one material, while the interface nodes belong to both and are shared by both regions5 .

As a crucial step, the shape function support and the resulting connectivity of each material point is truncated to include only those nodes that belong to the same material as the material point itself (see Figure 8(b)). This ensures that the kinematic interpolation at any material point does not explicitly depend on nodes from the other side of the interface. Since the interfacial nodes are shared by both regions and because of the weak Kronecker property, displacement continuity at the interface is automatically ensured. As a result, there is no smoothing nor spurious oscillations across the discontinuity, and discontinuous derivatives are captured accurately (Section 3.3 presents a numerical validation example). The simplicity and the fact that no additional computational costs are associated with this scheme is particularly advantageous for higher-dimensional problems involving multiple material interfaces such as polycrystals, high-volume-fraction composites (e.g., cermets or magneto-rheological elastomers), and fluidstructure interactions.

Note that in problems with complex geometry or large deformations, there may be curved material interfaces that give rise to non-convex sub-regions (see e.g. Figure 9). As a result, the convex hull of the nodal connectivity extends beyond the interface. However, as per the proposed scheme, nodes from the other side of the material interface are still excluded and the approximation at the material point remains unaffected by the kinematics of nodes across the interface. Analogous to the case of non-convex boundaries, the extraneous region of the convex hull will diminish with h-refinement or, alternatively, can be minimized by seamlessly increasing the locality of the shape functions near the interfaces.

Benchmark tests

We present a selection of benchmark simulations in three dimensions to study the performance of the proposed enhanced local max-ent scheme in problems involving large deformation. The first benchmark test simulates the inhomogeneous extension of a hyperelastic block to demonstrate the convergence rate and accuracy of the proposed scheme. The second benchmark simulates the torsion of a hyperelastic cube to illustrate its ability to simulate severe deformation and distortion. The third benchmark showcases the stretching of a periodic representative volume element (RVE) of a two-phase laminate composite to demonstrate the capability to accurately capture discontinuous derivatives using truncated shape function supports. The fourth benchmark simulates two spherical inclusions in a high-volume-fraction composite RVE undergoing such severe deformation that the initially distant inclusions get almost into Hertzian contact. The fifth and final benchmark simulates the elastoplastic response of a cylindrical specimen impacting a wall at high speed, classically known as the Taylor anvil test. Importantly, all five benchmarks demonstrate that the proposed adaptive use of anisotropic shape function supports successfully eliminates tensile instability in the updated-Lagrangian framework up to very large deformations. In fact, without the framework introduced here, all benchmark simulations become unstable even for small deformations. Even with the best chosen parameters, the non-adaptive updated-Lagrangian framework could only simulate no more than 10% of the deformations (in all benchmarks shown here) that our new framework can handle.

All simulations follow the same protocol to set up the initial configuration. Once the initial set of nodes is defined, an ad-hoc mesh is created by applying Delaunay triangulation to the initial node set. Next, a material point is inserted at the barycenter of each simplicial element with an initial volume V p 0 equal to the product of simplicial quadrature weight and volume of the containing simplex. The initial value of the locality parameter associated with each material point is set to β p 0 = γh -2 p I, where h p is the average distance between the associated material point and nodes of the surrounding element in the initial mesh6 , and γ is a user input, as discussed in Section 2.3. The mesh creation is a one-time process and the ad-hoc mesh is discarded after the initial material point setup. The initial connectivity of each material point is given by the cutoff ellipsoid associated with M p n = R -2 cut,p I, where R cut,p = (-log )h 2 p /γ, and = 10 -6 . Only one material point is used per tetrahedron in the initial auxiliary mesh. The proposed scheme to overcome tensile instability is, in principle, independent of the number of material points and certainly extensible to higher-order integration (i.e., to larger numbers of integration points per tetrahedron in the auxiliary mesh). We have experimented with using a higher number of material points per tetrahedron and, while minor improvements are seen, we could not draw any general conclusions. Our focus here is on overcoming the tensile instability (which is achieved with only a single integration point per initial tetrahedron in all examples shown). The effects of numerical quadrature were studied by [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] in the context of (isotropic) local maximum-entropy approximants, which also extend to our method.

As a representative example, all hyperelastic benchmark simulations (wherever applicable) use a compressible Neo-Hookean constitutive model to incorporate material non-linearity. The specific strain energy density is given by

W(F) = µ 2 tr F T F J -2/3 -3 + κ 2 (J -1) 2 , (50)
where µ and κ are the shear and bulk moduli, respectively, and J = det F. For simplicity, all subsequent simulations results are based on non-dimensional material and geometry parameters. We note that the presented max-ent and updated-Lagrangian schemes are sufficiently general to apply to, in principle, arbitrary material constitutive laws including elastoplasticity (see benchmark V in Section 3.5); and may make use of the extension in [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF] for general inelastic solids via variational constitutive updates [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF]. Newton-Raphson iteration is used to solve for the nodal displacements at each load step (or time step in case of the fifth, dynamic benchmark in Section 3.5 where explicit time integration is used). After every load/time step, the reference configuration is updated according to the updated-Lagrangian framework. Where applicable, obtained simulation results are compared to those of (total-Lagrangian) FEM overkill simulations (using 10-node tetrahedral elements). Since comparing quantities such as the strain energy entails ambiguity due to the different reference configurations in the total-vs. updated-Lagrangian settings, we use force-displacement curves as a comparison metric, because they are independent of the choice of reference configuration.

Benchmark I: extension of a hyperelastic block

The first benchmark test aims to reproduce the numerical experiment in the original formulation by [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] as a baseline. We simulate a hyperelastic Neo-Hookean cube that is stretched vertically with the top and bottom faces being fixed in all directions -due to symmetry it is sufficient to model one eighth of the block with the top face fixed, two adjacent side faces traction-free and symmetry boundary conditions applied on the rest of the faces, while the bottom face is vertically fixed and free to deform in-plane. The shear and bulk moduli are set to µ = 5 and κ = 10, respectively. The cube is deformed up to 100% extension in 50 load steps (see Video 1) and the total vertical reaction force on the top face is recorded. Figure 10(a) and 10(b) show the final deformed state obtained using local max-ent and FEM. For the purpose of a convergence study, the initial spatial discretization is a uniform grid of nodes with varying degrees of nodal spacing (see Figure 10(d)), and different values of γ to initialize the locality parameter β. Due to highly nonlinear deformations, an overkill FEM simulation with 20 × 20 × 20 elements is used as a proxy for the exact solution. Figure 10(c) illustrates the convergence in the reaction force with refinement of the spatial discretization. Three important observations are made: (i) The error produced by the local max-ent simulations is always smaller than that of FEM, implying improved accuracy, as similarly observed by [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF] for the same benchmark. (ii) The shown errors corresponding to larger values of γ approach the FEM errors, which is corroborated by the fact that larger γ implies higher locality in interpolation, thus tending towards the affine interpolation of FEM. (iii) More importantly, the error produced by the updated-Lagrangian formulation is lower than that obtained from total-Lagrangian calculations with the same γ, which is expected due to the large deformations involved.

Benchmark II: torsion of a hyperelastic cube

We simulate the torsion response of a hyperelastic Neo-Hookean cube with fixed top and bottom faces (this time the entire cube is modeled, see Figure 11(a)-11(d)). The cube is deformed in incremental torsion steps of 0.05 rad ≈ 2.86 • twist angle per step (see Video 2). For added complication, the shear and bulk moduli are set to µ = 5 and κ = 100, respectively. We note that [START_REF] Ortiz | Maximum-entropy meshfree method for compressible and near-incompressible elasticity[END_REF][START_REF] Ortiz | Maximum-entropy meshfree method for incompressible media problems[END_REF] have developed maximum-entropy meshfree methods suited to near-incompressible and incompressible materials. However, to avoid issues related to volumetric locking, we do not consider truly incompressible materials within the scope of this paper. The initial spatial discretization is a uniform grid of nodes with varying degrees of nodal spacing. The initial locality parameter for each material point is set isotropically using γ = 6.0. Figure 11(e) shows the total vertical reaction force on the top face in response to the applied torsion. Remarkably, the enhanced local max-ent scheme is able to simulate up to 200 • degrees of torsion, which is significantly higher than the maximum torsion of 144 • degrees achievable by our comparison FEM simulation. Moreover, local max-ent achieves a similar level of accuracy with the number of material points being only one-sixth the number of elements in FEM. These observations demonstrate the superior capability of the enhanced local max-ent scheme to simulate problems that are highly prone to mesh distortion and tensile instability. We simulate the response of an RVE made of a hyperelastic two-phase laminate composite subjected to uniaxial extension with traction-free lateral faces. As illustrated in Figure 12(a), the RVE consists of a material interface normal to the x-axis separating the two Neo-Hookean materials. To pronounce the discontinuity in derivatives across the interface, the elastic moduli are chosen such that one phase is highly compressible (µ = 1, κ = 1), whereas the other phase is relatively less compressible (µ = 5, κ = 100). Periodic boundary conditions are imposed on the RVE. The initial spatial discretization is a uniform grid of 9 × 9 × 9 nodes. The locality parameter for each material point is initialized using γ = 2.0. The normal strain along the x-axis is computed at each material point and node up to 25% extension. The specific geometry admits an analytical solution within finite elasticity, which predicts homogeneous deformation in each phase with a discontinuity in the deformation gradient at the material interface. A derivation of the (semi-)analytical solution is presented in C. Figure 12 shows the strain distribution when computed with and without shape function support truncation at the interface, respectively. With truncated support as discussed in Section 2.5, the discontinuity is accurately captured by the numerical solution (Figure 12(b)), whereas spurious oscillations are observed when the support truncation is deactivated (Figure 12(c)).

Benchmark III: hyperelastic laminate composite

Benchmark IV: high-volume-fraction composite with spherical inclusion

We consider a hyperelastic composite consisting of a periodic array of spherical inclusions embedded in a matrix. In order to accurately model the inclusions while approaching contact under straining, we simulate a cubic RVE containing two hemispheres of the inclusion phase, each of radius 0.45 with a spacing of 0.1 in between (see Figure 13(a)). We are interested in the stress distribution inside the RVE under uniaxial extension with traction-free lateral faces. Each inclusion (shown in red and green) is modeled by a stiff material (µ = 500, κ = 500) relative to a soft matrix (µ = 1, κ = 1). Periodic boundary conditions are imposed on the RVE. The locality parameter for each material point is initialized using γ = 6.0. The RVE is incrementally stretched as much as possible until the poles of the hemispheres almost come into contact. Figure 13 shows the deformation and reports the normal Cauchy stress distribution at different levels of extension. Remarkably, at 135% extension the inclusions almost come into contact, which is accompanied by a significant and localized increase in the compressive stress near the apparent point of contact (Figure 13(g)). Due to the relatively soft matrix, the stress distribution is akin to the Hertzian solution for elastic bodies in contact. Note that the chosen material models do not capture phenomena like inelasticity or damage near the point of contact; the purpose of this particular benchmark is to demonstrate the numerical capabilities of the method. It becomes apparent that the proposed framework is capable of simulating matrix material being squeezed out quasistatically between the inclusions to allow for contact to happen. By contrast, such simulation are not easily feasible with FEM due to the large mesh distortion and entanglement in the narrow region between the two inclusions.

Benchmark V: Taylor anvil test

As our final example, we depart from hyperelastic constitutive laws and quasistatic simulations. We simulate the elastoplastic dynamic response of a copper specimen impacting a wall at high speed (see Figure 14(a)), classically known as Taylor's anvil test (1948). The specimen is a homogeneous cylinder whose material response is described by von Mises (J 2 ) plasticity with linear isotropic hardening (the mass density is 8930 kg/m 3 , Young's modulus 117 GPa, Poisson's ratio 0.35, yield stress 400 MPa, and plastic modulus 100 MPa). The wall is assumed to be frictionless and oriented normal to the initial velocity of the specimen. A 3D simulation using explicit dynamics is performed. We note that, when using explicit time integration, the proposed method is similar to the Optimal Transportation Meshfree method of [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF] with the exception of the anisotropic shape functions and adaptivity schemes proposed here. The following two scenarios are simulated.

Case (A): a specimen of radius 3.2 mm and length 32.4 mm impacts the wall at a speed of 227 m/s (see Figure 14(a)). The cylinder is discretized into 11, 313 nodes and 52, 608 material points with the help of an initial tetrahedral mesh. A stable time step of 0.01 µs is chosen and the results at time 80 µs are compared to those obtained by FEM [START_REF] Kamoulakis | A simple benchmark for impact[END_REF][START_REF] Zhu | Unified and mixed formulation of the 4-node quadrilateral elements by assumed strain method: Application to thermomechanical problems[END_REF][START_REF] Camacho | Adaptive lagrangian modelling of ballistic penetration of metallic targets[END_REF], EFG [START_REF] Belytschko | Smoothing and accelerated computations in the element free Galerkin method[END_REF], and OTM [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF] in Table 1. Figures14(b) and 14(c) show the time history of the mushroom radius and height of the specimen, respectively (see [START_REF] Belytschko | Smoothing and accelerated computations in the element free Galerkin method[END_REF] for comparison). Figures14(d)-14(f) show the evolution of effective plastic strain at different times during the impact. The results closely match those reported in the literature and obtained using the aforementioned methods.

Case (B): a specimen of radius 3.2 mm and length 12.8 mm impacts the wall at a speed of 750 m/s (see Figure 15(a)). The cylinder is discretized into 5, 483 nodes and 26, 176 material points with the help of an initial tetrahedral mesh. A stable time step of 0.0025 µs is chosen, and the evolution of the effective plastic strain at different times during the impact is shown in Figure 15. While case (A) is a well-known benchmark with experimental validation, (B) serves as a numerical experiment without validation, proposed by [START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF], to push the robustness of the numerical scheme to extremes. In the latter case, due to considerably higher speed at impact, the specimen undergoes severe plastic flow and flattens out. Particularly in the mushroom region, the nodal spacing increases rapidly, which may cause tensile instability. In fact, without the anisotropic adaptive scheme outlined here, tensile instability emerges early on during the simulation at around 2.75 µs.

This last benchmark shows that the proposed method equally applies to inelastic constitutive laws and that the adaptivity scheme ensures the stability of the updated-Lagrangian kernel in scenarios exhibiting severe deformation such as during plastic flow. The benchmark further demonstrates that the adaptivity scheme is not limited to quasistatic settings but can be extended without significant modifications to dynamics.

Conclusions

We have presented a meshfree simulation framework with significantly enhanced stability, particularly for boundary value problems involving severe, finite deformations. We have formulated an updated-Lagrangian scheme based on an incremental deformation map and material point sampling. The scheme utilizes an improved local maximumentropy approximation that is extended to admit anisotropic shape function support. This is achieved by establishing an anisotropic compromise between minimal shape function widths and maximal information entropy via a tensorial locality parameter, in an extension of the original max-ent scheme of [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF]. We have introduced an adaptive scheme for the anisotropic evolution of shape function support and nodal connectivity, which jointly suppress tensile instability up to large deformations, as demonstrated by several benchmark tests. The weak Kronecker property of the original max-ent formulation is retained, which allows for the direct application of essential boundary conditions. We have also presented a scheme that utilizes truncated shape function connectivity to accurately capture discontinuous derivatives across material interfaces as found, e.g., in composite materials. The entire 24 simulation framework has been tested through several hyperelastic and elastoplastic benchmark problems involving severe deformations, and we observed faster convergence and the ability to simulate significantly larger distortions as compared to classical FEM and total-Lagrangian meshfree formulations. We also demonstrated the performance of the methodology in dynamic problems and reported excellent agreement of simulated Taylor impact results with those reported in the literature. The presented approximation scheme as well as the updated-Lagrangian setting are sufficiently general to apply to other problems beyond the scope of the present investigation, including multi-physics problems and incompressible materials, which will be presented in a subsequent study.

A. Derivation of solution (10)-( 12) to constrained optimization problem in (9)

The following derivation follows closely that presented by [START_REF] Arroyo | Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods[END_REF], here modified for the proposed scheme and filling in details. The first-order consistency condition in (2) can be reformulated by substituting the zeroth-order consistency condition as

n n a=1 N a (x)x a -        n n a=1 N a        x = n n a=1 N a (x)(x a -x) = 0.
(51)

The Lagrangian associated with the constrained optimization in (9) istherefore given by

L N β , λ 0 , λ = n n a=1 f (N a , β) + λ 0        n n a=1 N a -1        + λ •        n n a=1 N a (x a -x)        , (52)
where λ 0 ∈ R and λ ∈ R d are Lagrange multipliers, and f (N a , β) = N a xx a 2 β + N a ln N a (x).

(53)

Stationarity of ( 52) with respect to N a (x) for a = 1, . . . , n n yields

∂L ∂N a (x) = x -x a 2 β + 1 + ln N a + λ 0 + λ • (x a -x) = 0 ∀ a = 1, . . . , n n (54)
whose solution reads

N a (x) = exp -x -x a 2 β + λ • (x -x a ) exp(1 + λ 0 ) . ( 55
)
Note that the exponential form ensures the positivity constraint on the shape functions. The Lagrange dual function is given by g(λ 0 , λ) = inf

The function conjugate of f (N a , β) is obtained as

f * (Q a , β) = sup N a (Q a N a -f (N a , β)) . (57)
Stationarity with respect to N a requires

Q a -x -x a 2 β -1 -ln N a = 0 ⇒ N a = exp Q a -x -x a 2 β -1 . (58)
Therefore, the conjugate function f * (Q a , β) is given by Error: Insufficient nodal connectivity 13:

f * (Q a , β) = exp Q a -x -x a 2 β -1 . ( 59
end if 14: end for

Figure 1 :

 1 Figure 1: Local max-ent shape functions N a (x) evaluated at fixed interior point x (viz., the central node) of a 2D node set with different choices of the locality parameter β. That is, the graphs show the value of each node's shape function N a when evaluated at the central node, thus illustrating the contributions of nodal values to the interpolation at the central node. (We note that this visualization is different from that of Arroyo and Ortiz (2006) who plotted N a (x) as a function of x for fixed a). Cases (a,b,c) show the isotropic and (d,e) the anisotropic transition from global to local shape function support with increasing components of β.

Figure 2 :

 2 Figure 2: Local max-ent shape functions N a (x) evaluated at fixed point x: (a) on an edge, (b) on a corner of a 2D node set. As in Figure1, note that we do not plot a given shape function N a as a function of position x but we compute shape functions N a at a given point x.

Figure 4 :

 4 Figure 4: For all shape functions evaluated at a material point initially at x p 0 = 0.5 on a 1D domain Ω 0 = [0, 1] with constant locality parameter β p n = 1, (a) values N a of all a = 1, . . . , n n shape functions, plotted for different stretch ratios λ = |Ω n |; (b) associated internal nodal force kernels (F p n ∇ n N a V p 0 ) of all shape functions evaluated for different stretches λ; (c) internal nodal force kernels of selected nodes vs. strain ε = λ -1; (d) internal nodal force on selected nodes vs. strain ε = λ -1.

  Figure 5: (a) Internal nodal force kernels and (b) internal nodal forces of selected nodes vs. strain ε = λ -1, with an adaptively modified locality parameter β p n ∝ h -2 n where h n is the (uniform) nodal spacing and with an initial β p 0 = 1.

Figure 6 : 0 =Figure 7 :

 607 Figure 6: Shape function support of a material point initially located at x p 0 = (0.5, 0.5) T under homogeneous anisotropic deformation described by the total deformation gradient F 0→1 = 3 0 0 1 .

  Figure 8: Domains of nodal connectivity of four representative material points near an interface separating two different materials (indicated by the red and blue regions). Material points and nodes are denoted by open and filled dots, respectively. The hatched area around each material point indicates the cut-off region that encloses all nodes whose shape functions are contributing to the material point.

Figure 9 :

 9 Figure 9: Shape function support for a material point (blue open dot) in the proximity of a non-convex material interface. The convex hull (indicated by black solid line) extends beyond the interface. However, only nodes in the blue region (indicated by blue filled dots) are considered during approximation construction. Nodes from the red region (indicated by red filled dots) that lie inside the convex hull are excluded from the approximation.

  (a) Deformation at 100% extension using local max-ent (b) Deformation at 100% extension using FEM

  Figure 10: Benchmark I: extension of a hyperelastic block

  Figure 11: Benchmark II: torsion of a hyperelastic cube 20

  Figure 12: Benchmark III: hyperelastic laminate composite in a periodic RVE

Figure 13 :

 13 Figure 13: Benchmark IV: hyperelastic high-volume-fraction composite with periodic spherical inclusions. (a) Initial state of the RVE with two hemispherical inclusions (both nodes and material points are shown in red and green) embedded in a matrix (both nodes and material points in blue). (b, c), (d,e), (f,g) show the deformed configuration and the stress distribution at 45%, 90%, and 135% extension, respectively. The stresses shown in (c,e, g) are the Cauchy stress component along the x-direction interpolated at the central plane of the RVE denoted by the red rectangles in (b, d, f). Due to the large deformations involved, the stress plots use different scales at each step for improved visibility.

  Figure 14: Benchmark V-(A): Taylor anvil test, showing (a) the initial specimen, (b) the time history of the mushroom radius, (c) the time history of the specimen height, and (d,e,f) cut-section views of the specimen near the wall with effective plastic strains evaluated at material points (nodes are shown as gray dots).

  Figure 15: Benchmark V-(B): Taylor anvil test, showing (a) the initial specimen, and (b,c,d) the effective plastic strains the deformed sample, evaluated at the material points (nodes are shown as gray dots).

  N a (x)≥0, a=1,...,n n L N β .

) 0 = γ h p 0 - 2 I 0 = c h p 0 - 2 I 8 :

 0020028 Rewriting the Lagrange dual function (56) in terms of the function conjugate[START_REF] Bonet | Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations[END_REF] givesg(λ 0 , λ) = -λ 0 -n n a=1 f * (-λ 0 -λ.(x ax), β) . 0 -λ.(x ax) -xx a 2 β --x a 2 β + λ • (x -x a ) = Z(x, λ),(62)so that the Lagrange dual function reduces tog(λ * 0 , λ) = -λ * 0 -1 =ln Z(x, λ). (63)Finally, maximizing the reduced Lagrange dual function with respect to λ givesλ * = arg min λ∈R d ln Z(x, λ).(64)B. Pseudo-code for hyperelastic boundary value problemThe following algorithms summarize the numerical realization of the presented updated-Lagrangian meshfree framework for the example of a hyperelastic boundary value problem at finite strains.Algorithm 1 Example algorithm for a hyperelasticity simulation using anisotropic local max-ent 1: Initialize the spatial discretization (Algorithm 2) 2: for each material point p = 1, . . . , n p do 3:Set initial configuration as undeformed: F p 0 = I 4: end for 5: for n = 0, . . . , # of quasistatic load or time steps do increments ∆U n by assembly and iterative solution (Algorithm 4Initialize the spatial discretization Input: γ, c1: Create a simplicial mesh based on a given geometry 2: Initialize n n nodes at mesh vertices x a 0 (a = 1, . . . , n n ) 3: Initialize n p material points at the simplex quadrature points x p 0 (p = 1, . . . , n p ) 4: Initialize all material point densities: ρ p 0 = ρ 0 (x p 0 ) 5: Initialize all material point volumes: V p 0 = quadrature weight × simplex volume 6: Initialize the locality parameter for each material point: β p , where h p 0 is the mean distance from x p 0 to nodes of the surrounding elements in the mesh 7: Initialize cut-off region (ellipsoid) for each material point: M p Discard the mesh Algorithm 3 : Compute/update shape functions 1: for each material point p = 1, . . . , n p do 2: Find nodal connectivity based on cut-off region S p n = {x a n : (x a n -x p n ) • M p n (x a n -x p n ) ≤ 1} 3: if size(S p n ) ≥ (d + 1) and x p n ∈ conv(S p n )

One may also relax the first constraint and admit negative values of shape functions, see e.g.(Bompadre et al.,

2012), which however is not the focus here.2 Here and in the following, we use classical tensor and index notation common to continuum mechanics. Specifically, dots denote inner products of tensors of equal order, so v • u = v i u i and T • S = T i j S i j for vectors v, u ∈ R d and second-order tensors T, S ∈ GL(d). For tensors as linear mappings of vectors we write [Tv] j = T i j v j and anlogously [ST] i j = S ik T k j . Outer products are represented by [v ⊗ u] i j = v i u j .

Note that in subsequent benchmark examples the energy density is chosen to be quasiconvex, so that unique solutions may be found. However, the same numerical procedure can, in principle, be extended to non-convex potentials, in which the formation of microstructure is implied.

Note that[START_REF] Li | Optimal transportation meshfree approximation schemes for fluid and plastic flows[END_REF] focused on dynamic problems solved by explicit updates, whereas we are concerned with a quasistatic formulation. Therefore, the terminology of time vs. load steps is, in principle, interchangeable.

Note that the assignment of nodes to the different materials occurs only for purpose of shape function definitions, whereas the materials' constitutive models are evaluated only at the material points. This justifies the unique assignment of material points and the non-unique assignment of nodes to the materials

We note that all (structured and unstructured) node sets used in the following benchmark simulations have sufficient regularity for an isotropic initial β. In case of strongly anisotropic initial configurations, it is straightforward to adapt the definition of β p 0 accordingly.

Table 1: Benchmark V-(A): Comparison of results from Taylor anvil test

Acknowledgments

DMK and SK acknowledge support from the National Science Foundation through CAREER Award CMMI-1254424. KD acknowledges support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (Grant agreement no. 636903).

Algorithm 4 : Assembly and solver Input: ∆U n (initial guess = solution of previous step) 1: do 2:

Initialize global force vector and tangent matrix:

for each material point p = 1, . . . , n p do 4:

Compute incremental deformation gradient:

Compute total deformation gradient:

for all nodes a in S p n do 7:

Compute and assemble internal nodal forces:

for all nodes b in S p n do 9:

Compute and assemble tangent matrix: 10:

end for 12:

end for 13:

end for 14:

for each node a = 1, . . . , n n do 15:

Assemble external forces:

16:

for each node b = 1, . . . , n n do 17:

Assemble external tangent matrix:

end for 19:

end for 20:

Apply essential boundary conditions to force vector and tangent matrix 21:

Newton-Raphson step:

1: for each material point p = 1, . . . , n p do 2:

Compute incremental deformation gradient:

Compute total deformation gradient:

Update material point density:

Update material point volume:

Update locality parameter:

Update cut-off region:

Analytical solution of the finite-elasticity problem in benchmark III

Consider a hyperelastic cube of side length 1 composed of two perfectly bonded Neo-Hookean materials denoted by symbols A and B, which are arranged into a laminate whose material interface is given by the plane x = 0.5 in the initial configuration (see Figure 12). µ i and κ i denote the shear and bulk moduli, respectively, of material i. For the prescribed boundary conditions, the overall deformation gradient of the laminate is given by (66) where λ A 1 and λ B 1 are the stretch ratios for each material along the x-axis. Since the initial thickness of each laminate phase is the same, the overall stretch ratio λ 1 is related to λ A 1 and λ B 1 via

The overall strain energy density Π = W of the laminate block, using (50), is thus given by

) The principle of minimum potential energy requires that for a given applied λ 3

where the objective function being minimized is convex and admits a unique solution. Hence, the optimization problem in (69) can be easily solved numerically to yield the required stretch ratios.