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

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.



There are no comments yet.



Weak imposition of Signorini boundary conditions on the boundary element method

We derive and analyse a boundary element formulation for boundary condit...

Damping perturbation based time integration asymptotic method for structural dynamics

The light damping hypothesis is usually assumed in structural dynamics s...

Simple and robust element-free Galerkin method with interpolating shape functions for finite deformation elasticity

In this paper, we present a meshless method belonging to the family of e...

Mollified finite element approximants of arbitrary order and smoothness

The approximation properties of the finite element method can often be s...

Maximum Entropy Flow Networks

Maximum entropy modeling is a flexible and popular framework for formula...

The Projection Extension Method: A Spectrally Accurate Technique for Complex Domains

An essential ingredient of a spectral method is the choice of suitable b...
This week in AI

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

1 Introduction

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 . A variational approach to adapt 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 -rings () 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 .

The purpose of this work is to extend the CME approximants in 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 , based on the minimum relative entropy framework. In Section 3, we extend the Cell-based Maximum Entropy approximation in . 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 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 to the boundary of the support ring, , 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].

Figure 1: 1st and 2nd ring support for a field node (yellow) in (a) two and (b) three dimensional domain.

2.1 Minimum Relative Entropy Framework

For a scalar-valued function and a set of unstructured points , the MaxEnt approximation is given by


where are nodal coefficients, and are nonnegative basis functions that fulfill the zeroth and first order reproducing conditions:


The basis functions are defined from an optimization problem that is established at each evaluation point 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:


where 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:

subject to

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




is the partition function and

is the Lagrange multiplier vector. The Lagrange multiplier vector is obtained solving the unconstrained convex optimization problem:


where is the converged solution of at . The resulting basis functions from the solution of 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 evaluated on a point is given by:


where the superscript on any function indicates evaluation at , is the prior weight function for node , , and is given by:



is the identity matrix,

, and . A detailed description of the derivation of is given in Refs. [28, 34].

Figure 2: Cell-based Maximum Entropy (CME) basis function for a point located (a) at the center, with coordinates , and (b) at the right edge, with coordinates , of a two-dimensional rectangular domain.

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 -ring support domain (Figure 1). Any real-valued function , where , is considered a R-function if its sign is determined only by the sign of its arguments 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 from the boundary of its polygonal support domain is expressed as the composition of elementary R-functions corresponding to the piecewise linear segments of the polygonal curve . For each line segment belonging to with endpoints and , the signed distance function from a point to the line passing from and is defined by:


where is the length of the line segment. In addition, the line segment defines a disk of radius . This disk can be expressed by the trimming function:


where is normalized to first order and defines the disk with center . From Equations (10), (11) the signed distance function from a point to the line segment with endpoints and is given by the first order normalized function :


where is differentiable for any point that does not lay on the line segment. For line segments s.t. the normalized approximation, up to order , of the distance function for the node is given by the R-equivalence formula [41]:


The R-equivalence formula has the desirable properties of preserving the normalization up to the order and being associative, while the R-conjunction is normalized up to the order and it is not associative. Finally, the prior weight function for node is obtained by:


where are the indices of the -ring nodal neighbors of the point , is a smoothness modulation factor and is normalized s.t. . 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 , the bulging effect is reduced.

3 Extension Of Cell-based Maximum Entropy Approximants In E3

MaxEnt approximants (Equation (5)) are naturally generalized to , with . However, CME approximants have been limited to 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 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 of a finite domain in tetrahedra . The ring support domain for any vertex of is defined by the union of the piecewise triangular faces belonging to the ring of attached tetrahedra to the vertex (Figure 1). R-functions are employed to construct the approximation distance field to the boundary of the 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 and a trim volume defined by the planes , that pass through the edges of the triangle and are orthogonal to the carrier plane (Figure 3). The trim volume can be constructed joining the planes according to the -conjunction formula:


where , are constants controlling the shape of the R-function’s zero set. For points in that do not lay on the carrier plane the Equation (15) is , while for points on the carrier plane () it reduces to the standard -conjunction formula. Therefore, is given by:


where . The function defines a smooth trimming volume tangent to the three edges of the triangle. The coefficient controls the smoothness of the trimming volume, where for large the trimming volume flattens out and for it becomes identical to the unbounded prism constructed by the combination of the planes (Figure 3). While the unbounded prism trimming volume () does not maintain differentiability on the triangle edges, the smooth volume () is differentiable at any point except the triangle vertices where the bulging effect arises. Substituting and in Equation (12), the normalized distance function from any point to the triangle is approximated and the nodal prior weights are obtained by Equations (13), (14), similarly to the case.

Figure 3: (a) Prism trim volume for a triangle constructed by planes , , (). (b - d) Smooth trimming volume tangent to the edges of the triangle for . Parameter .

3.2 Derivation Of The CME Approximants Gradient In E3

To compute the gradient of the CME approximants (Equation (8)), we should first compute the gradient of the prior function . Following the abbreviation in [28], to derive , for a node at a point , we first derive the gradient of the normalized approximation to the distance function:


where is the distance field function of the triangular patch of the support domain’s boundary for node . The gradient of the distance field function is given by:


where is omitted for notation simplicity, is the gradient of the triangle’s carrier plane , which is identical to the plane’s normal vector, and is the gradient of the trim volume function:


where , is the gradient of the perpendicular plane to the edge of the triangle given by its normal vector and is given by:


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]:


where and are the Lamé parameters, is the first strain invariant and 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 , Young modulus , Poisson ratio , and density . It is subjected to unconstrained compression (Figure 4) by displacing the top surface by ( of the initial height). In this simple case, the vertical component of displacement () is given by the analytical solution , ( refers to the reference configuration).

Figure 4: Boundary conditions for unconstrained cube under compression. Face p1: , and unconstrained, face p2: , and unconstrained, face p3: , and unconstrained, and face p4: , and unconstrained.

Simulations are performed for two successively fined nodal distributions, consisting of and 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 ():


where is the number of nodes, is the -axis displacement component of the numerical solution, and the -axis displacement component of the analytical solution. Table 1 shows the computed for the two considered nodal distributions. Good agreement with the analytical solution is achieved even for the coarse grid ( nodes, quadrature points). Additionally, due to the weak Kronecker-delta property of the CME approximants, the essential boundary conditions defined on faces , , and are imposed exactly (absolute error is zero up to machine precision).

Nodes Quadrature points
152 1896 1.75E-3
4594 92876 3.85E-4
Table 1: Normalized Root Mean Square Error () for a cube under unconstrained compression.

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 . One face () of the cube is rigidly constrained, while the opposite face () is displaced (Figure 5). A maximum displacement loading of ( of the initial height) is smoothly applied. A hyper-elastic neo-Hookean material model with , and is chosen to capture the material response of soft tissue.

The cube is discretized using six successively denser nodal distributions (: , : ). 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 (), the number of execution steps (), and the total execution time () are evaluated for CME with and parameters ; ; ; . The simulations are also performed using the MMLS approximants with the SEBCIEM correction [18]. The dilatation coefficient (parameter controlling support size) is used in MMLS. The simulation characteristics are given in Table 2.

Figure 5: Steady state deformation under extension for (a) , and (b) discretizations. Nodes in the bottom surface are fixed and nodes at the top surface undergo a displacement in the -axis, CME, MMLS-SEBCIEM, FEM.

The for CME is found higher for all the parameter sets compared to MMLS-SEBCIEM. The gain is reduced for increasing and ; with increasing

the gradient becomes steeper, leading to larger eigenvalues in the spatial derivatives matrix (

). The CME performs better, in terms of computational efficiency, since it can reach steady state with less time steps than the MMLS-SEBCIEM.

Approx. (ms) (s)
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
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
  • CME smoothness modulation factor

  • CME R-function zero set shape factor or MMLS dilatation coefficient

Table 2: Simulation characteristics (critical explicit integration time step - , number of execution steps - , execution time - ) for cube under uniaxial extension compression for - discretization.

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 () as a convergence measure similar to [20]. We compute the for to discretization with respect to the (finest) discretization (Figure 6). The is obtained by:


where is the mean value of the strain energy density for the discretization, . The for CME is comparable with the for MMLS-SEBCIEM. While the MMLS-SEBCIEM lead to lower for higher nodal density (-), the CME results in lower for lower nodal density ( and ).

Figure 6: Signed relative error in strain energy density for a cube at: (a) uniaxial extension, (b) uniaxial compression.

For additional comparison, the for to discretization for a Finite Element Method (FEM) simulation is given. FEM simulations are performed with the FEBio v 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 for FEM simulations is found an order of magnitude higher from CME in extension.

Similar results are acquired for compression. The use of CME in MTLED leads to higher and lower compared to the MMLS-SEBCIEM (Table 2). The in compression is found lower for CME compared to MMLS-SEBCIEM and FEM for all the discretization levels (Figure 6). The set-up of and ; ; ; 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 - ). Moreover, the evaluation of the demonstrates the improved accuracy of the MTLED method over the FEM for large strain problems for both extension and compression conditions. Especially for severe distortion, such as in the 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).

Figure 7: Steady state deformation under compression for (a) , and (b) discretizations. Nodes in the bottom surface are fixed and nodes at the top surface undergo a displacement in the -axis, CME, MMLS-SEBCIEM, FEM.

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 and diameter , 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 ( 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 (; ; ). The cylindrical domain consists of nodes and quadrature points (4-point quadrature rule). The CME approximants are constructed using the set-up of and ; ; ; . 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.

Figure 8: (a) Fixed (blue) and displaced (red) boundary nodes for locally applied indentation on a cylinder. (b) Steady state deformation of the cylinder under locally applied indentation. The colormap represents the resulting displacements at Z-axis direction for 177372 quadrature points.

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 . 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].

Figure 9: Mean values of strain energy for increasing quadrature of a locally applied indentation to a cylinder. The strain energy reduces monotonically for increasing quadrature points. Mean strain energy values are normalized to the maximum value (1 quadrature point) to aid the results’ interpretation.

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 nodes and 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 -axis direction is simulated imposing a 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 ; ; 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.

Figure 10: (a) Nodes (blue) at the top and bottom faces of the cardiac bi-ventricular model are fixed in all directions. The nodes on a small patch located at the right ventricle (red) are displaced by . (b) Steady state deformation during right ventricle indentation computed with the CME-MTLED using ().

The indentation simulation is performed using both the CME (, ) and the MMLS-SEBCIEM. The characteristics , , 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 reduction in the computational time.

Approx. (ms) (h) Displaced EBC (m) Fixed EBC (m)
CME 0.636 36073 1.06 -0.01 1.0E-17 0.00 0.00
MMLS-SEBCIEM 0.455 49692 1.53 -0.01 3.3E-15 0.00 1.7E-15
  • Mean EBC values at steady state for fixed and displaced nodes with mean variation

Table 3: Simulation characteristics (critical explicit integration time step - , number of execution steps - , execution time - , Essential boundary conditions (EBC) values) for cardiac bi-ventricular model in indentation.

5 Concluding remarks

In the present study, we presented the extension of Cell-based Maximum Entropy (CME) approximants in the three-dimensional Euclidean space () 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 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 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.


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).


  • [1] Sahil Garg and Mohit Pant. Meshfree Methods: A Comprehensive Review of Applications. International Journal of Computational Methods, 15(04):1830001, jun 2018.
  • [2] Adam Wittek, Nicole M. Grosland, Grand Roman Joldes, Vincent Magnotta, and Karol Miller. From Finite Element Meshes to Clouds of Points: A Review of Methods for Generation of Computational Biomechanics Models for Patient-Specific Applications. Annals of Biomedical Engineering, 44(1):3–15, jan 2016.
  • [3] L. B. Lucy and L.B. A numerical approach to the testing of the fission hypothesis. The Astronomical Journal, 82:1013, dec 1977.
  • [4] R.A Gingold and J.J Monaghan. Kernel estimates as a basis for general particle methods in hydrodynamics. Journal of Computational Physics, 46(3):429–453, jun 1982.
  • [5] P. H. Wen and M. H. Aliabadi. An improved meshless collocation method for elastostatic and elastodynamic problems. Communications in Numerical Methods in Engineering, 24(8):635–651, feb 2007.
  • [6] X. Zhang, K. Z. Song, M. W. Lu, and X. Liu.

    Meshless methods based on collocation with radial basis functions.

    Computational Mechanics, 26(4):333–343, oct 2000.
  • [7] G.C. Bourantas, K.A. Mountris, V.C. Loukopoulos, L. Lavier, G.R. Joldes, A. Wittek, and K. Miller. Strong-form approach to elasticity: Hybrid finite difference-meshless collocation method (FDMCM). Applied Mathematical Modelling, 57:316–338, may 2018.
  • [8] S. N. Atluri and T. Zhu. A new Meshless Local Petrov-Galerkin (MLPG) approach in computational mechanics. Computational Mechanics, 22(2):117–127, aug 1998.
  • [9] T. Belytschko, Y. Y. Lu, and L. Gu. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering, 37(2):229–256, jan 1994.
  • [10] GR Liu and YT Gu. An introduction to meshfree methods and their programming. Springer Science & Business Media, 2005.
  • [11] Nicolas Libre, Arezoo Emdadi, Edward Kansa, Mohammad Rahimian, and Mohammad Shekarchi. A Stabilized RBF Collocation Scheme for Neumann Type Boundary Value Problems. CMES - Computer Modeling in Engineering and Sciences, 24(1):61–80, jan 2008.
  • [12] De’an Hu and Zhanhua Sun. The Meshless Local Petrov-Galerkin Method for Large Deformation Analysis of Hyperelastic Materials. ISRN Mechanical Engineering, 2011:1–11, sep 2011.
  • [13] Yumin Cheng, Funong Bai, Chao Liu, and Miaojuan Peng. Analyzing nonlinear large deformation with an improved element-free Galerkin method via the interpolating moving least-squares method. International Journal of Computational Materials Science and Engineering, 05(04):1650023, dec 2016.
  • [14] Ashley Horton, Adam Wittek, Grand Roman Joldes, and Karol Miller. A meshless Total Lagrangian explicit dynamics algorithm for surgical simulation. International Journal for Numerical Methods in Biomedical Engineering, 26(8):977–998, mar 2010.
  • [15] Mao Li, Karol Miller, Grand Roman Joldes, Ron Kikinis, and Adam Wittek. Biomechanical model for computing deformations for whole-body image registration: A meshless approach. International Journal for Numerical Methods in Biomedical Engineering, 32(12):e02771, dec 2016.
  • [16] Grand Joldes, George Bourantas, Benjamin Zwick, Habib Chowdhury, Adam Wittek, Sudip Agrawal, Konstantinos Mountris, Damon Hyde, Simon K. Warfield, and Karol Miller. Suite of meshless algorithms for accurate computation of soft tissue deformation for surgical simulation. Medical Image Analysis, 56:152–171, aug 2019.
  • [17] Grand Roman Joldes, Habibullah Amin Chowdhury, Adam Wittek, Barry Doyle, and Karol Miller. Modified moving least squares with polynomial bases for scattered data approximation. Applied Mathematics and Computation, 266:893–902, sep 2015.
  • [18] Grand Roman Joldes, Habib Chowdhury, Adam Wittek, and Karol Miller. A new method for essential boundary conditions imposition in explicit meshless methods. Engineering Analysis with Boundary Elements, 80:94–104, jul 2017.
  • [19] N. Sukumar. Construction of polygonal interpolants: a maximum entropy approach. International Journal for Numerical Methods in Engineering, 61(12):2159–2181, nov 2004.
  • [20] M. Arroyo and M. Ortiz. Local maximum-entropy approximation schemes: a seamless bridge between finite elements and meshfree methods. International Journal for Numerical Methods in Engineering, 65(13):2167–2202, mar 2006.
  • [21] E. T. Jaynes. Information Theory and Statistical Mechanics. Physical Review, 106(4):620–630, may 1957.
  • [22] A. Rosolen, C. Peco, and M. Arroyo. An adaptive meshfree method for phase-field models of biomembranes. Part I: Approximation with maximum-entropy basis functions. Journal of Computational Physics, 249:303–319, sep 2013.
  • [23] A. Rosolen and M. Arroyo. Blending isogeometric analysis and local maximum entropy meshfree approximants. Computer Methods in Applied Mechanics and Engineering, 264:95–107, sep 2013.
  • [24] Z. Ullah, C.E. Augarde, and W. M. Coombs. Local Maximum Entropy Shape Functions Based FE-EFGM Coupling. Technical report, University of Glasgow, College of Science and Engineering, 2013.
  • [25] David González, Elías Cueto, and Manuel Doblaré. A higher order method based on local maximum entropy approximation. International Journal for Numerical Methods in Engineering, 83(6):n/a–n/a, aug 2010.
  • [26] Daniel Millán, Adrian Rosolen, and Marino Arroyo. Nonlinear manifold learning for meshfree finite deformation thin-shell analysis. International Journal for Numerical Methods in Engineering, 93(7):685–713, feb 2013.
  • [27] Adrian Rosolen, Daniel Millán, and Marino Arroyo. On the optimum support size in meshfree methods: A variational adaptivity approach with maximum-entropy approximants. International Journal for Numerical Methods in Engineering, 82(7):n/a–n/a, may 2009.
  • [28] Daniel Millán, N. Sukumar, and Marino Arroyo. Cell-based maximum-entropy approximants. Computer Methods in Applied Mechanics and Engineering, 284:712–731, feb 2015.
  • [29] Vadim Shapiro. Theory of R-functions and Applications: A Primer. Cornell University, 1991.
  • [30] J. Shore and R. Johnson. Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy. IEEE Transactions on Information Theory, 26(1):26–37, jan 1980.
  • [31] Marino Arroyo and Michael Ortiz. Local Maximum-Entropy Approximation Schemes. In

    Meshfree Methods for Partial Differential Equations III

    , pages 1–16. Springer Berlin Heidelberg, Berlin, Heidelberg, 2007.
  • [32] N. Sukumar and R. J.-B. Wets. Deriving the Continuity of Maximum-Entropy Basis Functions via Variational Analysis. SIAM Journal on Optimization, 18(3):914–925, jan 2007.
  • [33] B. Li, F. Habbal, and M. Ortiz. Optimal transportation meshfree approximation schemes for fluid and plastic flows. International Journal for Numerical Methods in Engineering, 83(12):1541–1579, sep 2010.
  • [34] L. L. Yaw, N. Sukumar, and S. K. Kunnath. Meshfree co-rotational formulation for two-dimensional continua. International Journal for Numerical Methods in Engineering, 79(8):979–1003, aug 2009.
  • [35] A. Ortiz, M.A. Puso, and N. Sukumar. Maximum-entropy meshfree method for compressible and near-incompressible elasticity. Computer Methods in Applied Mechanics and Engineering, 199(25-28):1859–1871, may 2010.
  • [36] J.S. Hale and P.M. Baiz. A locking-free meshfree method for the simulation of shear-deformable plates based on a mixed variational formulation. Computer Methods in Applied Mechanics and Engineering, 241-244:311–322, oct 2012.
  • [37] K. Hormann and N. Sukumar. Maximum Entropy Coordinates for Arbitrary Polytopes. Computer Graphics Forum, 27(5):1513–1520, jul 2008.
  • [38] N. Sukumar. Quadratic maximum-entropy serendipity shape functions for arbitrary planar polygons. Computer Methods in Applied Mechanics and Engineering, 263:27–41, aug 2013.
  • [39] Keijo Nissen, Christian J. Cyron, Volker Gravemeier, and Wolfgang A. Wall. Information-flux method: a meshfree maximum-entropy Petrov–Galerkin method including stabilised finite element methods. Computer Methods in Applied Mechanics and Engineering, 241-244:225–237, oct 2012.
  • [40] C. T. Wu, D. L. Young, and H. K. Hong. Adaptive meshless local maximum-entropy finite element method for convection–diffusion problems. Computational Mechanics, 53(1):189–200, jan 2014.
  • [41] Arpan Biswas and Vadim Shapiro. Approximate distance fields with non-vanishing gradients. Graphical Models, 66(3):133–159, may 2004.
  • [42] Javier Bonet and Richard D. Wood. Nonlinear continuum mechanics for finite element analysis. Cambridge University Press, Cambridge, 2008.
  • [43] Ronald Cools and Philip Rabinowitz. Monomial cubature rules since “Stroud”: a compilation. Journal of Computational and Applied Mathematics, 48(3):309–326, nov 1993.
  • [44] Grand Roman Joldes, Adam Wittek, and Karol Miller. Computation of intra-operative brain shift using dynamic relaxation. Computer Methods in Applied Mechanics and Engineering, 198(41-44):3313–3320, sep 2009.
  • [45] Grand Roman Joldes, Adam Wittek, and Karol Miller. An adaptive dynamic relaxation method for solving nonlinear finite element problems. Application to brain shift estimation. International Journal for Numerical Methods in Biomedical Engineering, 27(2):173–185, feb 2011.
  • [46] Steve A. Maas, Benjamin J. Ellis, Gerard A. Ateshian, and Jeffrey A. Weiss. FEBio: Finite Elements for Biomechanics. Journal of Biomechanical Engineering, 134(1):011005, jan 2012.
  • [47] M. A. Puso and J. Solberg. A stabilized nodally integrated tetrahedral. International Journal for Numerical Methods in Engineering, 67(6):841–867, aug 2006.
  • [48] Grand Roman Joldes, Adam Wittek, and Karol Miller. Adaptive numerical integration in Element-Free Galerkin methods for elliptic boundary value problems. Engineering Analysis with Boundary Elements, 51:52–63, feb 2015.
  • [49] G Jackson, C R Gibbs, M K Davies, and G Y Lip. ABC of heart failure. Pathophysiology. BMJ (Clinical research ed.), 320(7228):167–16770, jan 2000.
  • [50] Atul Khasnis, Krit Jongnarangsin, George Abela, Srikar Veerareddy, Vivek Reddy, and Ranjan Thakur. Tachycardia-Induced Cardiomyopathy: A Review of Literature. Pacing and Clinical Electrophysiology, 28(7):710–721, jul 2005.
  • [51] Toshiaki Shishido, Masaru Sugimachi, Osamu Kawaguchi, Hiroshi Miyano, Toru Kawada, Wataru Matsuura, Yasuhiro Ikeda, Takayuki Sato, Joe Alexander, and Kenji Sunagawa. A new method to measure regional myocardial time-varying elastance using minute vibration. American Journal of Physiology-Heart and Circulatory Physiology, 274(4):H1404–H1415, apr 1998.
  • [52] D M Blair and H R Halperin. Hand-held, dynamic indentation system for measuring myocardial transverse stiffness. Biomedical instrumentation & technology, 30(6):517–25, 1996.
  • [53] Hang Si and Hang. TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator. ACM Transactions on Mathematical Software, 41(2):1–36, feb 2015.
  • [54] Wenjia Bai, Wenzhe Shi, Antonio de Marvao, Timothy J.W. Dawes, Declan P. O’Regan, Stuart A. Cook, and Daniel Rueckert. A bi-ventricular cardiac atlas built from 1000+ high resolution MR images of healthy subjects and an analysis of shape and motion. Medical Image Analysis, 26(1):133–145, dec 2015.