# Faster and objective Level Set-DEM mechanical simulations of discrete systems with convex particles from contact history and particle surface considerations

Jérôme Duriez, Cédric Galusinski

► To cite this version:

Jérôme Duriez, Cédric Galusinski. Faster and objective Level Set-DEM mechanical simulations of discrete systems with convex particles from contact history and particle surface considerations. Computer Physics Communications, 2025, pp.109803. 10.1016/j.cpc.2025.109803. hal-05213383

HAL Id: hal-05213383

https://hal.inrae.fr/hal-05213383v1

Submitted on 18 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

Contents lists available at ScienceDirect

Computer Physics Communications

journal homepage: www.elsevier.com/locate/cpc

<figure><img src="image_3.png" /></figure>

Computational Physics

Faster and objective level set-DEM mechanical simulations of discrete
systems with convex particles from contact history and particle surface
considerations

Jérôme Duriez<sup>a, ID, *</sup>, Cédric Galusinski<sup>b</sup>

ARTICLE INFO

The review of this paper was arranged by
Prof. Andrew Hazel

Keywords:

Level set discrete element method (LS-DEM)

Superellipsoids

Superquadrics

Convex shape

A B S T R A C T

Aiming for versatile simulations of the mechanics of discrete systems with arbitrary convex particles, an extended Level Set (LS) description of particle shape is proposed in the framework of the Discrete Element Method (DEM). The LS shape description as a discrete field of the signed distance function is first obtained with a faster initial generation and then proposed to directly output particle surface in a validated workflow. It mostly includes an innovative optimization for the surface nodes discretization which combines with the LS distance field for DEM contact treatment. In their optimized definition, surface nodes locate in a nearly uniform fashion over a particle and are handled in a sparse manner during contact treatment, thanks to an original consideration of contact history. As such, computation speed gains are reported with a factor of more than 3 during simulations of quasi-static mechanical loading. The proposed nodes definition is also shown to be instrumental to insure objectivity of the LS contact model and DEM simulations. The present LS approach is finally applied to a preliminary study of the residual shear strength of various packings made of superquadrics all sharing the same shape. After a justification of the chosen particle number to form a Representative Elementary Volume, the shear strength is shown to lack a clear relationship with Wadell's true sphericity as a shape index, while being possibly 73% higher than the one obtained with spherical particles.

1. Introduction

A number of physical processes spanning from the geotechnical
to the chemical engineering fields involve particulate materials, e.g.,
rock stones and pills, whose mechanics is dictated by a combination
of particle-scale, contact-scale and bulk-scale phenomena or properties
such as particle shape, contact friction or fabric. In line with the great
role of the latter bulk-scale geometric effects and thanks to its ability
to reproduce those in a straightforward manner, the Discrete Element
Method (DEM) [1] has shown a remarkable success in reproducing most
qualitative trends of granular mechanics even by just considering fric-
tionless spherical particles, e.g., the existence of a pressure-dependent
macroscopic shear strength [2,3]. However, truly predictive DEM simu-
lations still require a careful inclusion of the particle-scale and contact-
scale properties, e.g., [4–9] as it obviously appears on any mountainous
slope covered by centrimetric to metric rocks which would be impossi-

ble to observe with frictionless spheres, because of negligible bulk-scale
effects in this case of a thin cover.

While including appropriate contact-scale properties is mostly a mat-
ter of calibration against pertinent experiments, the due consideration
of the particle scale, i.e. a realistic shape with respect to imaging tech-
niques observations, is a technical challenge for the DEM method itself.
As a matter of fact, the theoretical and computational complexity of de-
tecting contact within the DEM cycle drastically increases as soon as
one leaves spherical shapes aside, which would be fully characterized
with just a scalar (radius or diameter) that is straightforward to compare
with interparticle distance. Such a convenient approach can be some-
what extended to more complex shapes after forming rigid multi-sphere
aggregates, i.e. clumps, following a number of available algorithms, see
e.g., [10–12]. In line with the existence of several such algorithms, the
creation of such a clump from a reference shape is nevertheless not a
straightforward technique in the general case [12]. It furthermore nec-

<footnote>a INRAE, Aix Marseille Univ, RECOVER, Aix-en-Provence, France</footnote>

<footnote>b IMATH, Université de Toulon, CS 60584, 83041 Toulon Cedex 9, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: jerome.duriez@inrae.fr (J. Duriez).</footnote>

<footnote>URL: https://eng-recover.paca.hub.inrae.fr/about-us/laboratory-s-members/social-pages/jerome-duriez (J. Duriez).</footnote>

<footnote>https://doi.org/10.1016/j.cpc.2025.109803</footnote>

<footnote>Received 4 March 2025; Received in revised form 29 July 2025; Accepted 31 July 2025</footnote>

<footnote>0010-4655/© 2025 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).</footnote>

<footnote>67</footnote>

<footnote>68</footnote>

<footnote>69</footnote>

<footnote>70</footnote>

<footnote>71</footnote>

<footnote>72</footnote>

<footnote>73</footnote>

<footnote>74</footnote>

<footnote>75</footnote>

<footnote>76</footnote>

<footnote>77</footnote>

<footnote>78</footnote>

<footnote>79</footnote>

<footnote>80</footnote>

<footnote>81</footnote>

<footnote>82</footnote>

<footnote>83</footnote>

<footnote>84</footnote>

<footnote>85</footnote>

<footnote>86</footnote>

<footnote>87</footnote>

<footnote>88</footnote>

<footnote>89</footnote>

<footnote>90</footnote>

<footnote>91</footnote>

<footnote>92</footnote>

<footnote>93</footnote>

<footnote>94</footnote>

<footnote>95</footnote>

<footnote>96</footnote>

<footnote>97</footnote>

<footnote>98</footnote>

<footnote>99</footnote>

<footnote>100</footnote>

<footnote>101</footnote>

<footnote>102</footnote>

<footnote>103</footnote>

<footnote>104</footnote>

<footnote>105</footnote>

<footnote>106</footnote>

<footnote>107</footnote>

<footnote>108</footnote>

<footnote>109</footnote>

<footnote>110</footnote>

<footnote>111</footnote>

<footnote>112</footnote>

<footnote>113</footnote>

<footnote>114</footnote>

<footnote>115</footnote>

<footnote>116</footnote>

<footnote>117</footnote>

<footnote>118</footnote>

<footnote>119</footnote>

<footnote>120</footnote>

<footnote>121</footnote>

<footnote>122</footnote>

<footnote>123</footnote>

<footnote>124</footnote>

<footnote>125</footnote>

<footnote>126</footnote>

<footnote>127</footnote>

<footnote>128</footnote>

<footnote>129</footnote>

<footnote>130</footnote>

<footnote>131</footnote>

<footnote>132</footnote>
essarily includes some approximation, in particular in terms of local roundness. As an alternative, a number of DEM implementations are available for polyhedral shapes, e.g., [13–16], while relying on different algorithms such as the Gilbert-Johnson-Keerthi algorithm (GJK) [17]. The consideration of such polyhedra has for instance been proven to be efficient for quantitative simulations of Toyoura sand in [7], without the need for enriching the contact parameters beyond the classical two stiffnesses and one Coulomb friction coefficient. As proposed in [18], the same GJK algorithm may also apply to superquadric shapes, i.e. superellipsoids, that are more generic than spheres but offer a somewhat similar analytical definition of their surfaces. Applying functional minimization techniques to superquadric surface equations led to alternate DEM implementations for such shapes in [19,20]. Among the remaining DEM approaches for complex shapes that have been more extensively reviewed in [21], a Level Set (LS) approach [22,4,23], also denoted as LS-DEM, appears of particular interest due to its versatility and possibly lighter computational complexity with respect to some polyhedra algorithms [24].

In a LS approach, particle surface is implicitly defined as the zero-level set of the signed (shortest) distance function (SDF) to the surface, with some similarity to [25,26] who considered other functions than the shortest distance. Contributing to the versatility of the LS method, the SDF is defined in a discretized fashion on a particle-attached grid, as one can numerically obtain for virtually any kind of surface thanks to, e.g., a Fast Marching Method (FMM) [27]. Contact detection for a pair of bodies 1-2 was then proposed in [22,4] to combine the distance data of 1 with an additional surface discretization in terms of nodes for 2, where distance queries would be carried out by looping over all nodes. The set of surface nodes obviously controls the precision and time costs of contact detection [23,24] and poses a number of practical issues, generally speaking. The question for instance arises whether to define contact forces at each contacting node or just at the one showing the greatest penetration depth, as discussed in [28]. While the consideration of every contacting node is essential to address non-convex cases where multiple contact points can be physically expected at a bodies pair, it may induce a magnification bias onto contact stiffness in case a tangent contact plane would intersect particle surface (and its node discretization) over a finite extent and/or when contact overlap would not be negligible, as addressed in [29]. In [29], a discretization-convergent solution was proposed via a continuum-based treatment of contacting nodes seen as contacting areas, while it has been proposed to replace nodes with either an application of minimization procedures on the SDF in [30] or a volume analysis of the overlap region in [31].

Restricting to the case of strictly convex particles (under DEM-usual stiff conditions), the definition of just one contact force at the most-penetrating node is however very much appropriate and discretization divergence such as discussed in [29] is no longer an issue. As such it is herein proposed for such cases to maintain in the Level Set approach a node discretization but to embed the latter with an innovative optimal surface coverage and new history-dependent properties for enabling more efficient LS-DEM simulations with reduced time costs. Applications are proposed for the case of superquadrics, as an example of strictly convex particles.

Section 2 first recalls the general definition and properties of superquadrics shapes and of a LS particle description, at the same time the latter is optimized for a faster initial generation and extended to a direct description of particle surface area. Section 3 then focuses on the surface nodes discretization of the LS approach, proposing optimized surface coverage and contact treatment for faster DEM simulations. Section 4 finally presents packing-scale simulations to quantify those time gains and illustrate how the method directly applies to studying the shear strength of packings of superquadrics. It also demonstrates on those packings how the LS approach outputs a consistently symmetric contact description in spite of its master-slave contact algorithm with different roles for bodies 1 and 2, as soon as appropriate surface nodes discretizations such as those herein proposed are used. All simulations

are run with the YADE code [32,33] whose source code is available at https://gitlab.com/yade-dev/trunk.<fnref n="1" />

## 2. Level set description and its application to superquadrics

### 2.1. General principles of a level set particle description in DEM

For the purpose of computing in DEM the relative movements of solid particles interacting through contact forces, it is recalled the LS approach combines two ingredients being defined at the scale of every particle [22,4,23].

First, a discrete set of values for the SDF $\phi(\bar{x})$ (with $\phi(\bar{x})$ being positive, resp. negative, for $\bar{x}$ outside, resp. inside of the particle) is defined on a particle-attached grid, after being for instance computed from a Fast Marching Method (FMM, [27]) which will be exposed in more details in the following Section 2.3.

Second, a set of so-called surface nodes discretize the particle surface itself and enable a precise contact detection and quantification, after a general broad phase of bounding boxes-based contact detection being systematic in YADE. Considering a pair of neighboring bodies as detected by the latter broad phase and denoting $S$, resp. $B$, the smaller, the resp. bigger, body in surface (extending [23] which was based on the volume), a master-slave algorithm actually checks every surface node $\bar{N}$ of $S$ in the distance field of $B$, $\phi_B$, to detect overlap through negative distance values. Interparticle overlap $u_n$ (chosen to be positive) and contact normal $\bar{n}$ are then deduced as follows:

$$u_n = - \min_{\bar{N} \in S} (\phi_B(\bar{N})) \qquad (1)$$

$$\bar{n} = \bar{\nabla} \phi_B \text{ at the corresponding } \bar{N} \qquad (2)$$

In Eq. (2), while a unit norm is directly obtained for an exact SDF, by virtue of the Eikonal equation $||\bar{\nabla}\phi|| = 1$, normalization is actually imposed in practice since the considered discrete distance fields are most often obtained from a numerical method such as the FMM (see the coming Section 2.3).

More details and new optimizations regarding the discrete SDF and surface nodes are discussed in the remainder of the manuscript. They will be illustrated on superquadric particles, whose definition is first recalled.

### 2.2. Superquadrics definition and properties

As proposed in [34,35], a superquadric, also known as superellipsoid, can be defined from five shape parameters, namely three “radii” $r_x, r_y, r_z > 0$ corresponding to the half-extents along the three axes of a local (inertial) frame and two exponents $\epsilon_e, \epsilon_n > 0$. The radii $r_x, r_y, r_z$, whenever different, classically distort shape length along the three local axes likewise to ellipsoids, while the exponent parameters $\epsilon_e, \epsilon_n$ enable one to extend beyond ellipsoids shapes, modifying cross-sections in planes including $\bar{z}$ or $(\bar{x}, \bar{y})$. Those parameters can be reduced to four whenever particle size does not matter after replacing $(r_x, r_y, r_z)$ with, e.g., $(r_x/r_z, r_y/r_z, 1)$. From these constant shape parameters and a pair of variable surface coordinates $(u, v)$ in $[-\pi; \pi] \times [-\pi/2; \pi/2]$ the following equations provide the locations of a surface point $(x, y, z)$ [35]:

$$x = r_x (\cos v)^{\epsilon_n} \text{sign}(\cos u) |\cos u|^{\epsilon_e} \qquad (3)$$

$$y = r_y (\cos v)^{\epsilon_n} \text{sign}(\sin u) |\sin u|^{\epsilon_e} \qquad (4)$$

$$z = r_z \text{sign}(\sin v) |\sin v|^{\epsilon_n} \qquad (5)$$

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

<footnote n="1">¹ Where examples/levelSet/lsNodeGeom.py illustrates the use of the §3.4 optimization, in particular.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Five possible superquadric shapes with parameters in Table 1.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Sharp variations of superquadric surface points (their (x, y) coordinates) as a function of (u, v) surface parameters in Eqs (3)-(5) for Shape B.</figcaption></figure>

<table><caption>**Table 1**

Superquadric parameters of the five shapes shown in Fig. 1.</caption><thead><tr><th rowspan="2">Shape</th><th colspan="3">Half-extents (length unit)</th><th colspan="2">Curvature exponents</th></tr><tr><th>r<sub>x</sub></th><th>r<sub>y</sub></th><th>r<sub>z</sub></th><th>ε<sub>e</sub></th><th>ε<sub>n</sub></th></tr></thead><tbody><tr><td>A</td><td>0.58</td><td>1</td><td>0.83</td><td>0.1</td><td>0.5</td></tr><tr><td>B</td><td>0.42</td><td>1</td><td>0.83</td><td>0.1</td><td>1</td></tr><tr><td>C</td><td colspan="3">same as Shape B</td><td>1</td><td>0.5</td></tr><tr><td>D</td><td>0.5</td><td>0.7</td><td>1</td><td>1.4</td><td>1.2</td></tr><tr><td>E</td><td>0.4</td><td>1</td><td>0.8</td><td>0.4</td><td>1.6</td></tr></tbody></table>

$$ \text{sign}(x) = \begin{cases} -1 & \text{if } x < 0 \\ 0 & \text{if } x = 0 \\ 1 & \text{if } x > 0 \end{cases} \tag{6} $$

which implies the following surface equation $f(x, y, z) = 0$:

$$ f(x, y, z) = \left( \left| \frac{x}{r_x} \right|^{\frac{2}{\epsilon_e}} + \left| \frac{y}{r_y} \right|^{\frac{2}{\epsilon_e}} \right)^{\frac{\epsilon_e}{\epsilon_n}} + \left| \frac{z}{r_z} \right|^{\frac{2}{\epsilon_n}} - 1 = 0 \tag{7} $$

Fig. 1 illustrates five possible convex superquadrics, with their parameters in Table 1.

The possibly sharp variations in the underlying analytical equations of superquadrics as functions of surface coordinates u, v, Eqs. (3)-(5), are illustrated in Fig. 2 with $\epsilon_e = 0.1$ and Shape B. Such variations certainly relate with possible convergence difficulties being often encountered while using an optimization approach based on Eq. (7) for computing the contact between superquadrics showing low values of $\epsilon_e, \epsilon_n$. As a matter of fact, optimization-based approaches were for instance carried out in, e.g., [19,20], adopting $\epsilon_e = \epsilon_n = 0.4$ in [19] or $\epsilon_e = \epsilon_n \ge 0.2$ in [20]. Here, the LS approach will be illustrated in the following paragraphs to be also somewhat adversely affected by extreme values for exponent shape parameters but to cope more easily with those since it accommodates $\epsilon_e$ values as low as 0.1 (Shapes A and B).

## 2.3. Level set particle description from signed distance function

Outside the very special case of degenerating into a sphere, superquadrics show no analytical expression of the shortest distance to their surface. As such, the distance field grounding the LS approach

is numerically obtained from the surface equation Eq. (7) and a Fast Marching Method (FMM, [27]), following [24]. The FMM solves the so-called Eikonal equation which is an inherent property of the SDF $\phi(\bar{x})$ as recalled in, e.g., [24]:

$$ ||\vec{\nabla}\phi|| = 1 \tag{8} $$

Eq. (8) actually serves as a non-linear PDE in space for $\phi$. It accepts as a boundary condition close to a superquadric surface the following distance approximation [24]:

$$ \phi(\bar{x}) = \frac{f(\bar{x})}{||\vec{\nabla}f(\bar{x})||} \quad \text{close to } \mathcal{S} \tag{9} $$

By construction, Eq. (9) indeed verifies distance properties since it obeys the Eikonal Eq. (8) close to the surface where $f(\bar{x}) = 0$, assuming a finite $\vec{\nabla}(1/||\vec{\nabla}f||)$. Considering a particle-attached discretization of the 3D space in the form of a rectilinear grid with an isotropic and uniform spacing $g_{grid}$, the FMM then provides distance values at every gridpoint after recursively propagating the boundary condition distance values along the surface normal, obeying Eq. (8).

Extending [24], a greater efficiency of the FMM is here obtained through the use of a heap-sort procedure for the required minimum search during recursive computations of distance that should apply to gridpoints with increasing distances. Reduced time costs are detailed in the Appendix A, with a possible speed-up greater than 100 with respect to the previous implementation of [24].

The particle inertial properties necessary for a DEM workflow directly derive from the distance field, e.g., the particle volume $V$ [4]:

$$ V = \iiint_{\text{Grid}} H(-\phi(\bar{x})) dV \tag{10} $$

with $H$ the Heaviside step function. In order to reduce discretization artifacts during the numerical estimation of Eq. (10) from discrete grid-point distance values $\phi_{ijk}$ and grid voxels volumes $g_{grid}^3$, the Heaviside step function is classically smoothened around $\phi = 0$ [4,31], i.e. within $[-\phi; \phi]$ with $\phi > 0$ a characteristic length, replacing Eq. (10) with:

$$ V = \iiint_{\text{Grid}} H_{sm} \left( \frac{-\phi(\bar{x})}{\phi} \right) dV \tag{11} $$
<figure><img src="image_6.png" /><figcaption>Fig. 3. Precision of the Level Set description of surface area S (for the superquadrics of Fig. 1) with respect to another estimation, $S^{ref}$, based on symbolic then numeric integration.</figcaption></figure>

while using $H_{sm}$ for a $C^2$ version of $H$:

$$H_{sm}(x) = \begin{cases} 0 & \text{if } x < -1 \\ 1 & \text{if } x > -1 \\ \frac{1}{2} \left( 1 + x + \frac{\sin(\pi x)}{\pi} \right) & \text{if } x \in [-1; 1] \end{cases} \qquad (12)$$

Considering that Eq. (11) is readily generalized to compute any volume included within a distance threshold $\phi = \epsilon$ (in addition to the particle volume corresponding to $\epsilon = 0$) as

$$V(\phi < \epsilon) = \iiint_{\text{Grid}} H_{sm} \left( \frac{\epsilon - \phi}{\varphi} \right) dV, \qquad (13)$$

it is herein proposed to make further use of the smooth function $H_{sm}$ to differentiate the above equation for a direct measure of the Level Set particle surface, $S$:

$$S = \lim_{\epsilon \to 0} \frac{V(\phi < \epsilon) - V(\phi < -\epsilon)}{2\epsilon} \qquad (14)$$

Considering the case of a sphere of radius $R$ with reference surface $S^{th} = 4\pi R^2$, it is logically observed that choosing $\varphi$ in the order of grid spacing, i.e. the characteristic length of discretization artifacts, gives satisfactory results where the relative error on $S$ versus $S^{th}$ drops at $10^{-3}$ for a grid resolution $2R/g_{grid} = 20$ and $\varphi = \sqrt{3}g_{grid}/2$ (a half-diagonal of grid cell) whereas $\varphi = 5\sqrt{3}g_{grid}/2$ or $\varphi = \sqrt{3}g_{grid}/10$ would keep that relative error above $10^{-2}$ or even $10^{-1}$. Other validations of Eq. (14) and of the LS particle description in general are proposed in the next subsection, with $\varphi = \sqrt{3}g_{grid}/2$ being now chosen for the remainder of the manuscript.

## 2.4. Particle-scale validation of the level set description

The precision of the present LS approach for particle description is naturally ruled by the grid spacing underpinning the above discrete distance field, for several reasons. First, $g_{grid}$ directly acts a spatial mesh size during the FMM for obtaining distance values, in a first-order fashion which has been discussed in [24]. Second, the grid spacing naturally affects the voxel-based estimation of volume and surface, Eqs. (11) and (14).

Regarding volume, the precision of the LS approach has previously been discussed in [24], where it was shown that using several tens of voxels per particle length led to satisfactory results, even with the non-smoothened Eq. (10).

Regarding surface, the Fig. 3 analyzes the grid resolution influence onto the LS evaluation in Eq. (14) for the superquadrics at hand

<caption>Table 2</caption>

Reference values of particle properties for the five superquadric shapes shown in Fig. 1. Surface area is computed as per the Appendix B, volume from exact formula [35] and Wadell's true sphericity [36] is deduced as $\psi = \sqrt{36\pi V^{ref}^2}/S^{ref}$.

<table><thead><tr><th>Shape</th><th>Surface area $S^{ref}$ (length<sup>2</sup> unit)</th><th>Volume $V^{ref}$ (length<sup>3</sup> unit)</th><th>True sphericity $\psi$</th></tr></thead><tbody><tr><td>A</td><td>12.228</td><td>3.353</td><td>0.886</td></tr><tr><td>B</td><td>8.501</td><td>1.852</td><td>0.858</td></tr><tr><td>C</td><td>8.598</td><td>1.914</td><td>0.867</td></tr><tr><td>D</td><td>5.636</td><td>1.093</td><td>0.910</td></tr><tr><td>E</td><td>5.970</td><td>1.086</td><td>0.856</td></tr></tbody></table>

(Fig. 1). Obtained surface area values are compared with another estimation based on symbolic then numeric, integration of surface equations (3)-(5), using `sympy` and `scipy` Python modules, that serves as a reference and which is itself precised and validated in Appendix B. Table 2 groups the corresponding reference values $S^{ref}$, with also the values for volume and Wadell's true sphericity [36] as a classical shape index, e.g. [37]. While being shape-dependent and not monotonously increasing with grid resolution, the precision in surface area also appears in the Fig. 3 as satisfactory for grids with tens of voxels per (smallest) particle length, i.e. $2\min(r_x, r_y, r_z)/g_{grid} \ge 20$, with a corresponding error not exceeding a couple of percent for the worst case (Shape A) and being much smaller in the other cases.

It is also to recall that surface nodes are located at $\phi = 0$ from a (tri-)linear interpolation between gridpoints of the SDF $\phi$ which is naturally again influenced by $g_{grid}$, similar to the estimation of normal vectors from the SDF-gradient.

The Fig. 4 then illustrates the spacing-dependent precision of the LS shape description for superquadric A (see also Appendix C for the other shapes), first for what concerns the position of surface nodes (Fig. 4a), as quantified looking at the corresponding value of the surface function $f(x, y, z)$, Eq. (7), which should be theoretically 0. Second, the precision in the normal vector at the same locations is also compared with an analytical reference given in [35] as a function of $u$ and $v$ surface parameters (Fig. 4b). It is to note such an analysis of normal orientation is somewhat biased by the impossibility to rigorously determine a $(u, v)$ pair for a $(x, y, z)$ point that is not perfectly located over the superquadric surface, when inverting Eqs. (3)-(5) and considering also the abrupt variations in, e.g., $x(u, v)$ shown in Fig. 2.

Anyway, the consideration of several tens of grid voxels per particle axis is here again shown to provide a satisfactory LS description of superquadric shapes, whereby using $g_{grid} = 2\min(r_x, r_y, r_z)/40$ induces a negligible error both on surface positions, with $f < 0.01$ for all cases, and surface normals, with 75% of surface points showing an orientation mismatch of less than one degree for normal vectors, with an upper bound of around 10 degrees for all cases. From a computational perspective, this corresponds to a total number of gridpoints between $1.9 \times 10^5$ and $3.7 \times 10^5$, depending on the considered shape and its aspect ratio, i.e. a minimum RAM requirement around 2 MB per particle if one just considers the discrete distance data, when implemented in a double floating-point precision.

One can also note the case of $g_{grid} = 2\min(r_x, r_y, r_z)/20$ which could serve as a good compromise for a (eight times) lower memory requirement since the latter induces median errors less than 0.01 for position and less than the degree for orientation.

## 3. History-dependent optimized contact treatment with a spatially uniform surface discretization

### 3.1. Definitions

In [23,24], surface nodes were computed through ray tracing in the above distance field, searching for $\phi = 0$ along a number of rays defined by their spherical angles $(\theta, \varphi) \in [0; \pi] \times [0; 2\pi]$, with $\theta$ being measured

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132
<figure><img src="image_7.png" /><figcaption>Fig. 4. Statistical properties of the errors in surface position or surface normal of the Level Set description of Shape A (other shapes are shown in Appendix C). Minimum values are outside the scales and do not appear on the graphs. Data have been obtained considering 4902 surface nodes as measurement points, being distributed as discussed in §3.3. Other kinds of precise surface sampling show the same trends, with some differences in extreme values.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Different possible distributions of points on a unit sphere (1602 here in both cases) and surface nodes spherical coordinates.</figcaption></figure>

from the $\vec{z}$ axis and $\varphi$ from $\vec{x}$ in the $(\vec{x}, \vec{y})$ plane (in a body-attached inertial frame). Logically, the number of surface nodes, $N_{nod}$, is a key parameter for the method precision and its computational cost [23,24]. It actually also combines with the spatial distribution of surface nodes since the latter should ideally induce an uniform surface coverage by each node, for a given total number. The possibility for such an ideal distribution is discussed in the following paragraphs, starting with the case of a LS description of a spherical shape as an idealized example with available theoretical results.

## 3.2. Optimizing the spatial distribution of nodes on a sphere

From a theoretical point of view, no perfectly uniform distribution can actually be drawn on a sphere beyond $N_{nod} = 12$ surface points, where the icosahedron as the fifth and last Platonic solid, showing 12 vertices and 20 equilateral triangular faces, would approximate a sphere as a regular mesh with uniform surface tiles and inter-vertex distances, as ideally desired from a LS-DEM point of view. For a greater $N_{nod}$, no regular mesh i.e. regular convex polyedra can be created strictly speaking, as imposed from a connectivity point of view and the Descartes-Euler polyhedral formula.

Approximating a uniform mesh on the sphere could still be realized by refining the icosahedron, introducing additional vertices at the middle of each edge followed by projections on the sphere, e.g. [38]. Such a construction nevertheless does not really enable one to control the number of obtained nodes as vertices, outside of choosing a refining level $k$ that will lead to $4^k \times 10 + 2$ vertices, e.g., resp. 42, 162, 642, 2562 and 10242 after refining resp. 1, 2, 3, 4 and 5 times.

Another interesting distribution of surface nodes, in an arbitrary $N_{nod}$ number and otherwise similar to the refined icosahedron (see below), can be obtained along a spiral path such as proposed in [39] and previously used in [23] i.e. along the following spherical coordinates $(\theta_i, \varphi_i)$, $i \in [0; N_{nod} - 1]$:

$$ \theta_i = \arccos \left( -1 + \frac{1+2i}{N_{nod}} \right) \qquad (15) $$

$$ \varphi_i = \pi(3 - \sqrt{5})i, \qquad (16) $$

as illustrated in Fig. 5a.

As it can be seen in Fig. 5, the above spiral distribution avoids to over-discretize in terms of nodes surface density towards the poles relative to the rest of the surface, as it would happen with the simplest distribution being implemented along meridians and parallels, i.e. through a rectangular partition of the spherical angle space $(\theta, \varphi) = [0; \pi] \times [0; 2\pi]$. For such a rectangular distribution, $N_{nod} = N_a^2 + 2$ surface nodes would be located over all possible pairs built from the $N_a \in \mathbb{N}^*$ different values of $\theta$, resp. $\varphi$, as per the following Eq. (17), resp. (18):

$$ \theta_i = i \frac{\pi}{N_a + 1} ; \quad i \in [1; N_a] \qquad (17) $$

$$ \varphi_j = j \frac{2\pi}{N_a} ; \quad j \in [0; N_a - 1] \qquad (18) $$

complemented with two nodes at the poles along the $\vec{z}$ axis:

$$ \theta = 0[\pi] \qquad (19) $$

The impact of such a choice of distribution onto the LS-DEM precision can be observed considering the reconstruction of a dense packing of 8000 spherical particles submitted to an isotropic stress $p_{ref} = 16.5$ kPa, following [23] where only the spiral distribution was considered.
<figure><img src="image_9.png" /><figcaption>Fig. 6. Role of the spatial distribution of surface nodes when reconstructing with the Level Set approach a packing of spheres under an isotropic stress $p_{ref}$, where the construction of a “Rectangular optimized” distribution is described in §3.3.</figcaption></figure>

Combined with contact stiffness $k_n = 6 \times 10^5$ N/m and mean particle diameter $D_{50} = 0.122$ m, such a stress state corresponds to a dimensionless stiffness ratio [40], basically the inverse of dimensionless interparticle overlap, $\kappa = k_n/(p_{ref} D_{50}) \approx 300$, which is consistent with the range of usual DEM applications [23]. The considered reconstruction namely uses known particles positions and diameters in a reference state obtained with classical DEM, in order to redefine them in a LS fashion with a very fine grid including 90 voxels per each particle diameter (that applies to all results discussed in this section) and a variable number of surface nodes being either distributed along the vertices of the refined icosahedron, or the spiral path of Eqs. (15)-(16) or in the rectangular fashion of Eqs. (17)-(19).

Fig. 6 illustrates the induced consequences on the obtained precision when re-calculating the obtained mean stress $p$ in a LS approach. As expected, the icosahedron construction appears as the most efficient for capturing the highest ratio $p/p_{ref}$, i.e. the highest proportion of expected contacts and overlaps. In addition, the spiral distribution is virtually as efficient, while also offering a precise control on the number of nodes. Quantitatively speaking, the spiral distribution enables one to achieve a same 85% or 90% precision (in terms of $p/p_{ref}$) with 30% less surface nodes with respect to the rectangular case, i.e. a time cost possibly 30% smaller since the latter is mostly proportional to the number of surface nodes (for such cases with a large nodes number) [24].

The least optimal nature of the rectangular path versus the spiral one (or the refined icosahedron distribution) is interpreted, similar to [41], by considering that the distribution $\{\min_j(d_{ij})\}$ of the closest distance $d_{ij}$ between a surface node $i$ and any other node $j$ shows the largest width for the former, with a number of surface nodes being significantly more isolated than others (Fig. 7).

Based on the same concept of internode distance as a metric for surface coverage (a choice which is further justified with respect to another possible metric in Appendix D), it is then to note how the same spiral path for surface nodes would induce widely different internode distances for more complex surfaces, e.g. Shapes A and E (Fig. 8), at the detriment of a LS-DEM application. This shows that another strategy, being more robust with respect to generic shapes, should be employed to attain almost uniformly distributed points.

### 3.3. Optimizing nodes distribution in the general case

An optimization procedure is then herein proposed whereby an initial spatial distribution of surface nodes (e.g., the rectangular one) is disturbed according to a fictitious repulsive force field between neighbors. Nodes being seen as particles of unit mass $m = 1$ (any consistent system of units is appropriate in this section), repulsive forces $\vec{F}_{ij}$, as exerted by a node $i$ onto a node $j$, apply with an electrostatic-like nature:

<figure><img src="image_10.png" /><figcaption>Fig. 7. Distributions of the smallest distance between a surface node and every other, for various choices of discretizations on a sphere surface. The Figure is constructed considering between 2502 (rectangular path) and 2562 (other cases) nodes and it is checked other nodes numbers show the same trends (see also Fig. D.23a).</figcaption></figure>

$$ \vec{F}_{ij} = \frac{k}{d_{ij}^3} \vec{d}_{ij} \qquad (20) $$

where $\vec{d}_{ij}$ is the branch vector from node $i$ to node $j$, of norm $d_{ij}$, and a unit $k$ is here considered. The consideration of other values for the stiffness coefficient $k$ or the mass $m$ would just induce a modified time scale for the procedure and would not change the finally obtained distribution. In order to limit the computational costs of the latter by restricting the number of interactions to the most influential ones (with the shortest distances), such a repulsive force field is considered to exist up to the arbitrary threshold of $d_{ij} = 5 \frac{L}{\sqrt{N_{nod}}} \quad \text{with } L = V^{1/3}$ a characteristic length of a particle of volume $V$.

DEM time integrations are then executed to update nodes positions, combining the effects of the repulsive interactions with additional backprojections of the nodes onto the surface at every time step, the latter being straightforward to apply thanks to the Level Set surface description.

While the principle of the procedure is to search for more uniform internodes distance through a balance of the distance-dependent repulsive forces, no convergence can actually be expected to a well-defined configuration with perfectly uniform internode distances, thinking that such an ideal case is even impossible to reach for the simplest case of a sphere beyond $N_{nod} = 12$, as mentioned in the above. As also rigorously analyzed in [42] in the context of minimal energy point distributions on manifolds, there is no stable equilibrium state for a large number of points subject to an antigravitational interaction on the sphere where any arrangement tends to be unstable, leading to persistent or divergent oscillations. For this reason the choice is made to just search for a local minimum of the distribution width of $\{\min_j(d_{ij})\}$, which is performed as follows.

After running a maximum of 1500 iterations while using a non-viscous “Cundall” damping [43,24] with coefficient $D = 0.2$ or $0.5$ (all present choices being arbitrary), the procedure enables one to explore a large variety of nodes distributions, whose most optimal one is eventually selected in a post-processing stage as the one showing the minimum width for the $\{\min_j(d_{ij})\}$ distribution, among all distributions obtained throughout the DEM iterations of every case.

Connected to the lack of convergence, optimized node distributions are not deemed to be unique. The fact that the dynamic behavior of the procedure is observed to be rather chaotic (with a salient sensitivity to initial conditions) also contributes to the non-unicity but it is anyway observed that different spatial distributions of nodes with similar optimized spreads of the $\{\min_j(d_{ij})\}$ set give very similar DEM results, as
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Spatial distribution along the spiral path of Eqs (15)-(16) of 1602 surface nodes over two possible superquadric shapes (with reference to Table 1), with an apparent non-uniform internode distance.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Surface nodes on a sphere as obtained after applying the repulsive-field-based optimization procedure to an initial rectangular distribution of 1602 nodes.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Spatial distribution of 1602 surface nodes over two possible superquadric shapes (see Table 1), as optimized with respect to Fig. 8 thanks to the herein proposed procedure.</figcaption></figure>

already suggested by the similar macro-scale performances of the different spiral and icosahedral distributions on the Fig. 6.

The outcome of such a procedure applied to a Level Set sphere with 1602 surface nodes initially distributed along the rectangular path is illustrated in Fig. 9. As it can be seen in Fig. 9a in comparison to Fig. 5b, the pole overdiscretization together with the initial rectangular pattern are clearly erased by the procedure. The distribution of shortest internode distance is also efficiently made much more compact, getting closer to the spiral distribution, see Fig. 9b vs 7 or D.23a. As a result, the new distribution is now virtually as efficient as the spiral one for capturing the reference mean stress during the isotropic reconstruction (Fig. 6), which demonstrates the interest of such an optimization procedure.

Its applicability to the case of superquadrics is illustrated in Fig. 10 which indeed displays more uniform internode distances with respect to previous Fig. 8, see also Fig. 11 for a quantitative analysis of the internode distance distributions.

Since the optimal coverage of a surface using points is of interest to countless scientific applications, be it for serving as a spatial discretization for meshless numerical methods or in computer graphics, a number of algorithms sharing the same objective have been previously proposed, starting with [44] for the spherical surface while considering similar repulsive interactions between discretization points. Numerous studies such as [45,46] have since been devoted to triangulated representations of generic surfaces, possibly helping to turn those meshes into a collection of equilateral triangles with uniform inter-vertices distances. They, however, inherently require connectivity information which are absent from the present points discretization of a LS shape. On the other hand, some more recent approaches directly handle a raw collection of points, e.g., [47,41,48]. In [41], the method builds upon a 2D parametrization of the surface to discretize and corresponding mappings are again absent for a general LS shape description, complicating any further comparison. The studies of [47,48] resemble much more closely the present approach, where discretization points are displaced towards a uniform
<figure><img src="image_14.png" /><figcaption>(a) Shape A before optimization</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Shape E before optimization</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) Shape A after optimization</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) Shape E after optimization</figcaption></figure>

<caption>Fig. 11. Distribution of internode distances in the cases shown in Figs. 8 and 10.</caption>

distribution based on repulsive interactions. In both those previous studies, the displacement algorithm is nevertheless unphysical since it does not follow dynamics law, making it for instance difficult to control the stability of their explicit algorithms. Looking, e.g., at the Shape E, the algorithm of [47], in its default formulation, is indeed observed to diverge when trying to build a distribution of 1484 surface nodes while successfully distributing 1596 nodes with a very similar span of {$min_j(d_{ij})$} to what is obtained here. As a matter of fact, the design of the present approach within a physical DEM framework gives both advantages of making it applicable in any DEM implementation with a Level Set shape description and of offering an easy control of its stability from classical considerations of a DEM critical time step. Namely, the iterations are run

with the variable, but automatically set, time step $\Delta t = 0.05\sqrt{\frac{md_{min}^3}{k}}$ that is obtained through the tangent stiffness of Eq. (20) as a function of the absolute minimum of internodes distance, $d_{min}$, and is small enough to avoid divergence of the explicit scheme.

### 3.4. History-dependent contact detection

Moreover the above optimization of the spatial distribution of nodes, a faster execution of the LS contact treatment for a given particle pair is also herein proposed, with a similar idea to the minimization procedure of [20,30]. In these studies, the contact point location at a previous DEM iteration serves as an initial guess in the iterative minimization procedure of distance that should eventually give the contact point at the current DEM iteration, leading to a reduced number of minimization iterations whenever past and present contact points are reasonably close to each other. In the present node-based implementation, a similar time-saving strategy is employed by limiting the loop at each LS-DEM iteration to a smaller subset of nodes. Instead of iterating over all surface nodes, the algorithm first focuses on the $N_n$ closest neighbors of the node previously identified as having the minimum interparticle distance, along with that node itself. This optimization is applied even in

the absence of contact such as when particles are separated by a small gap to minimize unnecessary looping over all nodes whenever possible.

Fig. 12 details the underlying procedure as a flowchart. In the most ideal case of the same surface node serving as contact point during several iterations, the procedure replaces a loop over a total number of $N_{nod}$ points with just one loop over $N_n + 1$ points. However, even in such a case, no $N_{nod}/(N_n + 1)$ speed-up can be expected for various reasons. First, the linear algorithmic complexity of contact treatment is lost for a small number of nodes [24], due to additional pre- and post-processing operations to the loop, e.g., computing and assigning contact normal in the current configuration once the contact-carrying node has been determined. Second, during the full consideration of all nodes, efficient operations prompting whether a given node of 2 enters the distance grid of 1 lead to quickly disregard a number of them being far away from the contact region, before actually computing any distance value. Since such a situation does not happen when considering just $N_n$ neighbors being in the contact region, the time cost per considered-node is necessarily higher in the latter case, preventing a $N_{nod}/(N_n + 1)$ speed-up ratio.

Thanks to the recursive consideration of possibly different subsets of neighbors, the procedure stays anyway robust for capturing the contact point also in case of important relative changes in the pair configuration, even though time gains are smaller in this case.

### 4. Case study: REV

#### 4.1. Computational comparison with classical DEM for spherical particles

The computational impact of the history-dependent contact treatment just proposed in §3.4 is quantified in the Fig. 13 for the realistic setting, previously considered in [23], of a quasi-static compression of an initially dense packing of 8,000 spherical particles for which a classical DEM approach can readily be executed. The example thus provides a reliable benchmark for comparing execution times of the classical DEM approach with those of LS-based methods employing various values for $N_{nod}$ and $N_n$.
<figure><img src="image_18.png" /><figcaption>Fig. 12. Flow chart for a history-dependent LS-DEM contact treatment.</figcaption></figure>

In more details, wall clock execution times *t* are measured when using a server machine with two Intel Xeon Platinum 8270, 2.7 GHz, processors with 26 cores, 52 threads, and 36 MB of cache memory each, together with 1.5 TB 2.9 GHz RAM. In order to account for the inherent variability of execution times, each LS case is executed five times in slightly different usage conditions of the server machine (in terms of total number of simulations running simultaneously, i.e. with different total memory requirements and cache access possibilities), while the classical DEM case is run twice, in a sequential execution for all of them. From the ten possible ratios $t(\text{LS-DEM})/t(\text{DEM})$, the Fig. 13 depicts the corresponding mean together with the standard deviation as an error bar, that should give an estimate, by default, of time variability that would arise from the cache memory access considerations, for instance. The Fig. 13 first recalls, as discussed in [23], how the versatile LS approach logically induces much higher time costs than classical DEM for these simple spherical shapes, in direct connection with the chosen total number of surface nodes $N_{\text{nod}}$. Even though an influence of $N_{\text{nod}}$ remains for such a realistic, packing-scale, simulation (unlike the ideal case of a fixed particle pair), it is also illustrated how the present optimization enables one to substantially decrease those overheads, where relative gains are barely affected by the chosen $N_n$ and the most interesting for a high $N_{\text{nod}}$, possibly dividing time costs by more than 3. As already mentioned in the §3.4, this speed-up ratio can not be expected to be in the order of $(N_{\text{nod}} + 1)/N_n$ since the Fig. 13b also shows how LS-DEM time costs, being optimized or not, include a constant (independent with respect to the total number of nodes and fairly common to the two LS-based variants) overhead with respect to classical DEM which is expected to come in particular from the contact post-processing operations and the higher memory footprint of a LS simulation that will affect cache misses.

<caption>Table 3</caption>

Obtained Cundall number (in DE.iterations/s) during sequential executions of different LS simulations (each case being repeated three times before an interval of extreme Cu values is formed) on the server hardware described in main text.

<table><thead><tr><th></th><th>Classical LS approach [23,24]</th><th>Present Optimized LS approach</th></tr></thead><tbody><tr><td>Quasi-static compression of a packing of 8,000 spheres with 2,500 surface nodes each</td><td>∈ [1785; 1790]</td><td>∈ [5540; 5575]</td></tr><tr><td>Dynamic discharge of 1,000 superquadrics with 2,000 surface nodes each</td><td>∈ [1955; 1972] [24]</td><td>∈ [2530; 2566]</td></tr></tbody></table>

Being significant anyway, the time gains of the present approach are accompanied with some additional memory (RAM) requirements for the purpose of storing lists of neighboring nodes, that adds to the LS required memory for storing the distance grids and the surface nodes of all particles, as discussed in [23,24]. In more details, the additional memory requirements increase with both $N_n$ and $N_{\text{nod}}$ but correspond to just an additional 13% for $N_n = 6$ and the precise case of $N_{\text{nod}} = 4000$ (with the present grid resolution of $D/g_{\text{grid}} = 20$).

Table 3 summarizes obtained time gains in the form of a Cundall Number Cu being defined from execution times *t* together with $N_{\text{DE}}$ as the total number of Discrete Elements (DE) handled in a given simulation, where $N_{\text{it}}$ DEM iterations have been executed:

$$Cu = \frac{N_{\text{DE}} N_{\text{it}}}{t} \qquad (21)$$

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132
<figure><img src="image_19.png" /><figcaption>(a) For various $N_n$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) For various $N_{nod}$</figcaption></figure>

<caption>**Fig. 13.** Reduced time overhead of the Level Set approach, relative to classical DEM, thanks to the herein proposed history-dependent contact treatment, unless $N_n = 0$ which stands for the classical Level Set approach; in the case of a compression of a packing of spheres and sequential executions.</caption>

In addition to the previous quasi-static compression of a tightly coordinated (at initial state, in particular) packing of 8,000 spheres where Cu can be brought up to approx. $5.5 \times 10^3$, the Table 3 also considers a more dynamic simulation of a discharge of 1,000 superellipsoids (of Shapes A to E in equal proportion) into a container, including a first phase of free fall of a cloud of particles, see [24] for a more detailed description. Obtained performance gains, up to $Cu \approx 2.5 \times 10^3$ are smaller for that second simulation, in line with more pronounced relative changes in pairs configurations during most of the discharge simulation.

In comparison with other previous studies based on optimization approaches instead of a Level Set one, such values are in an intermediate range. Indeed, it has been (either directly or indirectly) reported by resp. [49,50] and [20] values of $Cu \approx 2 \times 10^2$ and $Cu \approx 2 \times 10^5$, when using slightly faster processors with a base frequency of resp. 2.9 and 3.6 GHz for simulating resp. convex NURBS and convex superquadrics shapes, where it is recalled that the present LS approach also applies to every convex shape, beyond those cases with an analytical definition.

## 4.2. Residual shear strength of superquadrics

As an example of possible applications, the proposed implementation is applied for a preliminary study of the shear strength of mono-shaped granular packings, i.e. packings where particles all conform the same superquadric shape among Table 1. Particles are packed within frictionless wall boundaries, making the (X, Y, Z) axes coincide with principal directions of stress or fabric and their macroscopic shear strength is char-

<caption>**Table 4**</caption>

DEM parameters for the axisymmetric compressions on superquadrics packings.

<table><thead><tr><th colspan="3">Contact properties</th><th>Density</th><th>Timestep</th><th>Damping</th></tr><tr><th>$k_n$ (N/m)</th><th>$k_t/k_n$ (-)</th><th>$\mu$ (-)</th><th>$\rho$ (kg/m³)</th><th>$\Delta t$ (µs)</th><th>D (-)</th></tr></thead><tbody><tr><td>$5 \times 10^5$</td><td>0.8</td><td>$\tan(30^\circ)$ or 0 (surrounding walls)</td><td>2600</td><td>$0.4 \sqrt{\frac{m_{min}}{k_n}}$ $\in [9.6; 13.2]$</td><td>0.3</td></tr></tbody></table>

acterized during axisymmetric, so-called “triaxial”, compressions that are carried out as follows.

Starting from an initial cloud of particles (in a variable number $N_{DE}$ across different packings) with no contact in between (see Fig. 14a), an isotropic compression is applied until the stress state stabilizes at a $p = 100$ kPa isotropic state (Fig. 14b), which incidentally requires more iterations for Shapes A and B because of their sharper surface variations, with the lowest values of $\epsilon_e$ exponent in Table 1. Then, a shortening is imposed in the axial direction at a constant rate $\dot{\epsilon}_1 = \dot{\epsilon}_{YY} > 0$ (soil mechanics sign convention being applied here with positive compressive strains or stresses) while maintaining a constant confining stress $\sigma_3 = \sigma_{XX} = \sigma_{ZZ} = 100$ kPa in the radial directions (Fig. 14c). Aiming to characterize a behavior that would be close to quasi-static, the strain rate is chosen to correspond to an inertial number [3] $I = \dot{\epsilon}_1 D_{50} \sqrt{\rho/\sigma_3} = 10^{-4}$. While $I = 10^{-3}$ has sometimes been considered to be small enough for quasi-static conditions to hold [3,6,51], some non-negligible variations in the stress-strain response are here observed between $I = 10^{-4}$ (smaller values being not considered for computational cost reasons) and $10^{-3}$, see Fig. 15.

In terms of particle size, the latter conforms a narrow uniform (in count) distribution in [9.5 mm; 10.5 mm], i.e. around $D_{50} = 10$ mm, after being defined for those non-spherical shapes in a sieve-inspired manner, choosing as a particle characteristic length twice the greatest center-to-surface distance. That length is measured once for all outside of the triaxial simulations on a exhaustive collection of 4902 surface nodes distributed as per the §3.3 procedure. During the triaxial simulations themselves, the LS particles of Shapes A and B (with $\epsilon$ exponents as low as 0.1) require the same discretization with 4902 surface nodes for the packings and their contact networks to stabilize under the pre-shear $p = 100$ kPa isotropic stress state, while Shapes C to E (with $\epsilon$ exponents being greater than 0.5) accommodate using 3602 surface nodes being again distributed according to the §3.3.

Contact parameters, i.e. two linear stiffnesses $k_n, k_t$ and a friction coefficient $\mu$ as exactly defined in [24], are listed in Table 4. Restricting to a qualitative study with no experimental reference to match, the normal stiffness $k_n$ is just chosen such that inter-particle dimensionless overlaps are negligible as expected in DEM, with their inverse ranging in the order of $\kappa = k_n/(pD_{50}) = 500$, considering likewise to [52] that the later quantity applies with this meaning both to spherical [40] and non-spherical shapes. As for the tangential stiffness, the corresponding $k_t/k_n = 0.8$ ratio would correspond to a classical macroscopic Poisson's ratio around 0.3 in the initial elastic regime, according to [53].

Table 4 also includes other DEM parameters with a non-viscous damping coefficient $D$ that introduces additional dissipative forces in Newton-Euler equations of motion [43,24] and a time step estimation directly based on maximum contact stiffness and the minimum particle mass within a packing, $m_{min}$.

Towards the end of performing such a triaxial loading up to $\epsilon_1 = 0.5$, a classical critical state [54] is obtained with constant deviatoric stress $q = \sigma_1 - \sigma_3$ and porosity $n$, even though $\epsilon_1$ keeps increasing, where the final values of $q$ and $n$ do not depend on initial packing porosity (Fig. 16).

Those critical states are herein briefly analyzed in connection with particle shape for what concerns the residual shear strength of the material, i.e. the critical state value $\eta_{crit}$ of $\eta = q/p$. In order to account for the inherent fluctuations of the stress-strain curves for such granu-

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132
<figure><img src="image_21.png" /><figcaption>(a) Initial cloud of $N_{DE} = 1000$ particles (b) Compacted packing of $N_{DE} = 8000$ particles of Shape A under $p = 100$ kPa</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(c) Compressed packing of $N_{DE} = 16000$ particles of Shape D</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(a) Shape A</figcaption></figure>

<caption>(b) Shape E</caption>

<caption>Fig. 15. Influence of inertial number $I$ during axisymmetric compressions under $\sigma_3 = 100$ kPa for packing of 8000 particles of different shapes.</caption>

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132
<figure><img src="image_24.png" /><figcaption>Fig. 16. Selected macroscopic responses to the triaxial loadings on different shapes (in packings of 8000 particles with either dense or loose initial states).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 17. Residual shear strength $\eta_{crit}$ (its mean and extreme values on $\epsilon_1 \in [0.475; 0.5]$) for various packings with different numbers of particles of the same shape. Thicker symbols and error bars (only present for 1000, 8000 and 16000 particles to reduce the total number of simulations) correspond to packings with a same number of particles but a looser initial state; with data sets of those loose and dense packings being strictly superposed, as it is physically expected, in a number of cases (e.g., Shape C for 8000 particles).</figcaption></figure>

lar packings, in particular for a small number of particles $N_{DE}$, $\eta_{crit}$ is measured as an average over a final interval of $\epsilon_1 \in [0.475; 0.5]$ (that includes 10 measured datapoints) with also error bars denoting its extreme values on that interval, and depicted as such in Fig. 17 for the various shapes and $N_{DE}$.

Various values for $N_{DE}$ are tested since the appropriate number of superquadrics necessary to form a Representative Elementary Volume (REV) under rigid boundaries and obtain meaningful macroscopic quantities still seems to be a matter of choice, with $N_{DE}$ being chosen as 5000 in [55], 10000 in [51], or defined through convex sub-samples being extracted from a similarly cuboid packing of 1609 particles in [56]. One may actually expect the appropriate number of $N_{DE}$ to be shape-dependent, as suggested by the Fig. 14b where Shape A for instance appears to be more prone to boundary effects in terms of a favored orientation when touching the boundary walls.

Selecting for the scope of the present study the case of $N_{DE} = 8000$ that corresponds to a limited computational cost and consistent values of $\eta_{crit}$ (e.g., between the dense and loose cases) for a given shape, the Fig. 18 summarizes the critical shear strength of those mono-shaped

packings, in comparison with additional data for the case of a packing made of spheres which is sheared in the same conditions (same particle size distribution and other parameters of Table 4). Comparing with the sphere case, a obviously higher shear strength of superquadric particles appears. Expressing the latter through a Mohr-Coulomb internal friction angle $\phi_{crit} = \sin^{-1}(3\eta_{crit}/(6+\eta_{crit}))$, $\phi_{crit}$, the packing of spheres indeed shows $\phi_{crit} \approx 19^{\circ}$ while superquadric-based packings reach $\phi_{crit}$ values between 26$^{\circ}$ and 33$^{\circ}$. Within the present set of superquadric shapes, no exact decreasing trend of $\eta_{crit}$ or $\phi_{crit}$ with respect to $\psi$ can be observed though, prompting the need for further studies on the physics of granular packings with the present LS approach whose physical consistency is finally discussed in terms of objectivity.

### 4.3. Objectivity of the contact algorithm

By its master-slave nature appearing in Eqs. (1)-(2), the LS contact algorithm for computing an overlap and contact normal can not automatically entail the physically expected symmetry with respect to a possible swap in the choice of the S and B particles (e.g., for equal-sized

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132
<figure><img src="image_26.png" /><figcaption>Fig. 18. Obtained residual shear strength according to grain shape sphericity for the present mono-shaped packings of 8000 particles.</figcaption></figure>

pairs). The Eqs. (1)-(2) in particular do not directly imply the existence of a common normal which shall exist when the present rounded particles contact themselves along a unique tangent plane and which actually grounds other contact detection algorithms based on the idea of searching a common normal within the two bodies surfaces, e.g., [56,57].

However, one may again expect that a sufficient discretization in terms of LS surface nodes on strictly convex particles would lead to localize the single-point contact on a correct, unique, tangent plane and to indirectly obtain a symmetric contact description. Such an idea is investigated for some of the superquadric packings that were analyzed in the previous sections, namely those including 8000 particles of Shapes A or E in the pre-shear isotropic stress state of $p = 100$ kPa, but using now a variable number of nodes for each case.

The possible asymmetry of the LS contact algorithm is first quantified regarding the overlap $u_n$, where a dimensionless global error metric $\Delta^*(u_n)$ is defined from the average difference between the overlap $u_n$ deduced from Eq. (1) and another overlap $u'_n$ that would have been computed after exchanging the role of S and B in Eq. (1):

$$u'_n = -\min_{\vec{N} \in B} (\phi_S(\vec{N})) \qquad (22)$$

$$\Delta^*(u_n) = \frac{1}{u_n^\kappa N_c} \sum_c |u_n - u'_n| \qquad (23)$$

$$u_n^\kappa = \frac{p D_{50}^2}{k_n} \qquad (24)$$

In Eqs. (23)-(24), that error metric is made dimensionless after using the order of magnitude of interparticle overlap in a simulation, $u_n^\kappa$, instead of using relative differences between $u_n$ and $u'_n$ (with respect to, e.g., $u_n$), that would negate the physically greater role of contacts with a great overlap $u_n$, i.e., normal force [58]. By convention, Eq. (23) applies to the set of contacts $c$ (with $N_c$ for its total number) where either $u_n$ or $u'_n$ is positive, such that inconvenient cases where, e.g., Eq. (1) but not (22) would detect contact, are captured in the error definition. On the other hand, the Eq. (23) still excludes from $\{c\}$ the cases with both $u_n$ and $u'_n$ being negative, which do not contribute anyway to force transmission and are mechanically irrelevant. The data shown in the Fig. 19 then confirm how an increasing number of surface nodes leads to a greater consistency for overlap description in the LS algorithm, where $\Delta^*(u_n)$ reaches negligible values of 5% or below when using several thousands of surface nodes (3602 or 4902) as in previous section 4.2.

Regarding contact normal, the default estimation of the latter as $\vec{n}$, in Eqs. (1)-(2), is compared with the following Eq. (25) where master and slave bodies are exchanged to output $\vec{n}'$, forming another error metric $\delta(\vec{n})$ which is locally defined in Eq. (26) as the angle between $\vec{n}$ and $\vec{n}'$ for every interaction where actual contact is detected in either workflow, with both $u_n$ and $u'_n$ being positive:

$$\vec{n}' = -\vec{V}\phi_S \text{ at the } \vec{N} \in B \text{ corresponding to Eq. (22)} \qquad (25)$$

$$\delta(\vec{n}) = \arcsin(\|\vec{n} \wedge \vec{n}'\|) \in [0^{\circ}; 90^{\circ}] \qquad (26)$$

<figure><img src="image_27.png" /><figcaption>Fig. 19. Vanishing asymmetry of the LS contact algorithm regarding interparticle overlap, in terms of $\Delta^*(u_n)$, Eq. (23), as a function of used number of surface nodes.</figcaption></figure>

The distribution of $\delta(\vec{n})$ is analyzed in the Fig. 20 for the two packings presently at hand. While extreme values of several dozens of degrees for $\delta(\vec{n})$ are systematically possible for a small number of interactions, the data (Figs. 20a and 20b) confirm that an increasing number of nodes leads to a clear decreasing trend in normal orientation mismatch between $\vec{n}$ and $\vec{n}'$ for a vast majority (more than 75%) of particle pairs. For that majority of contacts, a negligible mismatch of few degrees at most is actually observed between $\vec{n}$ and $\vec{n}'$ when using several thousands of surface nodes as in previous section 4.2. Said otherwise, the proportion of contacts with a significant error on contact normal, i.e. greater than $15^{\circ}$ as an arbitrary threshold, clearly decreases with the number of surface nodes and reaches negligible values below 3% for the same choices of several thousands of surface nodes (Fig. 20c), confirming how the LS approach closely approximates the common contact normal conditions.

## 5. Conclusion

Restricting to the case of strictly convex particles, the Level Set-Discrete Element Method has been optimized and further discussed as a physically adequate modeling approach for discrete mechanics simulations, with no mechanically adverse effects of its master-slave contact algorithm.

Optimization has been proposed for what concerns time costs in several aspects, including the initial generation of the LS distance field and a double modification affecting the surface nodes discretization in terms of its definition and DEM treatment. Regarding the surface nodes definition, an original placement method has been applied for locating in a more uniform manner the nodes over the particle surface, where a direct LS-based estimation of the latter has also been proposed and validated. Regarding the surface nodes-based DEM contact treatment, the proposed consideration of contact history has enabled to dramatically reduce time costs by, e.g., 67% for quasi-static simulations.
<figure><img src="image_28.png" /><figcaption>(a) $\delta(\vec{n})$ distribution for Shape A (minima are at $10^{-6}\circ$ or below, out of graph scale)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(b) $\delta(\vec{n})$ distribution for Shape E (minima are around $10^{-4}\circ$ and out of scale)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(c) Proportion of contacts with $\delta(\vec{n}) \ge 15^\circ$</figcaption></figure>

<caption>**Fig. 20.** Level Set approximated description of a common contact normal for isotropically stressed ($\rho = 100$ kPa) packings of 8000 particles of either Shape A or E, in terms of statistical properties of the $\delta(\vec{n})$ distribution and of the proportion of contacts with $\delta(\vec{n}) \ge 15^\circ$.</caption>

Such improvements have led the consideration of several thousands of surface nodes per particle to be actually affordable, at the same time it has also been shown to be necessary for a sound, symmetric, LS contact description that closely matches the unique definitions of contact overlap and contact normal which are physically expected.

As an example of a possible application in future dedicated studies, the optimized LS approach has been applied to a preliminary analysis of the residual shear strength of packings made of convex superquadric shapes, whose definitions in terms of a Representative Elementary Volume have been justified. While the shear strength has been shown to lack a unique relationship with Wadell's true sphericity serving as a particle shape index, it has been quantified to be up to ~73% higher than the one of an otherwise similar packing made of spheres, illustrating once more the great influence of particle shape in discrete mechanics.

### CRediT authorship contribution statement

**Jérôme Duriez:** Writing – original draft, Software, Methodology, Investigation, Conceptualization.
**Cédric Galusinski:** Writing – review & editing, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

Other YADE developers, passed and present (within the yade-dev GitLab team), are acknowledged for setting up and maintaining the

<figure><img src="image_31.png" /><figcaption>**Fig. A.21.** Reduced time costs thanks to the use of a heap-sort procedure when executing a Fast Marching Method (FMM) on grids with various resolutions, the latter being expressed in terms of spacing $g_{grid}$ with respect to the particle size $D$, or in terms of the total number of gridpoints, $N_{gp}$. Each case is repeated twice.</figcaption></figure>

platform grounding the proposed implementation. The first author also acknowledges Dr Christine Valentin and Mrs Michèle Duriez for the fruitful working conditions he enjoyed during most of the manuscript construction.

## Appendix A. Computational analysis of the FMM computation of distance field

The effect of using a new heap-sort procedure during a FMM definition of distance fields in the §2.3 is illustrated in Fig. A.21 when

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132
<caption>Table B.5</caption>

Surface estimations (in square length unit, up to 8 digits in the fractional part) for two ellipsoids, including one spheroid.

<table><thead><tr><th></th><th>Exact</th><th>Present numeric<br>Python-based estimation</th><th>Ref [59]</th><th>Eq. (B.1) [61]</th></tr></thead><tbody><tr><td>Oblate spheroid<br>r<sub>x</sub> = r<sub>y</sub> = 1 &gt; r<sub>z</sub> = 0.5<br>(in length unit)</td><td>8.67188270</td><td>8.67188270</td><td>8.67188271</td><td>8.67263054</td></tr><tr><td>Ellipsoid<br>r<sub>x</sub> = 2; r<sub>y</sub> = 1; r<sub>z</sub> = 0.25<br>(in length unit)</td><td>NA</td><td>&isin; [13.69921079; 13.69921083]</td><td>13.69921081</td><td>13.71084885</td></tr></tbody></table>

computing, for illustration purpose, the distance field to a sphere of diameter *D* on a laptop with a 0.8-5.1 GHz Intel Xeon W-10855M CPU manually fixed at 4.0 GHz and a 3.2 GHz 32 GB RAM. Distance computations now stay around the second even for fine grids with hundred of voxels per diameter, which corresponds in those cases to a speed-up greater than 100 with respect to the previous implementation of [24]. On the other hand, FMM errors with reference to the theoretical distance field to a sphere are strictly unchanged with respect to [24] where the first-order nature of the numerical scheme was detailed.

## Appendix B. Validation of the numeric estimation of surface area for ellipsoids

While not including the Level Set-based approach discussed in the main text, §2.3, the Table B.5 compares different methods for estimating the surface of two ellipsoids, with $\epsilon_e = \epsilon_n = 1$ in Eq. (7), one of them being furthermore a spheroid with $r_x = r_y$. The comparison first includes, for the latter case of the spheroid only, the exact, closed-form, expression available in this case and recalled, e.g. in [59]. An original Python-based numeric estimation is also considered, whose implementation is available in the attached function `superquadricSurface.py`. It symbolically derives (with `sympy` module in its 1.9 version) the infinitesimal surface of a superquadric after differentiating Eqs. (3)-(5), similar to [60], whose integral is eventually numerically evaluated to the full surface through numerical quadrature (with `scipy` 1.8.0). That herein-proposed procedure is compared with the values proposed in [59] from alternate quadrature considerations presented in the form of a convergent series (truncated to its 40 first terms in [59]) and to the following closed-form approximation of Eq. (B.1) proposed in [61] (see also [62,63]) and used, e.g. in [64]:

$$S = 2\pi \left[ (r_x r_y)^p + (r_y r_z)^p + (r_x r_z)^p \right]^{1/p} \quad (B.1)$$

$$p = \log_2(3) \quad (B.2)$$

The results shown in Table B.5 validate the present Python-based surface estimation since the latter correctly estimates the true surface area of a spheroid up to the $10^{-8}$ digit, which means a relative error being strictly smaller than $5.8 \times 10^{-10}$. A just slightly greater relative error, around $1.2 \times 10^{-9}$, was obtained in [59] whereas it would be significantly more important, around $8 \times 10^{-5}$ and with discrepancies from the $10^{-3}$ digit, for the Eq. (B.1) taken from [61]. Similar observations hold for the other considered case of a true ellipsoid, with no known exact surface value, where the herein proposed estimation and the results of [59] approach one another within a relative variation of $1.5 \times 10^{-9}$ i.e. $1.5 \times 10^{-7}\%$, while Eq. (B.1) outputs a value which is $8.5 \times 10^{-2}\%$ different from that common estimate.

## Appendix C. Precision of the level set description of surface features for superquadric shapes B to E

The Fig. C.22 adds to the previous Fig. 4 for analyzing how the grid resolution affects the precision of the LS approach for describing superquadrics, be it in terms of surface points positions or surface

normal orientation. Confirming the observations made on Fig. 4 for shape A, it again appears that using a grid resolution, as defined by $2\min(r_x, r_y, r_z)/g_{grid}$ of at least 20 enables one to keep median errors on normal orientation below the degree and (dimensionless) errors on surface positions well below $10^{-2}$ for other shapes B to E.

## Appendix D. Quantifying the uniformity of nodes distribution

When interpreting in the § 3.2 the more optimal nature of the spiral path (or of the refined cosahedron) for distributing nodes over a sphere, the $\{\min_j(d_{ij})\}$ distributions of shortest internode distances have been considered in the Fig. 7 where the spiral case has shown to be more uniform, i.e. with a much smaller distribution width. The $\min_j(d_{ij})$ values in the distribution may nevertheless appear to be slightly greater for the spiral path versus the rectangular one (see also the Fig. D.23a), which is inconsistent with the idea of a better LS-DEM performance and this observation is herein interpreted as corresponding to a small bias of that metric.

Another metric could actually be considered after building for each node *i* a quadrilateral that includes *i* and three of its neighboring nodes and which is seen as corresponding to a hole in the discretization. A hole characteristic length $l_i$ is subsequently defined as the ratio between the quadrilateral surface area and its perimeter. Even though each set of four nodes is not necessarily coplanar, the surface area of each quadrilateral is computed by summing the surfaces of the two triangles subdividing the latter. For the same spiral and rectangular distributions of 1602 nodes on a sphere of diameter *D*, the Fig. D.23b depicts the distribution of $l_i/D$. Here, the interval of values for $l_i/D$ with the spiral path is strictly included into a wider one for the rectangular path, suggesting that second metric has a slightly better definition.

It has still been observed that results are virtually the same if it is attempted to minimize the distribution widths of either $\min_j(d_{ij})/D$ or $l_i/D$ (or the ratio between the standard deviation of $l_i/D$ divided by its mean), which are all considered to be appropriate to explain the greater LS-DEM performance of the spiral distribution with spherical particles (Fig. 6).

## Appendix E. Supplementary material

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.cpc.2025.109803.

## Data availability

The DEM code is available at https://gitlab.com/yade-dev/trunk. Proposed surface discretization nodes as well as a Python function for computing superquadric surfaces are also directly included.

## References

[1] P. Cundall, O. Strack, A discrete numerical model for granular assemblies, Géotechnique 29 (1979) 47–65, https://doi.org/10.1680/geot.1979.29.1.47.

[2] G. Combe, J.-N. Roux, Strain versus stress in a model granular material: a devil's staircase, Phys. Rev. Lett. 85 (2000) 3628–3631, https://doi.org/10.1103/PhysRevLett.85.3628.

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132
<figure><img src="image_32.png" /><figcaption>(a) Shape B: surface position</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(b) Shape B: normal vector</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(c) Shape C: surface position</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(d) Shape C: normal vector</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(e) Shape D: surface position</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(f) Shape D: normal vector</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(g) Shape E: surface position</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(h) Shape E: normal vector</figcaption></figure>

<caption>Fig. C.22. Statistical properties of the errors in surface position or surface normal of the Level Set description of Shapes B to E.</caption>
<figure><img src="image_40.png" /><figcaption>(a) Smallest distance between a surface node and all other nodes</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(b) Characteristic length of a hole defined for every node</figcaption></figure>

<caption>**Fig. D.23.** Distributions on the two considering geometric quantities for analyzing the surface coverage on a sphere when using either the rectangular or the spiral set of nodes locations (using 1602 nodes).</caption>

[3] G.D.R. MiDi, On dense granular flows, Eur. Phys. J. E 14 (2004) 341–365, https://doi.org/10.1140/epje/i2003-10153-0.

[4] R. Kawamoto, E. Andò, G. Viggiani, J.E. Andrade, Level set discrete element method for three-dimensional computations with triaxial case study, J. Mech. Phys. Solids 91 (2016) 1–13, https://doi.org/10.1016/j.jmps.2016.02.021.

[5] L. Sibille, P. Villard, F. Darve, R. Aboul Hosn, Quantitative prediction of discrete element models on complex loading paths, Int. J. Numer. Anal. Methods Geomech. 43 (5) (2019) 858–887, https://doi.org/10.1002/nag.2911.

[6] K. Karapiperis, J. Harmon, E. Andò, G. Viggiani, J.E. Andrade, Investigating the incremental behavior of granular materials with the level-set discrete element method, J. Mech. Phys. Solids 144 (2020) 104103, https://doi.org/10.1016/j.jmps.2020.104103.

[7] T. Mohamed, J. Duriez, G. Veylon, L. Peyras, DEM models using direct and indirect shape descriptions for Toyoura sand along monotonous loading paths, Comput. Geotech. 142 (2022) 104551, https://doi.org/10.1016/j.compgeo.2021.104551.

[8] T. Mohamed, J. Duriez, G. Veylon, L. Peyras, P. Soulat, A 3D-DEM model for tropical residual soils under monotonic and cyclic loadings, J. Geotech. Geoenviron. Eng. 149 (11) (2023) 1–13, https://doi.org/10.1061/JGEFK.GTENG-11323.

[9] S. Duverger, V. Angelidakis, S. Nadimi, S. Utili, S. Bonelli, P. Philippe, J. Duriez, Investigation techniques and physical aspects of the angle of repose of granular matter, Granul. Matter 26 (20) (2024), https://doi.org/10.1007/s10035-023-01378-z.

[10] T. Mede, G. Chambon, P. Hagenmuller, F. Nicot, A medial axis based method for irregular grain shape representation in DEM simulations, Granul. Matter 20 (1) (2018) 16, https://doi.org/10.1007/s10035-017-0785-7.

[11] V. Angelidakis, S. Nadimi, M. Otsubo, S. Utili, CLUMP: a code library to generate universal multi-sphere particles, SoftwareX 15 (2021) 100735, https://doi.org/10.1016/j.softx.2021.100735.

[12] H. Fathipour-Azar, J. Duriez, A comparative study of existing and new sphere clump generation algorithms for modeling arbitrary shaped particles, Arch. Comput. Methods Eng. (2025), https://doi.org/10.1007/s11831-025-10256-1.

[13] F. Dubois, Numerical modeling of granular media composed of polyhedral particles, in: F. Radjai, F. Dubois (Eds.), Discrete-Element Modeling of Granular Materials, ISTE-Wiley, 2011, pp. 233–262.

[14] A. Wachs, L. Girolami, G. Vinay, G. Ferrer, Grains3d, a flexible DEM approach for particles of arbitrary convex shape – part I: numerical model and validations, Powder Technol. 224 (2012) 374–389, https://doi.org/10.1016/j.powtec.2012.03.023.

[15] A.D. Rakotonirina, A. Wachs, Grains3d, a flexible DEM approach for particles of arbitrary convex shape - part II: parallel implementation and scalable performance, Powder Technol. 324 (2018) 18–35, https://doi.org/10.1016/j.powtec.2017.10.033.

[16] Y. Descantes, A new discrete element modelling approach to simulate the behaviour of dense assemblies of true polyhedra, Powder Technol. 401 (2022) 117295, https://doi.org/10.1016/j.powtec.2022.117295.

[17] E. Gilbert, D. Johnson, S. Keerthi, A fast procedure for computing the distance between complex objects in three-dimensional space, IEEE J. Robot. Autom. 4 (2) (1988) 193–203, https://doi.org/10.1109/56.2083.

[18] S. Zhao, J. Zhao, A poly-superellipsoid-based approach on particle morphology for DEM modeling of granular media, Int. J. Numer. Anal. Methods Geomech. 43 (13) (2019) 2147–2169, https://doi.org/10.1002/nag.2951.

[19] G. Lu, J. Third, C. Müller, Critical assessment of two approaches for evaluating contacts between super-quadric shaped particles in DEM simulations, Chem. Eng. Sci. 78 (2012) 226–235, https://doi.org/10.1016/j.ces.2012.05.041.

[20] A. Podlozhynyuk, S. Pirker, C. Kloss, Efficient implementation of superquadric particles in discrete element method within an open-source framework, Comput. Part. Mech. 4 (2017) 101–118, https://doi.org/10.1007/s40571-016-0131-6.

[21] J. Zhao, S. Zhao, S. Luding, The role of particle shape in computational modelling of granular matter, Nat. Rev. Phys. 5 (9) (2023) 505–525, https://doi.org/10.1038/s42254-023-00617-9.

[22] A.X. Jerves, R.Y. Kawamoto, J.E. Andrade, Effects of grain morphology on critical state: a computational analysis, Acta Geotech. 11 (3) (2016) 493–503, https://doi.org/10.1007/s11440-015-0422-8.

[23] J. Duriez, S. Bonelli, Precision and computational costs of level set-discrete element method (LS-DEM) with respect to DEM, Comput. Geotech. 134 (2021) 104033, https://doi.org/10.1016/j.compgeo.2021.104033.

[24] J. Duriez, C. Galusinski, A level set-discrete element method in YADE for numerical, micro-scale, geomechanics with refined grain shapes, Comput. Geosci. 157 (2021) 104936, https://doi.org/10.1016/j.cageo.2021.104936.

[25] G. Houlsby, Potential particles: a method for modelling non-circular particles in DEM, Comput. Geotech. 36 (6) (2009) 953–959, https://doi.org/10.1016/j.compgeo.2009.03.001.

[26] Z. Lai, J. Zhao, S. Zhao, L. Huang, Signed distance field enhanced fully resolved CFD-DEM for simulation of granular flows involving multiphase fluids and irregularly shaped particles, Comput. Methods Appl. Mech. Eng. 414 (2023) 116195, https://doi.org/10.1016/j.cma.2023.116195.

[27] J. Sethian, A fast marching level set method for monotonically advancing fronts, Proc. Natl. Acad. Sci. 93 (4) (1996) 1591–1595, https://doi.org/10.1073/pnas.93.4.1591.

[28] L. Li, E. Marteau, J.E. Andrade, Capturing the inter-particle force distribution in granular material using LS-DEM, Granul. Matter 21 (3) (2019) 43, https://doi.org/10.1007/s10035-019-0893-7.

[29] S. Feldfogel, K. Karapiperis, J. Andrade, D.S. Kammer, A discretization-convergent level-set-discrete-element-method using a continuum-based contact formulation, Int. J. Numer. Methods Eng. 125 (5) (2024) e7400, https://doi.org/10.1002/nme.7400.

[30] A.D. Davis, B.A. West, N.J. Frisch, D.T. O'Connor, M.D. Parno, ParticLS: object-oriented software for discrete element methods and peridynamics, Comput. Part. Mech. 9 (1) (2022) 1–13, https://doi.org/10.1007/s40571-021-00392-3.

[31] D.L. van der Haven, I.S. Fragkopoulos, J.A. Elliott, A physically consistent discrete element method for arbitrary shapes using volume-interacting level sets, Comput. Methods Appl. Mech. Eng. 414 (2023) 116165, https://doi.org/10.1016/j.cma.2023.116165.

[32] V. Smilauer, et al., Yade Documentation, 3rd ed., The Yade Project, 2021, http://yade-dem.org/doc/.

[33] V. Angelidakis, K. Boschi, K. Brzeziński, R.A. Caulk, B. Chareyre, C.A. del Valle, J. Duriez, A. Gladky, D.L. van der Haven, J. Kozicki, G. Pekmezi, L. Scholtès, K. Thoeni, YADE - an extensible framework for the interactive simulation of multiscale, multiphase, and multiphysics particulate systems, Comput. Phys. Commun. (2024) 109293, https://doi.org/10.1016/j.cpc.2024.109293.

[34] A.H. Barr, Superquadrics and angle-preserving transformations, IEEE Comput. Graph. Appl. 1 (1) (1981) 11–23, https://authors.library.caltech.edu/9756/1/BARieeeeca81.pdf.

[35] A.H. Barr, Rigid physically based superquadrics, in: D. Kirk (Ed.), Graphics Gems III, Academic Press, 1995, pp. 137–159.

[36] H. Wadell, Volume, shape, and roundness of rock particles, J. Geol. 40 (5) (1932) 443–451, https://doi.org/10.1086/623964.

[37] R. Rorato, M. Arroyo, E. Andò, A. Gens, Sphericity measures of sand grains, Eng. Geol. 254 (2019) 43–53, https://doi.org/10.1016/j.enggeo.2019.04.006.

[38] G. Mollon, J. Zhao, Generating realistic 3D sand particles using Fourier descriptors, Granul. Matter 15 (2013) 95–108, https://doi.org/10.1007/s10035-012-0380-x.

[39] E.A. Rakhmanov, E.B. Saff, Y.M. Zhou, Minimal discrete energy on the sphere, Math. Res. Lett. 1 (1994) 647–662, https://doi.org/10.4310/MRL.1994.v1.n6.a3.

[40] J.-N. Roux, G. Combe, Quasistatic rheology and the origins of strain, C. R. Phys. 3 (2) (2002) 131–140, https://doi.org/10.1016/S1631-0705(02)01306-3.

[41] B. Fornberg, N. Flyer, Fast generation of 2-D node distributions for mesh-free PDE discretizations, Comput. Math. Appl. 69 (7) (2015) 531–544, https://doi.org/10.1016/j.camwa.2015.01.009.

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132
[42] D.P. Hardin, E.B. Saff, Discretizing manifolds via minimum energy points, Not. Am. Math. Soc. 51 (10) (2004) 1186-1194.

[43] P.A. Cundall, Distinct element models of rock and soil structure, in: E.T. Brown (Ed.), Analytical and Computational Methods in Engineering Rock Mechanics, George Allen and Unwin, 1987, pp. 129-163.

[44] J. Thomson XXIV, On the structure of the atom: an investigation of the stability and periods of oscillation of a number of corpuscles arranged at equal intervals around the circumference of a circle with application of the results to the theory of atomic structure, Lond. Edinb. Dublin Philos. Mag. J. Sci. 7 (39) (1904) 237-265, https:// doi.org/10.1080/14786440409463107.

[45] P.J. Frey, H. Borouchaki, Geometric surface mesh optimization, Comput. Vis. Sci. 1 (3) (1998), https://doi.org/10.1007/s007910050011.

[46] P. Alliez, E. de Verdire, O. Devillers, M. Isenburg, Isotropic surface remeshing, in: Proceedings of the Shape Modeling International 2003, Seoul, Korea (South), 2003, pp. 49-58, https://graphics.stanford.edu/courses/cs164-10-spring/ Handouts/isotropic.pdf.

[47] P.-O. Persson, G. Strang, A simple mesh generator in Matlab, SIAM Rev. 46 (2) (2004) 329-345, https://doi.org/10.1137/S0036144503429121.

[48] R. Zamolo, E. Nobile, Two algorithms for fast 2D node generation: application to RBF meshless discretization of diffusion problems and image halftoning, Comput. Math. Appl. 75 (12) (2018) 4305-4321, https://doi.org/10.1016/j.camwa.2018.03.031.

[49] M.V. Craveiro, A.G. Neto, P. Wriggers, Contact between rigid convex NURBS particles based on computer graphics concepts, Comput. Methods Appl. Mech. Eng. 386 (2021) 114097, https://doi.org/10.1016/j.cma.2021.114097.

[50] M.V. Craveiro, A. Gay Neto, P. Wriggers, DEM simulations using convex NURBS particles, Comput. Part. Mech. 11 (3) (2024) 1087-1118, https://doi.org/10.1007/ s40571-023-00675-x.

[51] S. Zhao, J. Zhao, SudoDEM: unleashing the predictive power of the discrete element method on simulation for non-spherical granular particles, Comput. Phys. Commun. 259 (2021) 107670, https://doi.org/10.1016/j.cpc.2020.107670.

[52] H. Ouadfel, L. Rothenburg, 'Stress-force-fabric' relationship for assemblies of ellipsoids, Mech. Mater. 33 (4) (2001) 201-221, https://doi.org/10.1016/S0167- 6636(00)00057-0.

[53] C. Thornton, Granular Dynamics, Contact Mechanics and Particle System Simulations, Springer, 2015.

[54] K.H. Roscoe, A.N. Schofield, C.P. Wroth, On the yielding of soils, Géotechnique 8 (31) (1958) 22-53, https://doi.org/10.1680/geot.1958.8.1.22.

[55] S. Zhao, T.M. Evans, X. Zhou, Effects of curvature-related DEM contact model on the macro- and micro-mechanical behaviours of granular soils, Géotechnique 68 (12) (2018) 1085-1098, https://doi.org/10.1680/jgeot.17.P.158.

[56] C. Wellmann, C. Lillie, P. Wriggers, Homogenization of granular material modeled by a three-dimensional discrete element method, Comput. Geotech. 35 (3) (2008) 394-405, https://doi.org/10.1016/j.compgeo.2007.06.010.

[57] X. Lin, T.-T. Ng, Contact detection algorithms for three-dimensional ellipsoids in discrete element modelling, Int. J. Numer. Anal. Methods Geomech. 19 (9) (1995) 653-659, https://doi.org/10.1002/nag.1610190905.

[58] F. Radjai, D.E. Wolf, M. Jean, J.-J. Moreau, Bimodal character of stress transmission in granular packings, Phys. Rev. Lett. 80 (1) (1998) 61-64.

[59] S.R. Keller, On the surface area of the ellipsoid, Math. Comput. 33 (145) (1979) 310-314.

[60] M. Evans, Some properties of superellipsoid mixtures, in: 5th International Symposium on Geomechanics from Micro to Macro, IS-Grenoble 2024, Grenoble, 2024.

[61] K. Thomsen, Surface area of an ellipsoid, in: G.P. Michon (Ed.), Final Answers to Scientific Questions: Spheroids & Scalene Ellipsoids, 2004, https://www.numericana.com/answer/ellipsoid.htm#thomsen.

[62] M.S. Klamkin, Elementary approximations to the area of n-dimensional ellipsoids, Am. Math. Mon. 78 (3) (1971) 280-283, https://doi.org/10.2307/2317530, http:// www.jstor.org/stable/2317530.

[63] M.S. Klamkin, Corrections to “elementary approximations to the area of n-dimensional ellipsoids”, Am. Math. Mon. 83 (6) (1976) 478, https://doi.org/10. 2307/2318347, http://www.jstor.org/stable/2318347.

[64] M. Taylor, E. Garboczi, S. Erdogan, D. Fowler, Some properties of irregular 3-d particles, Powder Technol. 162 (1) (2006) 1-15, https://doi.org/10.1016/j.powtec.2005. 10.013.
<figure><img src="image_42.png" /></figure>

ELSEVIER

Contents lists available at ScienceDirect

Computer Physics Communications

journal homepage: www.elsevier.com/locate/cpc

<figure><img src="image_43.png" /></figure>

Graphical abstract

Faster and objective level set-DEM mechanical simulations of discrete systems with convex particles from contact history and particle surface considerations

Computer Physics Communications ••••, •••, •••

Jérôme Duriez<sup>a,∗</sup>, Cédric Galusinski<sup>b</sup>

<figure><img src="image_44.png" /></figure>

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

<footnote><sup>a</sup> INRAE, Aix Marseille Univ, RECOVER, Aix-en-Provence, France</footnote>

<footnote><sup>b</sup> IMATH, Université de Toulon, CS 60584, 83041 Toulon Cedex 9, France</footnote>