Cell-based Maximum Entropy Approximants for Three Dimensional Domains: Application in Large Strain Elastodynamics using the Meshless Total Lagrangian Explicit Dynamics Method

TL;DR

This paper introduces Cell-based Maximum Entropy (CME) approximants for 3D domains, enhancing meshless elastodynamics simulations by improving accuracy, stability, and boundary condition enforcement in large strain problems.

Contribution

The paper develops a novel meshfree approximation method combining maximum entropy properties with compact support, specifically tailored for 3D elastodynamics applications.

Findings

CME provides accurate boundary condition imposition without complex constraints.

The method maintains numerical stability for large time steps in explicit dynamics.

CME outperforms traditional meshless methods in large-scale 3D elastodynamics simulations.

Abstract

We present the Cell-based Maximum Entropy (CME) approximants in E3 space by constructing the smooth approximation distance function to polyhedral surfaces. CME is a meshfree approximation method combining the properties of the Maximum Entropy approximants and the compact support of element-based interpolants. The method is evaluated in problems of large strain elastodynamics for three-dimensional (3D) continua using the well-established Meshless Total Lagrangian Explicit Dynamics (MTLED) method. The accuracy and efficiency of the method is assessed in several numerical examples in terms of computational time, accuracy in boundary conditions imposition, and strain energy density error. Due to the smoothness of CME basis functions, the numerical stability in explicit time integration is preserved for large time step. The challenging task of essential boundary conditions imposition in non-interpolating meshless methods (e.g., Moving Least Squares) is eliminated in CME due to the weak Kronecker-delta property. The essential boundary conditions are imposed directly, similar to the Finite Element Method. CME is proven a valuable alternative to other meshless and element-based methods for large-scale elastodynamics in 3D.

Figures

Tables

Table 1. Table 1: Normalized Root Mean Square Error ( N ​ R ​ M ​ S ​ E 𝑁 𝑅 𝑀 𝑆 𝐸 NRMSE ) for a cube under unconstrained compression.

Nodes	Quadrature points	$NRMSE$
152	1896	1.75E-3
4594	92876	3.85E-4

Table 2. Table 2: Simulation characteristics (critical explicit integration time step - Δ ​ t c ​ r ​ i ​ t Δ subscript 𝑡 𝑐 𝑟 𝑖 𝑡 \Delta t_{crit} , number of execution steps - N e ​ x ​ e subscript 𝑁 𝑒 𝑥 𝑒 N_{exe} , execution time - t e ​ x ​ e subscript 𝑡 𝑒 𝑥 𝑒 t_{exe} ) for cube under 50 % percent 50 50\% uniaxial extension & \& compression for l ​ v ​ l 1 𝑙 𝑣 subscript 𝑙 1 lvl_{1} - l ​ v ​ l 6 𝑙 𝑣 subscript 𝑙 6 lvl_{6} discretization.

Approx.	$s$ $†$	$a$ $‡$	$\Delta t_{crit}$ (ms)	$N_{exe}$	$t_{exe}$ (s)
Extension
CME	2	0.2	2.173 - 0.260	2963 - 10455	1.4 - 3173.5
		1.6	1.972 - 0.243	2814 - 12707	1.6 – 3906.9
		2.0	1.918 - 0.238	3341- 13122	1.4 - 3911.0
	4	0.2	1.527 - 0.205	3396 - 16582	1.9 - 4652.4
		1.6	1.230 - 0.176	3529 - 20900	1.8 - 5711.7
		2.0	1.216 - 0.174	3509 - 21167	1.9 - 6070.6
MMLS-SEBCIEM	1.6	0.2	0.902 - 0.192	3208 - 16827	1.4 - 4603.5
Compression
CME	2	0.2	2.173 - 0.260	3172 - 10442	1.7 - 3184.9
		1.6	1.972 - 0.243	3093 - 17288	1.5 - 4630.3
		2.0	1.918 - 0.238	3117 - 16501	1.5 - 4414.7
	4	0.2	1.527 - 0.205	3279 - 31402	1.7 - 7231.9
		1.6	1.230 - 0.176	3664 - 24885	2.0 - 6297.4
		2.0	1.216 - 0.174	3735 - 18563	1.7 - 5557.1
MMLS-SEBCIEM	1.6	0.2	0.902 - 0.192	3975 - 18315	1.8 - 5710.6

Table 3. Table 3: Simulation characteristics (critical explicit integration time step - Δ ​ t c ​ r ​ i ​ t Δ subscript 𝑡 𝑐 𝑟 𝑖 𝑡 \Delta t_{crit} , number of execution steps - N e ​ x ​ e subscript 𝑁 𝑒 𝑥 𝑒 N_{exe} , execution time - t e ​ x ​ e subscript 𝑡 𝑒 𝑥 𝑒 t_{exe} , Essential boundary conditions (EBC) values) for cardiac bi-ventricular model in indentation.

Approx.	$\Delta t_{crit}$ (ms)	$N_{exe}$	$t_{exe}$ (h)	Displaced EBC$†$ (m)	Fixed EBC$†$ (m)
CME	0.636	36073	1.06	-0.01 $\pm$ 1.0E-17	0.00 $\pm$ 0.00
MMLS-SEBCIEM	0.455	49692	1.53	-0.01 $\pm$ 3.3E-15	0.00 $\pm$ 1.7E-15

Equations

$$u^h(\mathbf{x})=\sum _{a=1}^n{\varphi }_a(\mathbf{x})u_a$$

$${\varphi }_a(\mathbf{x})\ge 0,\sum _{a=1}^n{\varphi }_a(\mathbf{x})=1,\sum _{a=1}^n{\varphi }_a(\mathbf{x}){\mathbf{x}}_a=\mathbf{x}$$

$$H(\phi|w)=-\sum_{a=1}^{n}\phi_{a}(\bm{x})\ln\bigg{(}\frac{\phi_{a}(\bm{x})}{w_{a}(\bm{x})}\bigg{)}$$

$$(ME{)}_w\text{For fixed }\mathbf{x}\text{, maximize }$$

$${\varphi }_a(\mathbf{x})\ge 0(a=1,\dots ,n),$$

$$\sum _{a=1}^n{\varphi }_a(\mathbf{x})=1,$$

$$\sum _{a=1}^n{\varphi }_a(\mathbf{x}){\mathbf{x}}_a=\mathbf{x}$$

$${\varphi }_a(\mathbf{x})=\frac{w_a(\mathbf{x})\text{ exp}[\mathbf{\lambda }(\mathbf{x})\cdot (\mathbf{x}-{\mathbf{x}}_a)]}{Z(\mathbf{x},\mathbf{\lambda }(\mathbf{x}))}$$

$$Z(\mathbf{x},\mathbf{\lambda }(\mathbf{x}))=\sum _{b=1}^n{w}_b(\mathbf{x})\text{ exp}[\mathbf{\lambda }(\mathbf{x})\cdot (\mathbf{x}-{\mathbf{x}}_b)]$$

$${\mathbf{\lambda }}^{\ast }(x)={\text{arg min}}_{\mathbf{\lambda }\in E^d}\ln Z(\mathbf{x},\mathbf{\lambda })$$

$$\nabla {\varphi }_a^{\ast }=g_a^{\ast }\nabla w_a+{\varphi }_a^{\ast }\left[D{\mathbf{\lambda }}^{\ast }\cdot (\mathbf{x}-{\mathbf{x}}_a)-\sum _{b=1}^n{g}_b^{\ast }\nabla w_b\right]$$

$$D\bm{\lambda}^{*}=-\bigg{(}\sum_{a=1}^{n}g_{a}^{*}\nabla w_{a}\otimes(\bm{x}-\bm{x}_{a})+\bm{I}\bigg{)}(\bm{J}^{*})^{-1},$$

$$f\equiv f(\mathbf{x}):=\frac{(x-x_1)(y_2-y_1)-(y-y_1)(x_2-x_1)}{L}$$

$$t\equiv t(\mathbf{x}):=\frac{1}{L}\left[{\frac{L}{2}}^2-\parallel \mathbf{x}-{\mathbf{x}}_c\parallel \right]$$

$$\rho\equiv\rho(\bm{x}):=\sqrt{f^{2}+\frac{1}{4}\bigg{(}\sqrt{t^{2}+f^{4}}-t\bigg{)}^{2}}$$

$$d_a(l_1,l_2,\dots ,l_n):={\rho }_1\sim {\rho }_2\sim \dots \sim {\rho }_n=\frac{1}{\sqrt[m]{\frac{1}{{\rho }_1^m}+\frac{1}{{\rho }_2^m}+\dots +\frac{1}{{\rho }_n^m}}}$$

$$w_a(\mathbf{x})=\frac{d_a^s(\mathbf{x})}{{\sum }_{b\in N_x}{d}_b^s(\mathbf{x})}$$

$$p_1{\wedge }_f{p}_2\equiv p_1+p_2-\sqrt{p_1^2+p_2^2+\alpha f^k}$$

$$t\equiv (p_1{\wedge }_f{p}_2){\wedge }_f{p}_3\equiv p_1+p_2+p_3-G-\sqrt{(p_1+p_2-G{)}^2+p_3^2+\alpha f^k}$$

$$\nabla d_a(\mathbf{x})\equiv \nabla d_a=\frac{{\sum }_{i=1}^{n_a}{\rho }_i^{-(m+1)}\nabla {\rho }_i}{({\sum }_{i=1}^{n_a}{\rho }_i^{-m}{)}^{\frac{m+1}{m}}}$$

$$\nabla\rho=\frac{1}{\rho}\Bigg{(}f\nabla f+\left[\frac{\sqrt{t^{2}+f^{4}}-t}{4}\right]\left[\frac{t\nabla t+2f^{3}\nabla f}{\sqrt{t^{2}+f^{4}}}\right]\Bigg{)}$$

$$\nabla t$$

$$\nabla G$$

$$W=\frac{\mu }{2}(I_1-3)-\mu \ln J+\frac{\lambda }{2}(\ln J{)}^2$$

$$NRMSE=\frac{\frac{1}{N}{\sum }_{i=1}^N(u_{z_i}^h-u_{z_i}^{an}{)}^2}{u_{z_{max}}^{an}-u_{z_{min}}^{an}}$$

$$SR{E}_W=\frac{{\bar{W}}_{lv{l}_i}-{\bar{W}}_{lv{l}_6}}{{\bar{W}}_{lv{l}_6}}$$

Full text

K****eywords Meshless  $\cdot$ Cell-based Maximum Entropy  $\cdot$ Large strain elastodynamics  $\cdot$ Explicit time integration  $\cdot$ Weak Kronecker-delta  $\cdot$ Essential boundary conditions imposition

1 Introduction

Meshfree Methods (MMs) [1, 2] have been proposed as an alternative to the widely used Finite Element Method (FEM) for applications in solid mechanics, due to their attractive property to alleviate the mesh requirement for the approximation of an unknown field function. In mesh-based techniques, such as FEM, the approximation accuracy is deteriorated in problems that involve large deformations due to the mesh distortion. MMs are more suited for simulations involving large deformations. In MMs, the spatial domain is discretized by a set of nodes arbitrarily distributed without any interconnectivity. Since the introduction of Smoothed Particle Hydrodynamics (SPH) [3, 4], MMs proliferated with developments such as the Meshless Collocation Methods (MCMs) [5, 6, 7], the Meshless Local Petrov-Galerkin (MLPG) [8] and the Element-Free Galerkin (EFG) [9] methods. A detailed overview of meshfree methods and their advantages can be found in [1, 10].

In MCMs, the strong form solution of a PDE system is approximated using nodal collocation on the field nodes. Such methods demonstrate high efficiency, easy implementation, and direct Essential Boundary Conditions (EBC) imposition. However, MCMs often suffer from low accuracy when natural boundary conditions are involved due to the lower precision in the approximation of high-order derivatives [11]. In weak form MMs, such as the MLPG and EFG methods, the order of derivation is reduced. Natural boundary conditions are satisfied in a straightforward manner. Weak form MMs have demonstrated high numerical stability and accuracy in large strain problems [12, 13].

The Meshless Total Lagrangian Explicit Dynamics (MTLED) [14, 15, 16] is an efficient variant of the EFG method that was introduced to deal with large deformations. In MTLED, all calculations refer to the initial configuration of the continuum (total-Lagrangian formulation). MTLED implements a central difference scheme for explicit time integration, where basis functions and spatial derivatives are computed once, during the pre-processing stage. In addition, the explicit time integration demonstrates advantages over implicit integration schemes, such as (i) straightforward consideration of geometrical and material nonlinearity; (ii) easy implementation; (iii) suitability for massive parallelism [14]. Originally, basis functions approximation in MTLED was performed using the Moving Least Squares (MLS) scheme [14]. An improved version of the MLS with advanced approximation capability named as the Modified MLS (MMLS) has been introduced in [17]. The drawback of the approximation schemes of the MLS group is that the approximation of the field function does not possess the Kronecker-delta property and the imposition of Essential Boundary Conditions (EBC) requires special treatment.

Methods that modify the weak form to impose EBC (e.g., Lagrange multipliers, penalty methods) are not suitable for explicit time integration. Joldes et al. [18] have proposed the EBC imposition in explicit meshless (EBCIEM) method for EBC imposition in MTLED. In EBCIEM, a displacement correction applies by distributing force on EBC nodes through Finite Element shape functions. EBCIEM applies exact imposition of EBC up to machine precision. However, the need for a Finite Element layer on the EBC boundary imposes restrictions regarding the applicability of the method to arbitrary geometries. A simplified version (SEBCIEM) has been also proposed [18], where the distributed forces are lumped on the EBC nodes.

The EBC correction requirement can be alleviated if basis functions that possess the Kronecker-delta property are used. Approximation schemes of the Maximum Entropy (MaxEnt) group possess a weak Kronecker-delta property, which allows imposing EBC directly as in FEM. Different MaxEnt variations have been proposed up to date. Sukumar in [19] used the maximization of information entropy to formulate meshfree interpolants on polygonal elements. Aroyo and Ortiz in [20] used a Pareto compromise between locality of approximation and maximization of information entropy to create the Local Maximum Entropy (LME) approximants. LME approximants use generalized barycentric coordinates based on Jayne’s principle of maximum entropy [21] and provide a seamless transition between non-local approximation and simplicial interpolation on a Delaunay triangulation (linear interpolants in the context of FEM).

LME approximants have already found a large number of applications [22, 23, 24, 25, 26]. LME approximants are smooth and nonnegative approximants with local support that possess the weak Kronecker-delta property in the sense that basis functions on the boundary are not influenced by internal nodes. The LME support domain "locality" is controlled by a non-dimensional parameter $\gamma$. A variational approach to adapt $\gamma$ in the context of LME has been proposed in [27]. Cell-based Maximum Entropy (CME) approximants have been proposed as an alternative to LME in the Galerkin framework [28]. CME capitalizes on the background mesh connectivity only to define compact support for basis functions on $N$-rings ($N\geq 1$) of connected cells while the basis functions values are computed without using the connectivity information. Up to date, CME approximants have been implemented to approximate field functions in $E^{2}$.

The purpose of this work is to extend the CME approximants in $E^{3}$ and use them in the context of three-dimensional (3D) large-strain elastodynamic problems. Due to the previously described advantages of the explicit time integration over the implicit time integration for the solution of large-strain problems, we introduce the CME approximants in the MTLED method. While MTLED is an efficient method to address large deformation problems, the limitation of the special EBC treatment in MMLS can be restrictive in some scenarios. Our objective is to simplify the MTLED method and enable direct EBC imposition by replacing the MMLS with the CME. The rest of this paper is structured as follows. In Section 2, we provide a compact and cogent formulation of the Cell-based Maximum Entropy approximation in $E^{2}$, based on the minimum relative entropy framework. In Section 3, we extend the Cell-based Maximum Entropy approximation in $E^{3}$. In Section 4, we present benchmark examples, illustrating the accuracy and robustness of the proposed scheme in simple and complex geometries. Finally, in Section 5, we discuss our conclusions.

2 Cell-based Maximum Entropy Approximants Formulation In E2

CME approximants belong to the class of convex approximation schemes. They have compact support which leads to linear algebraic systems with small bandwidth and hence reduced memory requirements. CME are more time efficient compared to approximants with dense-nodal supports. Minimal supports are constructed for triangulated domains in $E^{2}$ by replacing Gaussian prior weights functions (commonly used in LME) with smooth distance approximation functions using the R-functions technology [29]. The smooth distance approximation of a point $\bm{x}\in\Omega^{2}$ to the boundary $\partial\Omega^{N_{R}}$ of the $N_{R}$ support ring, $N_{R}=1,2,...,n\text{ and }\partial\Omega^{N_{R}}\setminus\partial\Omega$, is computed to derive the prior weight functions (Figure 1). Following, we give a brief overview of the minimum relative entropy framework (Subsection 2.1) along with the nodal prior weight construction using R-functions (Subsection 2.2). For a detailed description of the CME formulation we refer the reader to Millán et al. [28].

2.1 Minimum Relative Entropy Framework

For a scalar-valued function $u(\bm{x})$ and a set of unstructured points $\{\bm{x}_{a}\}_{a=1}^{n}\subset E^{d}$, the MaxEnt approximation $u^{h}(\bm{x})$ is given by

[TABLE]

where $u_{a}$ are nodal coefficients, and $\{\phi_{a}(\bm{x})\}_{a=1}^{n}$ are nonnegative basis functions that fulfill the zeroth and first order reproducing conditions:

[TABLE]

The basis functions $\{\phi_{a}(\bm{x})\}_{a=1}^{n}$ are defined from an optimization problem that is established at each evaluation point $\bm{x}$ for the linear constraints given by Equation (2). Due to the nonnegative and partition of unity (zeroth order reproducing condition) properties, the basis functions can be interpreted as a discrete probability distribution. The informational entropy provides a canonical measure to the uncertainty associated with a discrete probability distribution. The least-biased approximation scheme that is consistent with the linear constraints is provided by the principle of maximum entropy. The formulation of MaxEnt approximants is derived by maximizing the Shannon-Jaynes entropy measure [21] when nodal prior weights are used:

[TABLE]

where $D(\phi|w)=-H(\phi|w)\geq 0$ is the relative entropy measure and the corresponding variational principle is given by the principle of minimum relative (cross) entropy [30]. The MaxEnt basis functions optimization problem can be stated in the variational formulation:

[TABLE]

Duality methods can be employed to solve the convex optimization problem in Equation (4) efficiently and robustly [30, 31]. The basis functions are then:

[TABLE]

where

[TABLE]

is the partition function and $\bm{\lambda}$ is the Lagrange multiplier vector. The Lagrange multiplier vector is obtained solving the unconstrained convex optimization problem:

[TABLE]

where $\bm{\lambda}^{*}$ is the converged solution of $\bm{\lambda}$ at $\bm{x}$. The resulting basis functions from the solution of $(ME)_{w}$ are noninterpolating, except on the boundary of the nodal set’s convex hull where the weak Kronecker-delta property holds (Figure 2). The weak Kronecker-delta property provides a strong advantage over other approximant types, such as the MLS, where special treatments are required for the imposition of EBC [18]. Moreover, the basis functions inherit the nodal prior weight functions smoothness [32, 31]. Various nodal prior weight functions, such as Gaussian weight function [20, 22, 33], quartic polynomial weight function[34, 35, 36], level set based nodal weight function [37, 38], exponential nodal weight function [39, 40], and approximate distance function to planar curves [28] have been used to construct maximum entropy approximants with specific desired properties. The expression of the gradient for the MaxEnt basis functions for a node $a$ evaluated on a point $\bm{x}\in E^{d}$ is given by:

[TABLE]

where the superscript $*$ on any function indicates evaluation at $\bm{\lambda}^{*}(\bm{x})$, $w_{a}$ is the prior weight function for node $a$, $g_{a}^{*}:=\phi_{a}^{*}/w_{a}$ , and $D\bm{\lambda}^{*}$ is given by:

[TABLE]

where $\bm{I}$ is the identity matrix, $\bm{J}^{*}:=\partial r/\partial\bm{\lambda}$, and $r:=\sum_{a=1}^{n}\phi_{a}(\bm{x}-\bm{x}_{a})$. A detailed description of the derivation of $\nabla\phi_{a}^{*}$ is given in Refs. [28, 34].

2.2 Nodal Prior Weights Construction For CME In E2

To satisfy the required properties of CME (i.e., smoothness, compact support, unimodality), smooth approximations of the distance function to planar curves [41] are considered in the construction of prior weight functions. The theory of R-functions [29] is used to approximate the distance function of each node in the nodal support of a basis function to the boundary polygon of the $N_{R}$-ring support domain (Figure 1). Any real-valued function $F(\omega_{1},\omega_{2},\ldots,\omega_{q})$, where $\omega_{i}(\bm{x}):E^{2}\rightarrow E$, is considered a R-function if its sign is determined only by the sign of its arguments $\omega_{i}$ and not their magnitude. Similarly to logical functions, R-functions can be written as a composition of elementary R-functions using operations such as negation, disjunction, conjunction, and equivalence [29, 41].

Using the composition property of R-functions, the approximate distance function of any given node $a$ from the boundary $\partial\Omega_{a}^{N_{R}}$ of its polygonal support domain $\Omega_{a}^{N_{R}}$ is expressed as the composition of elementary R-functions corresponding to the piecewise linear segments of the polygonal curve $\partial\Omega_{a}^{N_{R}}$. For each line segment belonging to $\partial\Omega_{a}^{N_{R}}$ with endpoints $\bm{x}_{1}\equiv(x_{1},y_{1})$ and $\bm{x}_{2}\equiv(x_{2},y_{2})$, the signed distance function from a point $\bm{x}$ to the line passing from $\bm{x}_{1}$ and $\bm{x}_{2}$ is defined by:

[TABLE]

where $L=\parallel\bm{x}_{2}-\bm{x}_{1}\parallel$ is the length of the line segment. In addition, the line segment defines a disk of radius $L/2$. This disk can be expressed by the trimming function:

[TABLE]

where $t$ is normalized to first order and $t\geq 0$ defines the disk with center $\bm{x}_{c}=(\bm{x}_{1}+\bm{x}_{2})/2$. From Equations (10), (11) the signed distance function from a point $\bm{x}$ to the line segment with endpoints $\bm{x}_{1}$ and $\bm{x}_{2}$ is given by the first order normalized function $\rho(\bm{x})$:

[TABLE]

where $\rho(\bm{x})$ is differentiable for any point that does not lay on the line segment. For $n$ line segments $l$ s.t. $\partial\Omega_{a}^{N_{R}}=l_{1}\bigcup l_{2}\bigcup\ldots\bigcup l_{n}$ the normalized approximation, up to order $m$, of the distance function $d_{a}$ for the node $a$ is given by the R-equivalence formula [41]:

[TABLE]

The R-equivalence formula has the desirable properties of preserving the normalization up to the $m^{th}$ order and being associative, while the R-conjunction is normalized up to the $(m-1)^{th}$ order and it is not associative. Finally, the prior weight function $w_{a}(\bm{x})$ for node $a$ is obtained by:

[TABLE]

where $N_{x}$ are the indices of the $N_{R}$-ring nodal neighbors of the point $\bm{x}$, $s\geq 2$ is a smoothness modulation factor and $w_{a}(\bm{x})$ is normalized s.t. $0\leq w_{a}(\bm{x})\leq 1$. It should be noted that the R-equivalence formula in Equation (13) does not preserve the normalization of the distance function on the joining vertices of the piecewise linear segments and hence introduces undesired bulging effect at these points. Increasing the smoothness modulation factor $s$, the bulging effect is reduced.

3 Extension Of Cell-based Maximum Entropy Approximants In E3

MaxEnt approximants (Equation (5)) are naturally generalized to $E^{d}$, with $d\geq 3$. However, CME approximants have been limited to $E^{2}$ due to the construction of nodal prior weight functions as smooth approximations to the distance function of 2D polygonal curves. In this section we extend the CME approximants to $E^{3}$ by constructing nodal prior weight functions as smooth approximations to the distance functions of polyhedral surfaces following the development of Biswas and Shapiro [41].

3.1 Nodal Prior Weights Construction For CME In E3

We consider the discretization $\Omega^{h}$ of a finite domain $\Omega\in E^{3}$ in tetrahedra $\tau$. The $N_{R}$ ring support domain for any vertex $a$ of $\tau\in\Omega^{h}$ is defined by the union of the piecewise triangular faces belonging to the $N_{R}$ ring of attached tetrahedra to the vertex $a$ (Figure 1). R-functions are employed to construct the approximation distance field to the boundary of the $N_{R}$ ring support domain by joining the fields of the individual boundary triangular faces.

The distance field approximation to any triangle is constructed considering the intersection of the triangle’s carrier plane $f(\bm{x})=0$ and a trim volume $t$ defined by the planes $p(\bm{x})_{i}$, $i=1,2,3$ that pass through the edges of the triangle and are orthogonal to the carrier plane $f(\bm{x})$ (Figure 3). The trim volume $t$ can be constructed joining the planes $p_{i}$ according to the $R_{f}$-conjunction formula:

[TABLE]

where $\alpha$, $k$ are constants controlling the shape of the R-function’s zero set. For points in $E^{3}$ that do not lay on the carrier plane the Equation (15) is $C^{\infty}$, while for points on the carrier plane ($f=0$) it reduces to the standard $R$-conjunction formula. Therefore, $t$ is given by:

[TABLE]

where $G:=\sqrt{p_{1}^{2}+p_{2}^{2}+\alpha f^{k}}$. The function $t$ defines a smooth trimming volume tangent to the three edges of the triangle. The coefficient $\alpha$ controls the smoothness of the trimming volume, where for large $\alpha$ the trimming volume flattens out and for $\alpha=0$ it becomes identical to the unbounded prism constructed by the combination of the planes $p_{i}$ (Figure 3). While the unbounded prism trimming volume ($\alpha=0$) does not maintain differentiability on the triangle edges, the smooth volume ($\alpha>0$) is differentiable at any point except the triangle vertices where the bulging effect arises. Substituting $f$ and $t$ in Equation (12), the normalized distance function from any point $\bm{x}$ to the triangle is approximated and the nodal prior weights $w_{a}(\bm{x})$ are obtained by Equations (13), (14), similarly to the $E^{2}$ case.

3.2 Derivation Of The CME Approximants Gradient In E3

To compute the gradient of the CME approximants $\nabla\phi_{a}(\bm{x})$ (Equation (8)), we should first compute the gradient of the prior function $\nabla w_{a}(\bm{x})$. Following the abbreviation in [28], to derive $\nabla w_{a}(\bm{x})$, for a node $a$ at a point $\bm{x}$, we first derive the gradient of the normalized approximation to the distance function:

[TABLE]

where $\rho_{i}$ is the distance field function of the $i^{th}$ triangular patch of the support domain’s boundary for node $a$. The gradient of the distance field function $\rho$ is given by:

[TABLE]

where $i$ is omitted for notation simplicity, $\nabla f=\widehat{n_{f}}=A_{f}\widehat{i}+B_{f}\widehat{j}+C_{f}\widehat{k}$ is the gradient of the triangle’s carrier plane $f$, which is identical to the plane’s normal vector, and $\nabla t$ is the gradient of the trim volume function:

[TABLE]

where $\nabla p_{j}$, $j=1,2,3$ is the gradient of the perpendicular plane to the $j^{th}$ edge of the triangle given by its normal vector $\widehat{n_{p_{j}}}$ and $\nabla G$ is given by:

[TABLE]

4 Evaluation Of The CME Approximants In The MTLED Method

We construct the matrix of the basis function derivatives using the CME approximation, instead of the MMLS, to derive the CME-MTLED method. We perform a series of numerical examples for three-dimensional large strain hyperelasticity at steady state to assess the efficiency of the CME approximants in the MTLED method. In all the examples, we use the hyper-elastic neo-Hookean material which is described by the hyper-elastic strain energy density function [42]:

[TABLE]

where $\lambda$ and $\mu$ are the Lamé parameters, $I_{1}$ is the first strain invariant and $J$ is the determinant of the deformation gradient. The evaluation of spatial integrals is applied on quadrature points, generated on a background unstructured tetrahedral mesh using a 4-point quadrature rule [43]. We use dynamic relaxation, described in detail in [44, 45], to ensure fast convergence for the explicit time integration scheme to the steady state solution. All the simulations are performed using a multi-threaded implementation of the MTLED algorithm on an Intel Core i7-8700K CPU using 12 threads and 16GB DDR4 RAM.

4.1 Cube Under Unconstrained Compression

We model the compression of an unconstrained cube with soft-tissue-like constitutive properties (hyper-elastic neo-Hookean material). The cube has side length $l=0.1m$, Young modulus $YM=3000Pa$, Poisson ratio $v=0.49$, and density $\rho=1000kg/m^{3}$. It is subjected to unconstrained compression (Figure 4) by displacing the top surface by $0.02m$ ($20\%$ of the initial height). In this simple case, the vertical component of displacement ($u_{z}$) is given by the analytical solution $u_{z}=-0.2z$, ($z$ refers to the reference configuration).

Simulations are performed for two successively fined nodal distributions, consisting of $152$ and $4594$ nodes, respectively. Computations are conducted using the CME-MTLED method. We evaluate the accuracy of the proposed scheme using the Normalized Root Mean Square Error ($NRMSE$):

[TABLE]

where $N$ is the number of nodes, $u_{z_{i}}^{h}$ is the $z$-axis displacement component of the numerical solution, and $u_{z_{i}}^{an}$ the $z$-axis displacement component of the analytical solution. Table 1 shows the computed $NRMSE$ for the two considered nodal distributions. Good agreement with the analytical solution is achieved even for the coarse grid ($152$ nodes, $1896$ quadrature points). Additionally, due to the weak Kronecker-delta property of the CME approximants, the essential boundary conditions defined on faces $p1$, $p2$, $p3$ and $p4$ are imposed exactly (absolute error is zero up to machine precision).

4.2 Cube Under Constrained Deformation

We further demonstrate the accuracy of the proposed scheme by considering the extension and compression of a cube with side length $l=0.1m$. One face ($z=0m$) of the cube is rigidly constrained, while the opposite face ($z=0.1m$) is displaced (Figure 5). A maximum displacement loading of $0.05m$ ($50\%$ of the initial height) is smoothly applied. A hyper-elastic neo-Hookean material model with $YM=3000Pa$, $v=0.49$ and $\rho=1000kg/m^{3}$ is chosen to capture the material response of soft tissue.

The cube is discretized using six successively denser nodal distributions ($lvl_{1}$: $152$, $lvl_{6}$: $50521$). The CME approximants are used to approximate the unknown field functions. Simulations are performed to test the convergence of the proposed scheme. The simulation characteristics such as the critical explicit integration time step ($\Delta t_{crit}$), the number of execution steps ($N_{exe}$), and the total execution time ($t_{exe}$) are evaluated for CME with $N_{R}=2$ and parameters $s=\{2,4\}$; $m=3$; $a=\{0.2,1.6,2.0\}$; $k=2$. The simulations are also performed using the MMLS approximants with the SEBCIEM correction [18]. The dilatation coefficient $a=1.6$ (parameter controlling support size) is used in MMLS. The simulation characteristics are given in Table 2.

The $\Delta t_{crit}$ for CME is found higher for all the parameter sets compared to MMLS-SEBCIEM. The $\Delta t_{crit}$ gain is reduced for increasing $s$ and $a$; with increasing $s$ the gradient becomes steeper, leading to larger eigenvalues in the spatial derivatives matrix ($\leftidx{{}^{t}_{0}}{\bm{B}}_{L0}^{t}$). The CME performs better, in terms of computational efficiency, since it can reach steady state with less time steps than the MMLS-SEBCIEM.

An analytical solution is not available for the given problem. Therefore, to evaluate the convergence and accuracy of the method we use the Signed Relative Error in strain energy density ($SRE_{W}$) as a convergence measure similar to [20]. We compute the $SRE_{W}$ for $lvl_{1}$ to $lvl_{5}$ discretization with respect to the $lvl_{6}$ (finest) discretization (Figure 6). The $SRE_{W}$ is obtained by:

[TABLE]

where $\bar{W}_{lvl_{i}}$ is the mean value of the strain energy density for the $lvl_{i}$ discretization, $i=1,\ldots,5$. The $SRE_{W}$ for CME is comparable with the $SRE_{W}$ for MMLS-SEBCIEM. While the MMLS-SEBCIEM lead to lower $SRE_{W}$ for higher nodal density ($lvl_{3}$-$lvl_{5}$), the CME results in lower $SRE_{W}$ for lower nodal density ($lvl_{1}$ and $lvl_{2}$).

For additional comparison, the $SRE_{W}$ for $lvl_{1}$ to $lvl_{5}$ discretization for a Finite Element Method (FEM) simulation is given. FEM simulations are performed with the FEBio v$2.5.2$ software [46] using isoparametric tetrahedral elements and the Newton-Raphson implicit time integration method. Linear tetrahedral elements are known for being prone to volumetric “locking”. For this reason, FEBio implements a nodally integrated tetrahedron with enhanced performance for finite deformation and near-incompressibility compared to the standard constant strain tetrahedron [47]. The $SRE_{W}$ for FEM simulations is found an order of magnitude higher from CME in $50\%$ extension.

Similar results are acquired for $50\%$ compression. The use of CME in MTLED leads to higher $\Delta t_{crit}$ and lower $t_{exe}$ compared to the MMLS-SEBCIEM (Table 2). The $SRE_{W}$ in $50\%$ compression is found lower for CME compared to MMLS-SEBCIEM and FEM for all the discretization levels (Figure 6). The set-up of $N_{R}=2$ and $s=2$; $m=3$; $a=2.0$; $k=2$ for CME is found to be a good trade-off between accuracy and efficiency in MTLED for both cases of extension and compression (Table 2). The use of CME instead of the MMLS-SEBCIEM in MTLED eliminates the necessity for EBC correction and the execution time is reduced. The execution time reduction is more evident for dense nodal discretization (up to $1295.9s$ - $22.7\%$). Moreover, the evaluation of the $SRE_{W}$ demonstrates the improved accuracy of the MTLED method over the FEM for large strain problems for both $50\%$ extension and compression conditions. Especially for severe distortion, such as in the $50\%$ compression case, the FEM simulation leads to poor results for coarse nodal discretization compared to MTLED using either CME or MMLS-SEBCIEM approximants (Figure 7).

4.3 Cylinder Under Locally Applied Indentation

We highlight the applicability of the proposed CME-MTLED method against extreme indentation and demonstrate its applicability beyond what is possible with FEM. In detail we consider indentation of a cylindrical domain with height $h=17mm$ and diameter $d=30mm$, modelled by the hyper-elastic neo-Hookean constitutive law as in [16]. The sub-region of the cylinder’s top surface (illustrated in Figure 8) undergoes a vertical displacement of $10mm$ ($60\%$ of the initial height) while the nodes belonging to the bottom surface are fixed to zero displacements at all directions. The assigned material properties are similar to the test cases in Subsections 4.1 and 4.2 ($YM=3000Pa$; $v=0.49$; $\rho=1000kg/m^{3}$). The cylindrical domain consists of $8223$ nodes and $177372$ quadrature points (4-point quadrature rule). The CME approximants are constructed using the set-up of $N_{R}=2$ and $s=2$; $m=3$; $a=2.0$; $k=2$. The deformation of the indented cylinder at steady state is given in Figure 8. The described displacements at the bottom and top surfaces are exactly satisfied.

The effect of the quadrature on the accuracy of the proposed method is evaluated. The mean strain energy density is measured for several simulations using 1, 4, 8, 16, and 32 quadrature points per integration cell. Quadrature with 8, 16, and 32 quadrature points is performed by using the 4-point quadrature rule after subdividing the cells of the original background mesh using 2, 4, and 8 subdivisions respectively. It is shown that for increasing the number of the quadrature points, the strain energy density is reduced monotonically (Figure 9). The difference of the mean strain energy density for quadrature using 4 and 32 points per integration cell is found $3.9\%$. Therefore, performing integration with the 4-point quadrature rule appears to be a good compromise between accuracy and efficiency. Additional improvement of the accuracy with minimum computational overhead could be achieved by using the adaptive integration method described in [48].

4.4 Cardiac Model Geometry Under Locally Applied Indentation

Changes in myocardium stiffness are related to different cardiac conditions or diseases (e.g., diastolic dysfunction [49], tachycardia-induced cardiomyopathy [50]). Evaluation of myocardial material properties under such circumstances is relevant for disease assessment and surgical planning. Material properties evaluation is based on the estimation of the force/displacement relation for a predefined displacement [51, 52]. In this context, we simulate the indentation of a 3D cardiac bi-ventricular geometry. Our objective is to demonstrate the applicability of the CME-MTLED method on domains with arbitrary shape and its ability to satisfy the prescribed displacement conditions on such domains. We use a mesh composed of $38764$ nodes and $186561$ tetrahedral cells which was generated from MRI cardiac segmentation images using TetGen [53]. The segmentation images were made available from [54]. Tissue indentation in the $x$-axis direction is simulated imposing a $u_{x}=0.01m$ displacement on a patch region at the right cardiac ventricle while the bottom and top faces of the ventricles are constrained in all directions (Figure 10). We use the neo-Hookean constitutive model with parameters $YM=3000Pa$; $v=0.49$; $\rho=1000kg/m^{3}$ to describe the mechanical response of the cardiac tissue. We assume that the cardiac bi-ventricular geometry is not beating. It should be noted that the selected essential boundary conditions and constitutive model are not based on available experimental data and do not represent neither in-vivo boundary conditions nor the realistic response of the heart. They are rather chosen aiming to demonstrate the capabilities of the CME-MTLED method. A realistic biomechanical simulation is out-of-scope of the current study.

The indentation simulation is performed using both the CME ($s=2$, $a=2.0$) and the MMLS-SEBCIEM. The characteristics $\Delta t_{crit}$, $N_{exe}$, $t_{exe}$ and the EBC values at steady state for the two approximation schemes are given in Table 3. It is shown that using CME, the EBC can be imposed exactly without the need of any special treatment such as SEBCIEM or EBCIEM. In addition, solution stability can be achieved using longer explicit integration time steps leading to a $27\%$ reduction in the computational time.

5 Concluding remarks

In the present study, we presented the extension of Cell-based Maximum Entropy (CME) approximants in the three-dimensional Euclidean space ($E^{3}$) using the R-functions technology based on the work of Millán et al [28], and Biswas and Shapiro [41]. The CME approximants were used specifically in the MTLED method, however they could be easily introduced in other solvers (explicit, implicit, semi-implicit). We evaluated the accuracy and efficiency of the CME in the Meshless Total Lagrangian Explicit Dynamics (MTLED) method through comparison with Modified Moving Least Squares (MMLS) and Finite Element Method (FEM). We presented two numerical examples to verify the capability of CME MTLED to generate accurate solutions to nonlinear equations of solid mechanics governing the behavior of soft, deformable tissues. We also verified the accuracy of the proposed scheme against extreme indentation, with the indentation depth reaching $60\%$ of the initial height of the sample. The evaluation was performed by analyzing the signed relative error in strain energy density for several combinations of parameters used in CME. Finally, we verified the accuracy of essential boundary conditions (EBC) imposition on domains of arbitrary shape in an indentation simulation of a cardiac bi-ventricular model derived from cardiac MRI segmentation images.

In all the numerical examples, the use of the CME approximants allowed longer explicit time integration step without compromising the numerical stability. A gain in computational time up to $27\%$ was achieved for CME, while demonstrating similar convergence to MMLS. The weak Kronecker-delta property of CME allowed to directly impose EBC, avoiding special EBC imposition treatments, which have been previously proposed to deal with the lack of the Kronecker delta property in MMLS. The smooth derivatives of CME led to derivatives matrix with smaller eigenvalues compared to MMLS and hence longer explicit time integration steps. CME was proven a valuable alternative to the group of MLS approximants in the context of MTLED and other similar EFG frameworks. A naive implementation of the CME approximants in $E^{3}$ is available to download at https://www.mountris.org/software/mlab/cme.

Acknowledgments

This work was supported by the European Research Council through grant ERC-2014-StG 638284, by MINECO (Spain) through project DPI2016-75458-R and by European Social Fund (EU) and Aragón Government through BSICoS group (T39_17R). Computations were performed by the ICTS NANBIOSIS (HPC Unit at University of Zaragoza). The authors wish to thank CONICET and grant PICTO-2016-0054 from UNCuyo-ANPCyT for funding the third author and the Department of Health, Western Australia, for its financial support to the fourth author through a Merit Award. The fifth and last authors acknowledge the funding from the Australian Government through the Australian Research Council ARC (Discovery Project Grants DP160100714, DP1092893, and DP120100402).