A Quick Guide to SudoDEM : A Discrete

Element Code for Non-spherical Particles

Shiwei Zhao1,2 and Jidong Zhao1

1 The Hong Kong University of Science and Technology
2 South China University of Technology

Updated on June 13, 2020, version 1.2.1

 A Quick Guide to SudoDEM

2
Contents

1 Introduction 5
1.1 Why SudoDEM . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Disclaimers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 Requirements 7
2.1 Basic Linux Commands . . . . . . . . . . . . . . . . . . . . . 7
2.2 Python 2.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Package Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3.1 Binary Installation . . . . . . . . . . . . . . . . . . . . 9
2.3.2 Non-root Compilation . . . . . . . . . . . . . . . . . . 10

3 Examples 17
3.1 Example 0: run a simulation . . . . . . . . . . . . . . . . . . . 17
3.2 Example 1: packing of super-ellipsoids . . . . . . . . . . . . . 20
3.3 Example 2: triaxial tests of super-ellipsoids . . . . . . . . . . . 28
3.4 Example 3: packing of GJKparticles . . . . . . . . . . . . . . 38
3.5 Example 4: packing of poly-superellipsoids . . . . . . . . . . . 46
3.6 Example 5: packing of superellipses . . . . . . . . . . . . . . . 50

4 Post-processing 55
4.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.2 Scene Visualization . . . . . . . . . . . . . . . . . . . . . . . . 56
4.2.1 SudoDEM3D . . . . . . . . . . . . . . . . . . . . . . . 56
4.2.2 SudoDEM2D . . . . . . . . . . . . . . . . . . . . . . . 60

5 Python Class Reference 63

5.1 SudoDEM3D . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.1.1 Basic classes . . . . . . . . . . . . . . . . . . . . . . . . 63
5.1.2 Module _superquadrics_utils . . . . . . . . . . . . . . 64
5.1.3 Module _gjkparticle_utils . . . . . . . . . . . . . . . . 66

3
CONTENTS A Quick Guide to SudoDEM

5.1.4 Module snapshot . . . . . . . . . . . . . . . . . . . . . 68

5.2 SudoDEM2D . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.2.1 Basic classes . . . . . . . . . . . . . . . . . . . . . . . . 69
5.2.2 Module _superellipse_utils . . . . . . . . . . . . . . . 70
5.2.3 Module _utils . . . . . . . . . . . . . . . . . . . . . . . 72
5.2.4 Module utils . . . . . . . . . . . . . . . . . . . . . . . . 73
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4
Chapter 1

Introduction

1.1 Why SudoDEM

A sound name is the beginning to get things started. The word SudoDEM is
coined as a combination of ‘sudo’ and ‘DEM’. The prefix ‘sudo’ is a program
for Unix-like computer operating systems that allows users to run programs
as the superuser, i.e., ‘super user do’, which here means a powerful and flex-
ible DEM simulator. On the other hand, the pronunciation of ‘sudo’ sounds
like ‘pseudo’, implying that SudoDEM does and will consist of features dis-
tinguishing from the conventional DEM codes.
SudoDEM is specifically designed for modeling non-spherical particles us-
ing discrete element method (DEM), which inherits a basic framework of an
open-source DEM code, YADE1 . A profound modification was performed on
YADE for higher efficiency in modeling non-spherical particles. The project
is hosted on ResearchGate with a goal of developing a robust 2/3D DEM
code for convex particles, e.g., super-ellipsoids, poly-superellipsoids, cylin-
ders, cones, polyhedrons. In SudoDEM, some general optimization algo-
rithms e.g., Levenberg-Marquardt and Nelder-Mead simplex are adopted for
contact detection of superquadric particles. The popular Gilbert-Johnson-
Keerthi (GJK) algorithm used in computer graphic simulation is employed
for convex and non-convex (under construction) polytopes.
The project of SudoDEM is hosted on an individual website, and the
synchronous update can be also found on the Researchgate Page.

1
[Link]

5
CHAPTER 1. INTRODUCTION A Quick Guide to SudoDEM

Figure 1.1: GUIs of SudoDEM 3D (left) and SudoDEM 2D (right).

1.2 Related Work

The users are appreciated to cite our work as, but not limited to, listed below.

(1) Zhao S., Zhao J. (2019). SudoDEM : an open-source discrete element

code for non-spherical particles, in preparation.

(2) Zhao S., Zhao J. (2019). A poly-superellipsoid-based approach on

particle morphology for DEM modeling of granular media. Interna-
tional Journal for Numerical and Analytical Methods in Geomechanics,
43(13): 2147–2169.

(3) Zhao S., Evans T. M., Zhou X. (2018). Effects of curvature-related

DEM contact model on the macro-and micro-mechanical behaviours of
granular soils. Gétechnique, 68(12): 1085–1098.

(4) Zhao S., Zhang N., Zhou X., Zhang L. (2017). Particle shape effects on
fabric of granular random packing. Powder technology, 310, 175–186.

1.3 Disclaimers
BECAUSE THE CODE IS FREE OF CHARGE, THERE IS NO WAR-
RANTY ‘AS IS’; NOT EVEN FOR MERCHANTABILITY OR FITNESS
FOR A PARTICULAR PURPOSE.

6
Chapter 2

Requirements

A Linux distribution (Ubuntu 14 to 18 suggested) with Python 2.7 and some

basic python packages (e.g., numpy) are assumed to have been installed at the
users’ side, but no other third-party packages (e.g., boost, qt) need installing.

2.1 Basic Linux Commands

Only several commands the users should know are enough to set off running a
program. A terminal (i.e., command window) is open by pressing a shortcut
‘CTRL+ALT+T’ or whatever other operations.

• sudo apt-get install package 1 install a package to the system with root
permission.

• pwd: show the current path.

• cd directory: change the directory to directory. The directory can be

a relative or absolute path. Two special notations ‘.’ and ‘..’ denote
the current directory and the parent directory, respectively, which are
useful to type a relative path. For example, ‘ cd ../sub’ will change
the directory to the subdirectory ‘sub’ of the parent directory of the
current path. ‘cd /home/user’ will change the directory to the absolute
path ‘/home/user’.

• prog or ./path/prog: run an executable named prog. If prog is ex-

posed to the terminal, i.e., prog is in the searching path of the system,
then just type the name prog like the command ‘cd’ at the last item;
otherwise, a path should be specified to locate prog.
1
or using sudo apt install package

7
 CHAPTER 2. REQUIREMENTS A Quick Guide to SudoDEM

• ‘CTRL+C’: terminate a program at the terminal.

1 $> sudo apt−get install python−numpy

2 $> pwd
3 $> cd /home/
4 $> sudodem3d

2.2 Python 2.7

Python is sensitive to indent which is used to identify a script block. ‘#’ is
used to comment a line which will not be executed. Here a piece of script is
shown to explain the primitives of Python:
1 import math # import t h e module math which i n c l u d e s some b a s i c
math f u n c t i o n s
2
3 # we have numbers
4 a = 2 # i n t e g e r number
5 b = 5
6 c = 2 . 0 # f l o a t number
7 d = a /b # t h e v a l u e i s an i n t e g e r
8 e = c /b # f l o a t number
9 p r i n t a , d , e # p r i n t t h e v a l u e s o f a , d and e t o t h e t e r m i n a l
10
11 #we have s t r i n g s
12 a = ' abc ' # we can a s s i g n any type o f v a l u e s t o any type o f
variables
13 b = " t h i s i s a ' Python ' s c r i p t . "
14 print a , b
15
16 #we have l i s t , t u p l e and d i c t t o s t o r e and manage t h e above b a s i c
primitives
17 c = [1 ,2 ,3 , ' ss ' ] # a l i s t
18 d = (1 ,2 ,3 , ' ss ' ) # a tuple
19 c [ 0 ] = 5 # change t h e f i r s t item i n t h e l i s t c , and t h e i n d e x
s t a r t s from 0 .
20 d [ 0 ] = 5 # f a i l e d ! a t u p l e i s c o n s t a n t , which i s t h e
d i s t i n g u i s h i n g p r o p e r t y from a l i s t
21 d = { 1 : 2 2 , 2 : 3 3 , 4 : ' s s s s ' , 'w ' : 1 2 3 } #a d i c t : 1 , 2 , 4 and 'w ' a r e k e y s
o f t h e d i c t , and 2 2 , 3 3 , ' s s s s ' and 123 a r e t h e c o r r e s p o n d i n g
values .
22 p r i n t d [ 1 ] , d [ 4 ] , d [ 'w ' ]
23

24 # we i n i t i a l i z e two o b j e c t s c and d f o r s t o r i n g data a t a f o r

loop l a t e r
25 c = l i s t ( ) # an empty l i s t

8
 A Quick Guide to SudoDEM CHAPTER 2. REQUIREMENTS

26 d = d i c t ( ) # an empty d i c t
27 #t e s t e x p r e s s i o n
28 i f 1 >2:
29 p r i n t " 1 i s g r e a t e r than 2 , r e a l l y ? "
30 else :
31 i f 2 in [ 1 , 2 , 3 ] : # i f 2 i s in the l i s t [ 1 , 2 , 3 ]
32 #ok , l e t ' s s t a r t a f o r l o o p
33 f o r i in range (5) : # equal to f o r i in [ 0 , 1 , 2 , 3 , 4 ]
34 p r i n t i # you can s e e what ' s g o i n g on i n t h e f o r l o o p
35 # we append some data t o a l i s t c and a d i c t d
36 c . append ( math . c o s ( i ) ∗10 + 2 . 0 ∗ ∗ 5 ) # we append t h e v a l u e o f
10∗ c o s ( i ) + 2 . 0 ^ 5 t o t h e l i s t c . We u s e t h e math f u n c t i o n c o s
from t h e module math .
37 i f i not i n d . k e y s ( ) : # i f i i s not a key o r i n d e x o f t h e
dict d.
38 d[ i ] = str ( i )
39
40 #we can c a p s u l e our commands by a f u n c t i o n
41 d e f GoodJob ( input_number , keyword1 = 1 , keyword2 = True ) :
42 i f keyword2 : # t h a t i s i f keyword2 == True
43 p r i n t " t h e i n p u t number i s " , i n p u t
44 input_number += keyword1 # i . e . , t h e v a l u e o f input_number i s
updated t o input_number + keyword1
45 r e t u r n input_number # we can r e t u r n t h e v a l u e
46
47 # c a l l the f u n c t i o n
48 a = GoodJob ( 2 ) # a = 3 , with p r i n t i n f o ' t h e i n p u t number i s 3 '
49 b = GoodJob ( 3 , keyword1 = 5 , keyword2 = F a l s e ) # b = 8 , w i t h o u t
print info .

Ok, that is all you need to know prior to running a simulation. There are also
other further techniques you might be interested in, e.g., how to open/read-
/write/close a file, and how to write a Python script with the object-oriented
programming, which will be roughly introduced in the examples of this tu-
torial.

2.3 Package Setup

The project of SudoDEM is hosted on the website. Two ways to get the
package installed are provided:

2.3.1 Binary Installation

Step 1: Get the package2
2
The binary is only compiled for Ubuntu 14 to 18.

9
CHAPTER 2. REQUIREMENTS A Quick Guide to SudoDEM

Download the compiled package from the download page and unzip it to
anywhere (e.g., /home/xxx/). If you just get a package of the main program
which does not include the third-party libraries, then you need to download
the package "[Link]" for the third-party libraries and extract it into
the subfolder "lib". After that, make sure you get the following structure of
folders:
/home/xxx/SudoDEM/
bin/
sudodem3d

lib/
3rdlibs/

sudodem/

share/
Step 2: Set the environmental variable PATH
Append the path of the executable sudodem3d, i.e., ‘/home/xxx/Sudo-
DEM/bin’ to the environmental variable PATH by adding the following line
to the config file ‘/home/xxx/.bashrc’:
1 PATH=${PATH}:/home/xxx/SudoDEM/bin

2.3.2 Non-root Compilation

(1) Directory trees
After unzipping the source package, you have the following directory tree:
/home/xxx/SudoDEM
SudoDEM2D
SudoDEM3D
scripts

INSTALL

LICENSE

[Link]
You may add subfolders, e.g.,‘3rdlib’, ‘build2d’, ‘build3d’, and ‘sudo-
deminstall’, inside the parent folder ‘SudoDEM’. Thus, the directory tree

10
A Quick Guide to SudoDEM CHAPTER 2. REQUIREMENTS

looks like
/home/xxx/SudoDEM
SudoDEM2D
SudoDEM3D
3rdlib
build2d
build3d
sudodeminstall
scripts

INSTALL

LICENSE

[Link]
The subfolder ‘3rdlib’ is for the third-party libraries. Note that the sub-
folder ‘3rdlib’ hosts the header files and the compiled dynamic libraries, while
for launching SudoDEM you need to copy all dynamic libraries ‘*.so.*’ files
into the install directory, e.g., ‘lib/3rdlibs/’. The compilation of the 2D and
3D versions will be done within the subfolders ‘build2d’ and ‘build3d’, re-
spectively. The compiled outputs (binary and library) will be installed in the
subfolder ‘sudodeminstall’.

(2) Compilation of third-party libraries

/home/xxx/SudoDEM/3rdlib
boost-1_6_7

eigen3.3.5

libQGLViewer-2.6.3

minieigen
Four major libraries:

• Boost-1.67
Download the source from its official page ([Link]
users/history/version_1_67_0.html) or the github repo (https:
//[Link]/SwaySZ/boost-1_6_7).

11
CHAPTER 2. REQUIREMENTS A Quick Guide to SudoDEM

1 cd boost−1\_6\_7
2 ./[Link] −−prefix=$PWD/../boost167 −−with−
libraries=python,thread,filesystem,iostreams,regex,
serialization,system,date_time link=shared runtime−link=
shared −−without−icu
3 ./b2 −j3
4 ./b2 install

• Eigen-3.3.5: no need to compile but all head files will be included.

You may download the source from the repo ([Link]
SwaySZ/Eigen-3.3.5).

• MiniEigen
Download the source from the repo ([Link]
minieigen). Set the directory of your compiled boost in the file ‘CMake-
[Link]’:
1 s e t (BOOST_ROOT "/home/ swayzhao / s o f t w a r e /DEM/3 r d l i b / b o o s t 1 6 7
")

Then,
1 mkdir build
2 cd build

3 cmake ../

4 make

After compilation, you will get the shared library ‘[Link]’. Copy
it to the folder ‘lib/3rdlibs/py’.

• LibQGLViewer-2.6.3
Download the source from the repo ([Link]
libQGLViewer-2.6.3)
1 cd QGLViewer
2 qmake
3 make

Tools and dependencies:

• biuld-essential, cmake

• freeglut3-dev

12
A Quick Guide to SudoDEM CHAPTER 2. REQUIREMENTS

• zlib1g-dev (boost)

• python-dev (boost)

• pyqt4-dev-tools

• qt4-default

• python-numpy python-tk

• libbz2-dev

• python-xlib python-qt4

• ipython3.0 python-matplotlib

• libxi-dev

• libglib2.0-dev

• libxmu-dev

(3) Compilation of SudoDEM main programs

Here we show the example to compile SudoDEM3D as follows. In the folder
‘SudoDEM3D’, you have the directory tree:
/home/xxx/SudoDEM/SudoDEM3D
cMake
[Link]
core
doc
gui

lib
pkg
py
Prior to compiling, you may need to edit the file ‘[Link]’, in
which you may change the installation directory, paths of library header
files, etc.
(a) The installation path is set by

13
 CHAPTER 2. REQUIREMENTS A Quick Guide to SudoDEM

1 SET(CMAKE_INSTALL_PREFIX " $ {CMAKE_CURRENT_SOURCE_DIR} / . . /

s u d o d e m i n s t a l l /SudoDEM3D" )

(b) The root directory of BOOST library is set by

1 s e t (BOOST_ROOT " $ {CMAKE_CURRENT_SOURCE_DIR} / . . / 3 r d l i b / b o o s t 1 6 7

")

(c) For the QGLViewer, the paths for both header files and library are set
by the following two lines:

1 s e t (QGLVIEWER_INCLUDE_DIR " $ {CMAKE_CURRENT_SOURCE_DIR} / . . / 3 r d l i b /

libQGLViewer − 2 . 6 . 3 / " )
2 s e t (QGLVIEWER_LIBRARIES " $ {CMAKE_CURRENT_SOURCE_DIR} / . . / 3 r d l i b /
libQGLViewer − 2 . 6 . 3 / QGLViewer/ libQGLViewer . s o " )

(d) Copy your system-installed numpy and pyqt4 to the path ’./lib/3rdlibs/py’.
You may also change the search path of numpy in the file [Link]
under the folder ‘cMake’:

1 s e t (DPDIR " $ {CMAKE_CURRENT_SOURCE_DIR} / . . / d e m 2 d i n s t a l l /SudoDEM/

l i b /3 r d l i b s /py / " )

Compile and install SudoDEM3D:

1 cd build3d
2 cmake ../SudoDEM3D
3 make −j3

4 make install

The binary ‘sudodem3d’ will be installed in the path ‘/home/xxx/Sudo-

DEM/sudodeminstall/SudoDEM3D/bin/’. Append the path to the environ-
mental variable PATH by adding the following line to the config file ‘/home-
/xxx/.bashrc’:

1 PATH=${PATH}:/home/xxx/SudoDEM/sudodeminstall/
SudoDEM3D/bin

You are expected to get the following directory tree after installation:

14
A Quick Guide to SudoDEM CHAPTER 2. REQUIREMENTS

/home/xxx/SudoDEM/sudodeminstall
bin
lib
3rdlibs
*.so
py
[Link]

numpy

PyQt4

other Python modules

sudodem
share
Note: the installation procedure of SudoDEM will overwrite the above
directories except ‘3rdlibs’. You need to copy all the 3rd-party libraries com-
piled above (boost, libQGLViewer, minieigen) to ‘3rdlibs’. If you encounter
any problems about the 3rd-libraries, you may check the directory tree for
the binary package of SudoDEM.
IMPORTANT: Add rpath to all 3rd-party libraries under ‘3rdlibs’ by
executing the script ‘[Link]’ in the folder ‘scripts’.
1 cd 3rdlibs
2 cp /home/xxx/SudoDEM/scripts/[Link] .
3 chmod +x [Link]

4 ./[Link]

Note: ‘[Link]’ will add rpath ‘$ORIGIN’ to all dynamic libraries

under ‘3rdlibs’ and rpath ‘$ORIGIN/../../’ to all dynamic libraries under
‘3rdlibs/py/PyQt4’. An rpath would make a dynamic library to search its
dependencies (other dynamic libraries) from the rpath-specified path in the
first place, which helps to avoid invoking other versions of dependencies in
the system path.

15
CHAPTER 2. REQUIREMENTS A Quick Guide to SudoDEM

16
 Chapter 3

Examples

We prepared several examples to show a quick start of using SudoDEM, and

the corresponding script for each example is available on the website or the
Github repository ([Link]

3.1 Example 0: run a simulation

This example will show the ingredients of a simulation by simulating a sphere
free falling onto a ground. We prepare a Python script named ‘[Link]’
at our work directory (e.g., /home/xxx/example’). Open a terminal, and
type ’sudodem3d [Link]’ if the work directory is at ‘/home/xxx/ex-
ample/’ otherwise providing a relative or absolute path of ‘[Link]’.
Here is the script of ‘[Link]’:
1 # A s p h e r e f r e e f a l l i n g on a ground under g r a v i t y
2
3 from sudodem import u t i l s # module u t i l s has some a u x i l i a r y
functions
4
5 #1 . we d e f i n e m a t e r i a l s f o r p a r t i c l e s
6 mat = RolFrictMat ( l a b e l="mat1" ,Kn=1e8 , Ks=7e7 , f r i c t i o n A n g l e=math .
atan ( 0 . 0 ) , d e n s i t y =2650) # m a t e r i a l 1 f o r p a r t i c l e s
7 wallmat1 = RolFrictMat ( l a b e l=" wallmat1 " ,Kn=1e8 , Ks=0, f r i c t i o n A n g l e
=0.) # m a t e r i a l 2 f o r w a l l s
8
9 # we add t h e m a t e r i a l s t o t h e c o n t a i n e r f o r m a t e r i a l s i n t h e
simulation .
10 # O can be r e g a r d e d a s an o b j e c t f o r t h e whole s i m u l a t i o n .
11 # We can s e e t h a t m a t e r i a l s i s a p r o p e r t y o f O, s o we u s e a dot
o p e r a t i o n (O. m a t e r i a l s ) t o g e t t h e p r o p e r t y ( m a t e r i a l s ) .
12 O. m a t e r i a l s . append ( mat ) # m a t e r i a l s i s a l i s t , i . e . , a c o n t a i n e r
f o r a l l m a t e r i a l s i n v o l v i n g in the s i m u l a t i o n s

17
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

13 O. m a t e r i a l s . append ( wallmat1 )
14 # adding a m a t e r i a l t o t h e c o n t a i n e r ( m a t e r i a l s ) s o t h a t we can
e a s i l y a c c e s s i t by i t s l a b e l . See t h e d e f i n i t i o n o f a w a l l
below .
15
16 #c r e a t e a s p h e r i c a l p a r t i c l e
17 sp = s p h e r e ( ( 0 , 0 , 1 0 . 0 ) , 1 . 0 , m a t e r i a l = mat )
18 O. b o d i e s . append ( sp )
19 # c r e a t e t h e ground a s an i n f i n i t e p l a n e w a l l
20 ground = u t i l s . w a l l ( 0 , a x i s =2, s e n s e =1, m a t e r i a l = ' wallmat1 ' ) #
u t i l s . w a l l d e f i n e s a w a l l with normal p a r a l l e l t o one a x i s o f
t h e g l o b a l C a r t e s i a n c o o r d i n a t e system . The f i r s t argument (
h e r e = 0 ) s p e c i f i e s t h e l o c a t i o n o f t h e w a l l , and with t h e
keyword a x i s ( h e r e = 2 , a x i s has a v a l u e o f 0 , 1 , 2 f o r x , y ,
z a x i e s r e s p e c t i v e l y ) we know t h a t t h e normal o f t h e w a l l i s
a l o n g z a x i s with c e n t e r i n g a t z = 0 ( 0 , s p e c i f i e d by t h e
f i r s t argument ) ; t h e s e c o n d keyword s e n s e ( h e r e = 1 , s e n s e
can be −1, 0 , 1 f o r n e g a t i v e , both , and p o s i t i v e s i d e s
r e s p e c t i v e l y ) s p e c i f i e s that the p o s i t i v e s i d e ( i . e . , the
normal p o i n t s t o t h e p o s i t i v e d i r e c t i o n o f t h e a x i s ) o f t h e
w a l l i s a c t i v a t e d , a t which we w i l l compute t h e i n t e r a c t i o n
o f a p a r t i c l e and t h i s w a l l ; t h e l a s t keyword m a t e r i a l s
s p e c i f i e s a m a t e r i a l t o t h i s w a l l , and we a s s i g n a s t r i n g ( i .
e . , t h e l a b e l o f wallmat1 ) t o t h e keyword , and d i r e c t l y
a s s i g n i n g t h e o b j e c t o f t h e m a t e r i a l ( i . e . , wallmat1 r e t u r n e d
by RolFrictMat ) i s a l s o a c c e p t a b l e .
21 O. b o d i e s . append ( ground ) # add t h e ground t o t h e body c o n t a i n e r
22

23
24 # d e f i n e a f u n c t i o n i n v o k e d p e r i o d i c a l l y by a PyRunner
25 d e f record_data ( ) :
26 #open a f i l e
27 f o u t = open ( ' data . dat ' , ' a ' ) #c r e a t e / open a f i l e named ' data .
dat ' i n append mode
28 #we want t o r e c o r d t h e p o s i t i o n o f t h e s p h e r e
29 p = O. b o d i e s [ 0 ] # t h e s p h e r e has been appended i n t o t h e l i s t
O. b o d i e s , and we can u s e t h e c o r r e s p o n d i n g i n d e x t o a c c e s s i t
. We have added two b o d i e s i n t o t a l t o O. b o d i e s , and t h e
s p h e r e i s t h e f i r s t one , i . e . , t h e i n d e x i s 0 .
30 # p i s the o b j e c t of the sphere
31 # a s a body o b j e c t , p has s e v e r a l p r o p e r t i e s , e . g . , s t a t e ,
m a t e r i a l , which can be l i s t e d by p . d i c t ( ) i n t h e commandline .
32 pos = p . s t a t e . pos # g e t t h e p o s i t i o n o f t h e s p h e r e
33 p r i n t >>f o u t , O. i t e r , pos [ 0 ] , pos [ 1 ] , pos [ 2 ] # we p r i n t t h e
i t e r a t i o n number and t h e p o s i t i o n ( x , y , z ) i n t o t h e f i l e f o u t .
34 # then , c l o s e t h e f i l e and r e l e a s e r e s o u r c e
35 fout . close ()
36
37 # c r e a t e a Newton e n g i n e by which t h e p a r t i c l e s move under Newton

18
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

' s Law
38 newton=N e w t o n I n t e g r a t o r ( damping = 0 . 1 , g r a v i t y = ( 0 . , 0 . 0 , − 9 . 8 ) , l a b e l
=" newton " )
39 # we s e t a l o c a l damping o f 0 . 1 and g r a v i t a t i o n a l a c c e l e r a t i o n
−9.8 m/ s 2 a l o n g z a x i s .
40
41 # the engine c o n t a i n e r i s the k e r n e l part of a simulation
42 # During each time s t e p , each e n g i n e w i l l run one by one f o r
c o m p l e t i n g a DEM c y c l e .
43 O. e n g i n e s =[
44 F o r c e R e s e t t e r ( ) , # f i r s t , we need t o r e s e t t h e f o r c e c o n t a i n e r
t o s t o r e upcoming c o n t a c t f o r c e s
45 # next , we conduct t h e broad phase o f c o n t a c t d e t e c t i o n by
comparing a x i s −a l i g n e d bounding boxs (AABBs) t o r u l e out
t h o s e p a i r s t h a t a r e d e f i n i t e l y not c o n t a c t i n g .
46 I n s e r t i o n S o r t C o l l i d e r ( [ Bo1_Sphere_Aabb ( ) , Bo1_Wall_Aabb ( ) ] ,
v e r l e t D i s t =0.2∗0.01) ,
47 # then , we e x e c u t e t h e narrow phase o f c o n t a c t d e t e c t i o n
48 I n t e r a c t i o n L o o p ( # we w i l l l o o p a l l p o t e n t i a l l y c o n t a c t i n g
pairs of p a r t i c l e s
49 [ Ig2_Sphere_Sphere_ScGeom ( ) , Ig2_Wall_Sphere_ScGeom ( ) ] , # f o r
d i f f e r e n t p a r t i c l e shapes , we w i l l c a l l t h e c o r r e s p o n d i n g
f u n c t i o n s t o compute t h e c o n t a c t geometry , f o r example ,
Ig2_Sphere_Sphere_ScGeom ( ) i s used t o compute t h e c o n t a c t
geometry o f a s p h e r e −s p h e r e p a i r .
50 [ Ip2_RolFrictMat_RolFrictMat_RolFrictPhys ( ) ] , # we w i l l
compute t h e p h y s i c a l p r o p e r t i e s o f c o n t a c t , e . g . , c o n t a c t
s t i f f n e s s and c o e f f i c i e n t o f f r i c t i o n , a c c o r d i n g t o t h e
physical properties of contacting p a r t i c l e s .
51 [ R o l l i n g R e s i s t a n c e L a w ( u s e _ r o l l i n g _ r e s i s t a n c e=F a l s e ) ] # Next
, we compute c o n t a c t f o r c e s i n terms o f t h e c o n t a c t law with
t h e i n f o ( c o n t a c t g e o m e t r i c and p h y s i c a l p r o p e r t i e s ) computed
a t l a s t two s t e p s .
52 ),
53 newton , # f i n a l l y , we compute a c c e l e r a t i o n and v e l o c i t y o f
each p a r t i c l e and update t h e i r p o s i t i o n s .
54 #a t each time s t e p , we have a PyRunner f u n c t i o n h e l p t o hack
i n t o a DEM c y c l e and do whatever you want , e . g . , c h a n g i n g
p a r t i c l e s ' s t a t e s , and s a v i n g some data , o r j u s t s t o p p i n g t h e
program a f t e r a c e r t a i n r u n n i n g time .
55 PyRunner ( command= ' record_data ( ) ' , i t e r P e r i o d =10000 , l a b e l= '
r e c o r d ' , dead = F a l s e )
56 ]
57
58 # we need t o s e t a time s t e p
59 O. dt = 1 e−5
60
61 # c l e a n t h e f i l e data . dat and w r i t e a f i l e head
62 f o u t = open ( ' data . dat ' , 'w ' ) # we open t h e f i l e i n a w r i t e mode s o

19
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

t h a t t h e f i l e w i l l be c l e a n
63 p r i n t >>f o u t , ' i t e r , x , y , z ' # w r i t e a f i l e head
64 fout . c l o s e ( ) # c l o s e the f i l e
65

66 # run t h e s i m u l a t i o n
67 # you have two c h o i c e s :
68 # 1 . g i v e a command here , l i k e ,
69 # O. run ( )
70 # you can s e t how many i t e r a t i o n s t o run
71 O. run ( 1 6 0 0 0 0 0 )
72 # 2 . c l i k e t h e run button a t t h e GUI c o n t r o l e r
Run the simulation by typing the following command in a terminal:
1 sudodem3d [Link]

For multi-threads, e.g., using 2 threads 1

1 sudodem3d −j2 [Link]

The recorded data ’[Link]’ looks like below:

1 iter , x, y, z
2 10000 0.0 0.0 9.95588676912
3 20000 0.0 0.0 9.82357353912
4 30000 0.0 0.0 9.60306030912
5 40000 0.0 0.0 9.29434707912
6 50000 0.0 0.0 8.89743384912
7 60000 0.0 0.0 8.41232061912
8 70000 0.0 0.0 7.83900738912
We can take a quick view on the recorded data using gnuplot. Install gnuplot
by typing the following command in a terminal:
1 sudo apt install gnuplot

Then, open gnuplot and plot the data as shown in Fig. 3.1:
1 $> gnuplot
2 gnuplot> plot '[Link]' using ($1/10000):4; set xlabel 'iterations
[x10000]'; set ylabel 'position z [m]'

3.2 Example 1: packing of super-ellipsoids

The surface function of a superellipsoid in the local Cartesian coordinates
can be defined as
 x  2  y  2  1  z  2
  1 +   1 2 +   2 = 1 (3.1)
rx ry rz
1
Note: it is not acceptable to use multi-threads for the exemplified case.

20
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

10
'[Link]' using ($1/10000):4

6
position z [m]

0
0 20 40 60 80 100 120 140 160
iterations [x10000]

Figure 3.1: The z-position of a sphere free falling onto a ground.

where rx , ry and rz are referred to as the semi-major axis lengths in the

direction of x, y, and z axies, respectively; and i (i = 1, 2) are the shape pa-
rameters determining the sharpness of particle edges or squareness of particle
surface. Varying i between 0 and 2 yields a wide range of convex-shaped
superellipsoid. Two functions for generating a superellipsoid are highlighted
below (see Sec. 5.1.2):

• NewSuperquadrics2(rx , ry , rz , 1 , 2 , mat, rotate, isSphere)

• NewSuperquadrics_rot2(rx , ry , rz , 1 , 2 , mat, qw , qx , qy , qz , isSphere)

Here we give an example of a simulation of superellipsoids free falling into

a cubic box.
1 ###############################################
2 # g r a n u l a r p a c k i n g o f super−e l l i p s o i d s
3 # We p o s i t i o n p a r t i c l e s a t a l a t t i c e g r i d ,
4 # then p a r t i c l e s f r e e f a l l i n t o a c u b i c box .
5 ###############################################
6 # import some modules
7 from sudodem import _ s u p e r q u a d r i c s _ u t i l s
8
9 from sudodem . _ s u p e r q u a d r i c s _ u t i l s import ∗

21
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

2.0
1.5
2.0
1.0 1.5
1.0
0.8
0.8
0.2 0.2
Figure 3.2: Superellipsoids.

10 import math
11 import random a s rand
12 import numpy a s np
13
14 ########d e f i n e some p a r a m e t e r s#################
15 i s S p h e r e=F a l s e
16
17 num_x = 5
18 num_y = 5
19 num_z = 20
20 R = 0.1
21
22 #######d e f i n e some a u x i l i a r y f u n c t i o n s#########
23

24 #d e f i n e a l a t i c e g r i d
25 #R: d i s t a n c e between two n e i g h b o r i n g nodes
26 #num_x : number o f nodes a l o n g x a x i s
27 #num_y : number o f nodes a l o n g y a x i s
28 #num_z : number o f nodes a l o n g z a x i s
29 #r e t u r n a l i s t o f t h e p o s t i o n s o f a l l nodes
30 d e f G r i d I n i t i a l (R, num_x=10 ,num_y=10 ,num_z=20) :
31 pos = l i s t ( )
32 f o r i i n r a n g e (num_x) :

22
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

33 f o r j i n r a n g e (num_y) :
34 f o r k i n r a n g e (num_z) :
35 x = i ∗R∗ 2 . 0
36 y = j ∗R∗ 2 . 0
37 z = k∗R∗ 2 . 0
38 pos . append ( [ x , y , z ] )
39 r e t u r n pos
40
41 #g e n e r a t e a sample
42 d e f GenSample ( r , pos ) :
43 f o r p i n pos :
44 e p s i l o n 1 = rand . uniform ( 0 . 7 , 1 . 6 )
45 e p s i l o n 2 = rand . uniform ( 0 . 7 , 1 . 6 )
46 a = r
47 b = r ∗ rand . uniform ( 0 . 4 , 0 . 9 )#a s p e c t r a t i o
48 c = r ∗ rand . uniform ( 0 . 4 , 0 . 9 )
49
50 body = NewSuperquadrics2 ( a , b , c , e p s i l o n 1 , e p s i l o n 2 , p_mat , True ,
isSphere )
51 body . s t a t e . pos=p
52 O. b o d i e s . append ( body )
53

54 #########s e t u p a s i m u l a t i o n####################
55 # material for p a r t i c l e s
56 p_mat = SuperquadricsMat ( l a b e l="mat1" ,Kn=1e5 , Ks=7e4 , f r i c t i o n A n g l e
=math . atan ( 0 . 3 ) , d e n s i t y =2650 , betan =0, b e t a s =0)
57 # betan and b e t a s a r e c o e f f i c i e n t s o f v i s c o u s damping a t c o n t a c t ,
no v i s c o u s damping with 0 by d e f a u l t .
58 # material f o r the s i d e w a l l s
59 wall_mat = SuperquadricsMat ( l a b e l=" wallmat " ,Kn=1e6 , Ks=7e5 ,
f r i c t i o n A n g l e =0.0 , betan =0, b e t a s =0)
60 # m a t e r i a l f o r t h e bottom w a l l
61 wallmat_b = SuperquadricsMat ( l a b e l=" wallmat " ,Kn=1e6 , Ks=7e5 ,
f r i c t i o n A n g l e=math . atan ( 1 ) , betan =0, b e t a s =0)
62 # add m a t e r i a l s t o O. m a t e r i a l s
63 O. m a t e r i a l s . append (p_mat)
64 O. m a t e r i a l s . append ( wall_mat )
65 O. m a t e r i a l s . append ( wallmat_b )
66
67 # c r e a t e t h e box and add i t t o O. b o d i e s
68 O. b o d i e s . append ( u t i l s . w a l l (−R, a x i s =0, s e n s e =1, m a t e r i a l = wall_mat
) )#l e f t w a l l a l o n g x a x i s
69 O. b o d i e s . append ( u t i l s . w a l l ( 2 . 0 ∗R∗num_x−R, a x i s =0, s e n s e =−1,
m a t e r i a l = wall_mat ) )#r i g h t w a l l a l o n g x a x i s
70 O. b o d i e s . append ( u t i l s . w a l l (−R, a x i s =1, s e n s e =1, m a t e r i a l = wall_mat
) )#f r o n t w a l l a l o n g y a x i s
71 O. b o d i e s . append ( u t i l s . w a l l ( 2 . 0 ∗R∗num_y−R, a x i s =1, s e n s e =−1,
m a t e r i a l = wall_mat ) )#back w a l l a l o n g y a x i s
72 O. b o d i e s . append ( u t i l s . w a l l (−R, a x i s =2, s e n s e =1, m a t e r i a l =wallmat_b

23
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

) )#bottom w a l l a l o n g z a x i s
73
74 # create a l a t t i c e grid
75 pos = G r i d I n i t i a l (R, num_x=num_x, num_y=num_y, num_z=num_z) # g e t
p o s i t i o n s o f a l l nodes
76 # c r e a t e p a r t i c l e s a t each nodes
77 GenSample (R, pos )
78
79 # create engines
80 # d e f i n e a Newton e n g i n e
81 newton=N e w t o n I n t e g r a t o r ( damping = 0 . 1 , g r a v i t y = ( 0 . , 0 . , − 9 . 8 ) , l a b e l=
" newton " , i s S u p e r q u a d r i c s =1) # i s S u p e r q u a d r i c s : 1 f o r
superquadrics
82
83 O. e n g i n e s =[
84 F o r c e R e s e t t e r ( ) , I n s e r t i o n S o r t C o l l i d e r ( [ Bo1_Superquadrics_Aabb
( ) , Bo1_Wall_Aabb ( ) ] , v e r l e t D i s t = 0 . 2 ∗ 0 . 1 ) ,
85 InteractionLoop (
86 [ Ig2_Wall_Superquadrics_SuperquadricsGeom ( ) ,
Ig2_Superquadrics_Superquadrics_SuperquadricsGeom ( ) ] ,
87 [ Ip2_SuperquadricsMat_SuperquadricsMat_SuperquadricsPhys ( )
] , # c o l l i s i o n " physics "
88 [ SuperquadricsLaw ( ) ] # c o n t a c t law
89 ),
90 newton
91
92 ]
93 O. dt=5e−5
With a GUI controller, we can click ’show 3D’ button at the panel ’Sim-
ulation’ to display the simulating scene. On the ’Display’ panel, select the
render ’Gl1_Superquadrics’ to configure the resolution of solid or wireframe
particles, as shown in Fig. 3.3.
Figs. 3.4 - 3.6 show simulations of packing of superellipsoids at different
states for the presented example.
Several properties and methods for a shape Superquadrics are listed be-
low:

• Properties:
• rx , ry , rz : float, semi-major axis lengths in x, y, and z axies.
• 1 , 2 : float, shape parameters.
• isSphere: bool, particle is spherical or not.
• Methods:
• getVolume(): float, return particle’s volume.

24
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

Figure 3.3: GUI for display configuration of superellipsoids.

• getrxyz(): Vector3, return particle’s rxyz.

• geteps(): Vector2, return particle’s eps1 and eps2.

Here is an example to get the volume of a particle with an id of 100:

1 volume = O. b o d i e s [ 1 0 0 ] . shape . getVolume ( )

A general method to list all properties of an object is as follows:

1 O. b o d i e s [ 1 0 0 ] . d i c t ( )

The output is like follows:

1{'bound': <Aabb instance at 0x56072ec83da0>,
2 'chain': −1,
3 'clumpId': −1,

4 'flags': 3,

5 'groupMask': 1,

6 'id': 100,

7 'iterBorn': 0,

8 'material': <SuperquadricsMat instance at 0x56072e6ed5f0>,

25
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

Figure 3.4: Initial packing of superellipsoids.

9 'shape': <Superquadrics instance at 0x56072ec83b60>,

10 'state': <State instance at 0x56072ec839e0>,
11 'timeBorn': 0.0}

We can see that the properties are printed in a dict form. Note that the value
enclosed by pointy brackets is another object (actually, instance of an object
with a specified memory address), which means that we can further access it
using the ".dict()" method. Again, the properties of the object State can be
accessed by:
1 O. b o d i e s [ 1 0 0 ] . s t a t e . d i c t ( )

The output is as follows:

1 {'angMom': Vector3(0,0,0),
2 'angVel': Vector3(0,0,0),
3 'blockedDOFs': 0,

4 'densityScaling': 1.0,

5 'inertia': Vector3(0.0045389863901528615,0.007287422614611933,

6 0.0067053718092146865),
7 'isDamped': True,

8 'mass': 3.225715855242557,

9 'refOri': Quaternion((1,0,0),0),

10 'refPos': Vector3(0,0,0),

26
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

Figure 3.5: Packing of superellipsoids during falling.

11 'se3': (Vector3(0,0.8,3),
12 Quaternion((0.3970010520699211,0.5339678220536823,
13 −0.7465042060609055),2.1754719537556495)),
14 'vel': Vector3(0,0,0)}

For an object Shape, we list all properties which may vary for different child
classes.
1 O. b o d i e s [ 1 0 0 ] . shape . d i c t ( )

The output is as follows:

1{'color': Vector3(0.6407663308273844,0.07426042997942327,
2 0.05872324344642611),
3 'eps': Vector2(1.3975848801318045,1.5376309088540607),

4 'eps1': 1.3975848801318045,

5 'eps2': 1.5376309088540607,

6 'highlight': False,

7 'isSphere': False,

8 'rx': 0.1,

9 'rxyz': Vector3(0.1,0.06469580193313924,0.07572423080016706),

10 'ry': 0.06469580193313924,

11 'rz': 0.07572423080016706,

12 'wire': False}

27
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

Figure 3.6: Final packing of superellipsoids.

Here we give an example to list ids and positions of all superellipsoids:

1 f o r b i n O. b o d i e s : # l o o p a l l b o d i e s i n t h e simulation
2 # we have t o check i f t h e p r e s e n t body i s a s u p e r e l l i p s o i d or
others , e . g . , a wall
3 i f i s i n s t a n c e ( b . shape , S u p e r q u a d r i c s ) : # if it is a
s u p e r e l l i p s o i d , then t o do
4 pid = b . id # get the id of p a r t i c l e b
5 pos = b . s t a t e . pos # g e t t h e p o s i t i o n o f particle b
6 p r i n t pid , pos

3.3 Example 2: triaxial tests of super-ellipsoids

This example shows a triaxial compression simulation of superellipsoids using
an Engine CompressionEngine. The users can also customized their engines
or compression procedures easily via a PyRunner Engine. Here is a basic
setup of CompressionEngine for consolidation and compression below.

• Consolidation
1 t r i a x=CompressionEngine (
2 g o a l x = 1 . 0 e5 , # c o n f i n i n g s t r e s s 100 kPa
3 g o a l y = 1 . 0 e5 ,

28
A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

4 g o a l z = 1 . 0 e5 ,
5 savedata_interval = 10000 ,
6 echo_interval = 2000 ,
7 continueFlag = False ,
8 max_vel = 1 0 0 . 0 ,
9 gain_alpha = 0 . 5 ,
10 f _ t h r e s h o l d = 0 . 0 1 , #1− f / g o a l
11 f m i n _ t h r e s h o l d = 0 . 0 1 , #t h e r a t i o o f d e v i a t o r i c s t r e s s
t o mean s t r e s s
12 unbf_tol = 0 . 0 1 , #unblanced f o r c e r a t i o
13 f i l e = ' r e c o r d −con '
14 )

• goalx = goaly = goalz = confining stress

• savedata_interval : the interval of DEM iterations for saving some
data (Axial Strain, Volumetric Strain, Mean Stress and Devia-
torStress) into a file specified by file. Deprecated for consolida-
tion.
• echo_interval : the interval of DEM iterations for outputting in-
formation (Iterations, Unbalanced force ratio, stresses in x, y, z
directions) in the terminal
• continueFlag: reserved keyword for continuing consolidation or
compression.
• max_vel : the maximum velocity of a wall, which is useful to con-
strain the movement of a wall during confining.
• gain_alpha α: guiding the movement of the confining walls.
fw − fg
vw = α P (3.2)
∆t kn
where vw is the velocity of a confining wall at the next time step;
fw and fg are the present Pand goal confining forces, respectively;
∆t is the time step, and kn is the summation of normal contact
stiffness from all particle-wall contacts on this wall.
• f_threshold αf : a target ratio of the deviation of the present stress
from the goal to the goal, i.e.,
|fw − fg |
αf = (3.3)
fg
• fmin_threshold βf : a target ratio of the deviatoric stress q to the
mean stress p, i.e.,
r
q σii 3 0 0
βf = , p = , q= σ σ (3.4)
p 3 2 ij ij

29
CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

1 X cc 0
σij = f l, σij = σij − pδij (3.5)
V c∈V i j
0
where σij is the stress tensor with an assembly, and σij is its cor-
responding deviatoric part; V is the total volume of the specimen;
f c and lc are the contact force and branch vector at contact c,
respectively, and δij is the Kronecker delta.
• unbf_tol γf : a target unbalanced force ratio defined as

fjc + fjb |
P P
|
Bi ∈V j∈Bi
γf = P (3.6)
2 |f c |
c∈V

Note: the consolidation ends up with that all αf , βf and γf reach

the targets.
• file: the file name for the recording file

• Compression
1 t r i a x=CompressionEngine (
2 z_servo = F a l s e ,
3 g o a l x = 1 . 0 e5 ,
4 g o a l y = 1 . 0 e5 ,
5 goalz = 0.01 ,
6 ramp_interval = 1 0 ,
7 ramp_chunks = 2 ,
8 savedata_interval = 1000 ,
9 echo_interval = 2000 ,
10 max_vel = 1 0 . 0 ,
11 gain_alpha = 0 . 5 ,
12 f i l e=' sheardata ' ,
13 target_strain = 0.4 ,
14 gain_alpha = 0 . 5
15 )

• z_servo: a False value will trigger a strain (or displacement) con-

trol in the z direction, and in this case goalz needs a strain rate.
• ramp_interval : iterations needed for reaching a constant strain
rate. This keyword is designed to mimic a continuing loading at
laboratory, but it seems not too valued and setting a small value
to kindly deprecate it.

30
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

• ramp_chunks: number of intervals to reach the constant strain

rate combined with the keyword ramp_interval. It will be depre-
cated as does ramp_interval.
• target_strain: at this level of strain the shear procedures ends.
The strain is ogarithmic strain.
Note: to ensure quasi-static behavior during shear, the shear
strain rate should be sufficiently small. Specifically, the inertial
number I should be maintained less than 10−3

dˆ
I = ˙1 p (3.7)
σ0 /ρ

where ˙1 = d1 /dt is the major principal strain rate; dˆ and ρ are
the average diameter and material density of particles respectively;
σ0 is the confining stress.

As a demonstration, we first give an example of consolidation as follows:

1 ##################################################
2 from sudodem import p l o t , _ s u p e r q u a d r i c s _ u t i l s
3 from sudodem . _ s u p e r q u a d r i c s _ u t i l s import ∗
4
5 import math
6 import random a s rand
7 rand . s e e d ( 1 1 2 4 5 6 4 2 )
8 ####################################################
9 epsilon = 1.4
10 isSphere = False
11 b o x s i z e = [ 1 6 . 2 e −3 ,16.2 e −3 ,16.2 e −3]
12 ############################
13 #some a u x i l i a r y f u n c t i o n s
14 #g e n e r a t e a sample o f m o n o d i s p e r s e s u p e r e l l i p s o i d s with random
o r i e n t a i o n s and p o s i t i o n s
15 d e f gen_sample ( width , depth , h e i g h t , r_max = 1 . 0 e −3) :
16 b_num_tot = 5000
17 f o r j i n r a n g e ( b_num_tot ) :
18 rb = r_max ∗ 0 . 5
19 x = rand . uniform ( rb , width−rb )
20 y = rand . uniform ( rb , depth − rb )
21 z = rand . uniform ( rb , h e i g h t − rb )
22 r = r_max ∗ 0 . 5
23 b = NewSuperquadrics2 ( r ∗ 1 . 2 5 , r , r , e p s i l o n , e p s i l o n , mat , True ,
isSphere )
24 b . s t a t e . pos=(x , y , z )
25 O. b o d i e s . append ( b )
26 #M a t e r i a l s d e f i n i t i o n

31
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

27 mat = SuperquadricsMat ( l a b e l="mat1" ,Kn=3e4 , Ks=3e4 , f r i c t i o n A n g l e=

math . atan ( 0 . 1 ) , d e n s i t y =2650 e4 ) #d e f i n e M a t e r i a l with d e f a u l t
values
28 wallmat1 = SuperquadricsMat ( l a b e l=" wallmat1 " ,Kn=1e6 , Ks =0. ,
f r i c t i o n A n g l e =0.) #d e f i n e M a t e r i a l with d e f a u l t v a l u e s
29 wallmat2 = SuperquadricsMat ( l a b e l=" wallmat2 " ,Kn=1e6 , Ks =0. ,
f r i c t i o n A n g l e =0.) #d e f i n e M a t e r i a l with d e f a u l t v a l u e s
30 O. m a t e r i a l s . append ( mat )
31 O. m a t e r i a l s . append ( wallmat1 )
32 O. m a t e r i a l s . append ( wallmat2 )
33 #############################################
34 #####c r e a t e c o n f i n i n g w a l l s
35 ############################################
36 O. b o d i e s . append ( u t i l s . w a l l ( 0 , a x i s =0, s e n s e =1, m a t e r i a l = ' wallmat1
' ) )#l e f t x
37 O. b o d i e s . append ( u t i l s . w a l l ( b o x s i z e [ 0 ] , a x i s =0, s e n s e =−1, m a t e r i a l =
' wallmat1 ' ) )#r i g h t x
38 O. b o d i e s . append ( u t i l s . w a l l ( 0 , a x i s =1, s e n s e =1, m a t e r i a l = ' wallmat1
' ) )#f r o n t y
39 O. b o d i e s . append ( u t i l s . w a l l ( b o x s i z e [ 1 ] , a x i s =1, s e n s e =−1, m a t e r i a l =
' wallmat1 ' ) )#back y
40

41 O. b o d i e s . append ( u t i l s . w a l l ( 0 , a x i s =2, s e n s e =1, m a t e r i a l = ' wallmat2

' ) )#bottom z
42 O. b o d i e s . append ( u t i l s . w a l l ( b o x s i z e [ 2 ] , a x i s =2, s e n s e =−1, m a t e r i a l =
' wallmat2 ' ) )#top z
43 #####c r e a t e p a r t i c l e s
44

45 gen_sample ( b o x s i z e [ 0 ] , b o x s i z e [ 1 ] , b o x s i z e [ 2 ] )
46
47 ####f u n c t i o n f o r s a v i n g data ( e . g . , s i m u l a t i o n s t a t e s h e r e )
48 def savedata () :
49 O. s a v e ( s t r (O. time )+ ' . xml . bz2 ' )
50

51 #########################################
52 ####t r i a x e n g i n e
53 #CAUTION: t h e CompressionEngine r e g a r d s t h e f i r s t s i x b o d i e s a s
c o n f i n i n g w a l l s by d e f a u l t . So t h e w a l l s s h o u l d be g e n e r a t e d
first .
54 t r i a x=CompressionEngine (
55 g o a l x = 1 . 0 e5 , # c o n f i n i n g s t r e s s 100 kPa
56 g o a l y = 1 . 0 e5 ,
57 g o a l z = 1 . 0 e5 ,
58 savedata_interval = 10000 ,
59 echo_interval = 2000 ,
60 continueFlag = False ,
61 max_vel = 1 0 . 0 ,
62 gain_alpha = 0 . 5 ,
63 f _ t h r e s h o l d = 0 . 0 1 , #1− f / g o a l

32
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

64 f m i n _ t h r e s h o l d = 0 . 0 1 , #t h e r a t i o o f d e v i a t o r i c s t r e s s t o mean
stress
65 unbf_tol = 0 . 0 1 , #unblanced f o r c e r a t i o
66 f i l e = ' r e c o r d −con ' ,
67 )
68
69 ############################
70
71
72

73 newton=N e w t o n I n t e g r a t o r ( damping = 0 . 3 , g r a v i t y = ( 0 . , 0 . , 0 . ) , l a b e l="

newton " , i s S u p e r q u a d r i c s =1)
74 O. e n g i n e s =[
75 ForceResetter () ,
76 I n s e r t i o n S o r t C o l l i d e r ( [ Bo1_Superquadrics_Aabb ( ) , Bo1_Wall_Aabb
( ) ] , v e r l e t D i s t =0.1 e −3) ,
77 InteractionLoop (
78 [ Ig2_Wall_Superquadrics_SuperquadricsGeom ( ) ,
Ig2_Superquadrics_Superquadrics_SuperquadricsGeom ( ) ] ,
79 [ Ip2_SuperquadricsMat_SuperquadricsMat_SuperquadricsPhys ( )
],
80 [ SuperquadricsLaw ( ) ]
81 ),
82 triax ,
83 newton ,
84 PyRunner ( command= ' q u i e t _ b a l l ( ) ' , i t e r P e r i o d =2, i t e r L a s t = 5 0 0 0 0 ,
l a b e l= ' calm ' , dead = F a l s e ) ,
85 PyRunner ( command= ' s a v e d a t a ( ) ' , r e a l P e r i o d =7200.0 , l a b e l= '
s a v e d a t a ' , dead = F a l s e )#s a v i n g data e v e r y two h o u r s
86 ]
87
88
89 #d e l e t e t h e b a l l s o u t s i d e t h e box
90 def del_balls () :
91 num_del = 0
92 #w a l l p o s i t i o n
93 l e f t = O. b o d i e s [ 0 ] . s t a t e . pos [ 0 ]#x
94 r i g h t = O. b o d i e s [ 1 ] . s t a t e . pos [ 0 ]
95 f r o n t = O. b o d i e s [ 2 ] . s t a t e . pos [ 1 ]#y
96 back = O. b o d i e s [ 3 ] . s t a t e . pos [ 1 ]
97 bottom = O. b o d i e s [ 4 ] . s t a t e . pos [ 2 ]#z
98 top = O. b o d i e s [ 5 ] . s t a t e . pos [ 2 ]
99 f o r b i n O. b o d i e s :
100 i f i s i n s t a n c e ( b . shape , S u p e r q u a d r i c s ) :
101 ( x , y , z ) = b . s t a t e . pos
102 i f x > r i g h t or x < l e f t or y < f r o n t or
y > back o r z < bottom o r z > top :
103 O. b o d i e s . e r a s e ( b . i d )
104 num_del += 1

33
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

105 p r i n t s t r ( num_del )+" p a r t i c l e s a r e d e l e t e d ! "

106
107
108

109 def quiet_ball () :

110 g l o b a l calm_num
111
112 newton . q u i e t _ s y s t e m _ f l a g = True
113 i f calm_num > 2 0 0 0 :
114 O. e n g i n e s [ − 2 ] . i t e r P e r i o d= 5
115 calm_num += 5
116 else :
117 calm_num += 2
118 i f calm_num > 4 0 0 0 0 :
119 i f calm_num < 4 0 0 1 0 :
120 p r i n t " calm p r o c e d u r e i s o v e r "
121 del_balls ()
122 i f calm_num > 5 0 0 0 0 :
123 O. e n g i n e s [ − 2 ] . dead = True
124 O. s a v e ( " i n i t _ a s s e m b l y . xml . bz2 " )
125 #c o n s o l i d a t i o n b e g i n s
126 t r i a x . dead=F a l s e
127 print " consolidation begins "
128
129
130 O. dt=1e−4
131
132

133 t r i a x . dead=True
134 calm_num = 0

The output in the terminal is akin to the follows:

1 calm procedure is over
2 0 particles are deleted!
3 consolidation begins

4 Iter 50000 Ubf 0.483161, Ss_x:9.67986e−05, Ss_y:4.15364e−05,

Ss_z:6.1184e−05, e:0.992618
5 Iter 52000 Ubf 0.461741, Ss_x:1.2412, Ss_y:1.06812, Ss_z:1.10116,

e:0.854445
6 Iter 54000 Ubf 0.263699, Ss_x:1.75334, Ss_y:1.98567, Ss_z

:1.62181, e:0.757469
7 Iter 56000 Ubf 0.118329, Ss_x:3.06302, Ss_y:3.54818, Ss_z

:3.22986, e:0.679106
8 Iter 58000 Ubf 0.0110747, Ss_x:22.0177, Ss_y:21.3303, Ss_z

:21.3485, e:0.614824
9 Iter 60000 Ubf 0.000577066, Ss_x:91.4199, Ss_y:89.6039, Ss_z

:89.1057, e:0.58984

34
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

10 Iter 62000 Ubf 0.000128293, Ss_x:100.012, Ss_y:99.7251, Ss_z

:99.6123, e:0.587807
11 consolidation completed!

Figure 3.7: Packing of superellipsoids after consolidation.

After consolidation (see the snapshot in Fig. 3.7), we obtain a state file
by
1 O. s a v e ( " f i n a l _ c o n . xml . bz2 " )

Then, the consolidated assembly is subjected to triaxial compression with

the following scripts.
1 ##################################################
2 from sudodem import p l o t , _ s u p e r q u a d r i c s _ u t i l s
3
4 from sudodem . _ s u p e r q u a d r i c s _ u t i l s import ∗
5
6 import math
7
8 saving_num = 0
9 ####################################################
10 def savedata () :
11 g l o b a l saving_num
12 saving_num += 1
13 #p r i n t " t e s t i n g , no data saved "

35
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

14 O. s a v e ( ' s h e a r '+s t r ( saving_num )+ ' . xml . bz2 ' )

15
16 d e f quiet_system ( ) :
17 f o r i i n O. b o d i e s :
18 i . state . vel =(0. ,0. ,0.)
19 i . s t a t e . angVel = ( 0 . , 0 . , 0 . )
20 s u f f i x=" "
21 O. l o a d ( " f i n a l _ c o n "+s u f f i x+" . xml . bz2 " )
22 t r i a x = O. e n g i n e s [ 3 ]
23 O. dt = 5 e−5
24 f r i c t i o n = 0.5
25 O. m a t e r i a l s [ 0 ] . f r i c t i o n A n g l e=math . atan ( f r i c t i o n )
26 i f 1:
27 f o r i t r i n O. i n t e r a c t i o n s :
28 i d = min ( i t r . id1 , i t r . i d 2 )
29 i f id > 5:
30 i t r . phys . t a n g e n s O f F r i c t i o n A n g l e = f r i c t i o n
31 #i t r . phys . betan = 0 . 5
32
33 #########################
34 #s h e a r b e g i n s
35 #########################
36 def shear () :
37 quiet_system ( )
38 t r i a x . z_servo = F a l s e
39 t r i a x . ramp_interval = 10
40 t r i a x . ramp_chunks = 2
41 t r i a x . f i l e =" s h e a r d a t a "+s u f f i x
42 triax . goalz = 0.01
43 t r i a x . g o a l x = 1 e5
44 t r i a x . g o a l y = 1 e5
45 t r i a x . s a v e d a t a _ i n t e r v a l = 2500
46 t r i a x . e c h o _ i n t e r v a l = 2000
47 triax . target_strain = 0.4
48
49 shear ()
50
51 O. e n g i n e s [ − 1 ] . i t e r P e r i o d = 2500
52 O. e n g i n e s [ − 1 ] . r e a l P e r i o d = 0 . 0

The compression ends up with an axial strain of 40%, and the final con-
figuration of particles is shown in Fig. 3.8.
The macroscopic mechanical response of a granular material during shear-
ing can be quantified by the stress tensor that is calculated from discrete
measurements with the following definition:

1 X cc
σij = f l (3.8)
V c∈V i j

36
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

Figure 3.8: Packing of superellipsoids after shear.

where V is the total volume of the assembly, f c is the contact force at the
contact c, and lc is the branch vector joining the the centres of the two
contacting particles at contact c. The mean stress p and the deviatoric stress
q are defined as: r
1 3 0 0
p = σii , q = σ σ (3.9)
3 2 ij ij
0
where σij is the deviatoric part of stress tensor σij .
Given that the cubical specimens are confined by rigid walls, the axial
strain z and the volumetric strain v can be approximately calculated from
the positions of the boundary walls, i.e.,
Z H Z V
dh H0 dv V0
z = = − ln , v = x + y + z = = − ln (3.10)
H0 h H V0 v V

where H and V are the height and volume of the specimen during shearing,
respectively, and H0 and V0 are their initial values before shear. Positive
values of volumetric strain represent dilatancy.
The recorded data ’sheardata’ looks like below:
1 A x i a l S t r a i n V o l S t r a i n m ea n S tr e s s d e v i a t o r S t r e s s
2 0.000511000266 0.000490480823 106090.98 15931.56
3 0.003011001516 0.002450117535 131061.41 85718.91

37
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

4 0.005511002766 0.003867756447 148076.65 131980.95

5 0.008011004016 0.004825115794 159924.28 163706.90
6 0.010511005266 0.005409561745 168606.27 186520.14
7 0.013011006516 0.005673777079 175131.53 203243.28
8 0.015511007766 0.005650846281 180293.19 216248.17
9 0.018011009016 0.005393122254 184581.23 227010.23

Then we can plot the stress-strain variation and the loading path as shown
in Figs. 3.9 and 3.10.

1.50 15

1.25
Deviatoric stress ratio, q/p []

Volumetric strain, εv [%]

10
1.00

0.75
5
0.50

0.25 q/p
εv 0
0.00
0 10 20 30 40
Axial strain, [%]

Figure 3.9: Variation of deviatoric stress ratio q/p and volumetric strain v
with axial strain during shearing.

During shearing, a series of states (i.e., files "[Link].bz2") has been

saved sequentially so that it is readily to access other data for post-processing,
e.g., variation of mean coordination number, anisotropy and so forth.

3.4 Example 3: packing of GJKparticles

There are five kinds of shapes, namely sphere, polyhedron, cone, cylinder
and cuboid, in SudoDEM, for which the GJK algorithm of contact detection
is performed 2 , so that we coined a name GJKparticle for grouping these
particle shapes. In the Python module _gjkparticle_utils (see Sec. 5.1.3),
several construction functions are provided including:
2
GJKparticle is just experimental, and has not been rigorously tested yet.

38
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

300

250
Deviatoric stress, q [kPa]

200

150

100

50

0
0 50 100 150 200 250 300
Mean stress, p [kPa]

Figure 3.10: Loading path during shearing.

• GJKSphere(radius,margin,mat)

• GJKPolyhedron(vertices,extent,margin,mat, rotate)

• GJKCone(radius, height, margin, mat, rotate)

• GJKCylinder(radius, height, margin, mat, rotate)

• GJKCuboid(extent, margin,mat, rotate)

Note: a margin is used to expand a particle for achieving a better compu-

tational efficiency, which however should be large enough for ensuring that
collision occurs in this small region but small enough for shape fidelity.
Similar to Example 1, we conduct some simulations of particles free falling
into a cubic box. We first import some modules:
1 from sudodem import p l o t , _ g j k p a r t i c l e _ u t i l s
2 from sudodem import ∗
3

4 from sudodem . _ g j k p a r t i c l e _ u t i l s import ∗

5 import random a s rand
6 import math

Next, define some constants:

39
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

1 cube_v0 =
[[0 ,0 ,0] ,[1 ,0 ,0] ,[1 ,1 ,0] ,[0 ,1 ,0] ,[0 ,0 ,1] ,[1 ,0 ,1] ,[1 ,1 ,1] ,[0 ,1 ,1]]

2 cube_v = [ [ i ∗ 0 . 0 1 f o r i i n j ] f o r j i n cube_v0 ]
3
4 i s S p h e r e=F a l s e
5
6 num_x = 5
7 num_y = 5
8 num_z = 20
9 R = 0.01

Define a lattice grid:

1 #R: d i s t a n c e between two n e i g h b o r i n g nodes
2 #num_x : number o f nodes a l o n g x a x i s
3 #num_y : number o f nodes a l o n g y a x i s
4 #num_z : number o f nodes a l o n g z a x i s
5 #r e t u r n a l i s t o f t h e p o s t i o n s o f a l l nodes
6 d e f G r i d I n i t i a l (R, num_x=10 ,num_y=10 ,num_z=20) :
7 pos = l i s t ( )
8 f o r i i n r a n g e (num_x) :
9 f o r j i n r a n g e (num_y) :
10 f o r k i n r a n g e (num_z) :
11 x = i ∗R∗ 2 . 0
12 y = j ∗R∗ 2 . 0
13 z = k∗R∗ 2 . 0
14 pos . append ( [ x , y , z ] )
15 r e t u r n pos

Set properties of materials:

1 p_mat = GJKParticleMat ( l a b e l="mat1" ,Kn=1e4 , Ks=7e3 , f r i c t i o n A n g l e=
math . atan ( 0 . 5 ) , d e n s i t y =2650 , betan =0, b e t a s =0) #d e f i n e M a t e r i a l
with d e f a u l t v a l u e s
2 #mat = GJKParticleMat ( l a b e l ="mat1 " ,Kn=1e6 , Ks=1e5 , f r i c t i o n A n g l e
=0.5 , d e n s i t y =2650) #d e f i n e M a t e r i a l with d e f a u l t v a l u e s
3 wall_mat = GJKParticleMat ( l a b e l=" wallmat " ,Kn=1e4 , Ks=7e3 ,
f r i c t i o n A n g l e =0.0 , betan =0, b e t a s =0) #d e f i n e M a t e r i a l with
default values
4 wallmat_b = GJKParticleMat ( l a b e l=" wallmat " ,Kn=1e4 , Ks=7e3 ,
f r i c t i o n A n g l e=math . atan ( 1 ) , betan =0, b e t a s =0) #d e f i n e M a t e r i a l
with d e f a u l t v a l u e s
5
6 O. m a t e r i a l s . append (p_mat)
7 O. m a t e r i a l s . append ( wall_mat )
8 O. m a t e r i a l s . append ( wallmat_b )

Create the cubic box:

1 # c r e a t e t h e box and add i t t o O. b o d i e s

40
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

2 O. b o d i e s . append ( u t i l s . w a l l (−R, a x i s =0, s e n s e =1, m a t e r i a l = wall_mat

) )#l e f t w a l l a l o n g x a x i s
3 O. b o d i e s . append ( u t i l s . w a l l ( 2 . 0 ∗R∗num_x−R, a x i s =0, s e n s e =−1,
m a t e r i a l = wall_mat ) )#r i g h t w a l l a l o n g x a x i s
4 O. b o d i e s . append ( u t i l s . w a l l (−R, a x i s =1, s e n s e =1, m a t e r i a l = wall_mat
) )#f r o n t w a l l a l o n g y a x i s
5 O. b o d i e s . append ( u t i l s . w a l l ( 2 . 0 ∗R∗num_y−R, a x i s =1, s e n s e =−1,
m a t e r i a l = wall_mat ) )#back w a l l a l o n g y a x i s
6 O. b o d i e s . append ( u t i l s . w a l l (−R, a x i s =2, s e n s e =1, m a t e r i a l =wallmat_b
) )#bottom w a l l a l o n g z a x i s
Define engines and set a fixed time step:
1 newton=N e w t o n I n t e g r a t o r ( damping = 0 . 3 , g r a v i t y = ( 0 . , 0 . , − 9 . 8 ) , l a b e l=
" newton " , i s S u p e r q u a d r i c s =4)
2
3 O. e n g i n e s =[
4 ForceResetter () ,
5 I n s e r t i o n S o r t C o l l i d e r ( [ Bo1_GJKParticle_Aabb ( ) , Bo1_Wall_Aabb ( )
] , v e r l e t D i s t =0.2∗0.01) ,
6 InteractionLoop (
7 [ Ig2_Wall_GJKParticle_GJKParticleGeom ( ) ,
Ig2_GJKParticle_GJKParticle_GJKParticleGeom ( ) ] ,
8 [ Ip2_GJKParticleMat_GJKParticleMat_GJKParticlePhys ( ) ] , #
c o l l i s i o n " physics "
9 [ GJKParticleLaw ( ) ] # c o n t a c t law −− apply f o r c e s
10 ),
11 newton
12 ]
13 O. dt=1e−5
Polyhedral particles are generated as follows:
1 #g e n e r a t e a sample
2 d e f GenSample ( r , pos ) :
3 f o r p i n pos :
4 body = GJKPolyhedron ( [ ] , [ 1 . e −2 ,1. e −2 ,1. e − 2 ] , 0 . 0 5 ∗ 1 e −2,p_mat ,
False )
5 body . s t a t e . pos=p
6 O. b o d i e s . append ( body )
7 O. b o d i e s [ − 1 ] . shape . c o l o r =(rand . random ( ) , rand . random ( ) , rand .
random ( ) )
8
9 # create a l a t t i c e grid
10 pos = G r i d I n i t i a l (R, num_x=num_x, num_y=num_y, num_z=num_z) # g e t
p o s i t i o n s o f a l l nodes
11 # c r e a t e p a r t i c l e s a t each nodes
12 GenSample (R, pos ) #
For cones, we change the GenSample functions as follows:
1 #g e n e r a t e a sample

41
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

Figure 3.11: Final packing of polyhedrons without random initial orientations

(the box not shown).

Figure 3.12: Final packing of cones without random initial orientations (the
box not shown).

2 d e f GenSample ( r , pos ) :
3 f o r p i n pos :
4 body = GJKCone ( 0 . 0 0 5 , 0 . 0 1 , 0 . 0 5 ∗ 1 e −2,p_mat , F a l s e )
5 body . s t a t e . pos=p
6 O. b o d i e s . append ( body )
7 O. b o d i e s [ − 1 ] . shape . c o l o r =(rand . random ( ) , rand . random ( ) , rand .

42
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

random ( ) )

Fig. 3.12 shows final packing of cones without random initial orientations
(i.e., the argument rotate is False). It can be seen that particles still align in
the lattice grid due to no disturbance in the tangential direction.
Set the argument rotate to True as follows:
1 body = GJKCone ( 0 . 0 0 5 , 0 . 0 1 , 0 . 0 5 ∗ 1 e −2,p_mat , True )

Figure 3.13: Final packing of cones with random initial orientations (the box
not shown).

Rerun the simulation, and the initial particles have random orientations.
After free falling under gravity, it can be seen that the lattice alignment
corrupts, referring to Fig. 3.13
For cylindrical particles, the function GenSample is replaced as follows:
1 d e f GenSample ( r , pos ) :
2 f o r p i n pos :
3 body = GJKCylinder ( 0 . 0 0 5 , 0 . 0 1 , 0 . 0 5 ∗ 1 e −2,p_mat , F a l s e )
4 body . s t a t e . pos=p
5 O. b o d i e s . append ( body )
6 O. b o d i e s [ − 1 ] . shape . c o l o r =(rand . random ( ) , rand . random ( ) , rand .
random ( ) )

Rerun the simulation, we obtain final packing of cylindrical particles as

shown in Fig. 3.14, where it is evident that all particles stay in a lattice grid

43
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

Figure 3.14: Final packing of cylinders without random initial orientations

(the box not shown).

Figure 3.15: Final packing of cylinders with random initial orientations (the
box not shown).

after free falling. Again, we reset the argument rotate to True so that all
initial particles will have random orientations.
1 body = GJKCylinder ( 0 . 0 0 5 , 0 . 0 1 , 0 . 0 5 ∗ 1 e −2,p_mat , True )

44
A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

Rerun the simulation again, the column of cylindrical particles collapse

into the cubic box. The final state is visualized in Fig. 3.15.

Figure 3.16: Final packing of cuboids without random initial orientations

(the box not shown).

Figure 3.17: Final packing of cuboids with random initial orientations (the
box not shown).

45
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

We then continue two more simulations of cuboids. Change the function

GenSample as follows:
1 #g e n e r a t e a sample
2 d e f GenSample ( r , pos ) :
3 f o r p i n pos :
4 body = GJKCuboid ( [ 0 . 0 0 5 , 0 . 0 0 5 , 0 . 0 0 5 ] , 0 . 0 5 ∗ 1 e −2,p_mat , F a l s e )
5 body . s t a t e . pos=p
6 O. b o d i e s . append ( body )
7 O. b o d i e s [ − 1 ] . shape . c o l o r =(rand . random ( ) , rand . random ( ) , rand .
random ( ) )

We obtain the final packing of cuboids without initial random orientations

as shown in Fig. 3.16. Again, the alignment of particles remains a lattice
form. Rerun the simulation with the argument rotate setting to True as
follows:
1 body = GJKCuboid ( [ 0 . 0 0 5 , 0 . 0 0 5 , 0 . 0 0 5 ] , 0 . 0 5 ∗ 1 e −2,p_mat , True )

The column of particles collapses but not too much due to small gaps between
particles.

3.5 Example 4: packing of poly-superellipsoids

Poly-superellipsoid has a capability of capturing major features (e.g., flatness,
asymmetry, and angularity) of realistic particles such as sands in nature. A
poly-superellipsoid is constructed by assembling eight pieces of superellipoids,
governed by the following surface function at a local Cartesian coordinate
system:
 x  2  y  2  1  z  2
  1 +   1 2 +   2 = 1 (3.11)
rx ry rz
with

rx = rx+ if x ≥ 0 else rx− (3.12a)

ry = ry+ if y ≥ 0 else ry− (3.12b)
rz = rz+ if z ≥ 0 else rz− (3.12c)

where rx+ , ry+ , rz+ and rx− , ry− , rz− are the principal elongation along positive and
negative directions of x, y, z axes, respectively; 1 , 2 control the squareness or
blockiness of particle surface, and their possible values are in (0, 2) for convex
shapes. Fig. 3.18 intuitively shows how 1 and 2 affect particle surface.
The module _superquadrics_utils provides two functions to generate a poly-
superellipsoid (see Sec. 5.1.2):

46
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

• NewPolySuperellipsoid(eps, rxyz, mat, rotate, isSphere, inertiaScale =

1.0)

• NewPolySuperellipsoid_rot(eps, rxyz, mat, qw , qx , qy , qz , isSphere, in-

ertiaScale = 1.0)

(a) (b) (c)

Figure 3.18: Poly-superellipsoids with rx+ = 1.0, rx− = 0.5, ry+ = 0.8, ry− =
0.9, rz+ = 0.4, rz− = 0.6 and (a) 1 = 0.4, 2 = 1.5, (b) 1 = 2 = 1.0, (c)
1 = 2 = 1.5.

Here we give an exemplified script for packing of poly-superellipsoids as

follows:
1 #
2 from sudodem . _ s u p e r q u a d r i c s _ u t i l s import ∗
3 import math
4 import random
5 import numpy a s np
6
7 a =1.1734 e−2
8 i s S p h e r e=F a l s e
9 ap =0.4
10 e p s =0.5
11 num_s=200
12 num_t=8000
13 t r i a l s =0
14
15
16 w a l l b o x s i z e =[10 e − 1 , 0 , 0 ] #s i z e o f t h e l a y e r f o r t h e p a r t i c l e s
generation
17 b o x s i z e=np . a r r a y ( [ 7 e −1 ,7 e −1 ,20 e −1])
18
19 d e f down ( ) :
20 f o r b i n O. b o d i e s :
21 i f ( i s i n s t a n c e ( b . shape , P o l y S u p e r e l l i p s o i d ) ) :
22 i f ( l e n ( b . i n t r s ( ) )==0) :
23 b . s t a t e . v e l =[0 ,0 , −2]
24 b . s t a t e . angVel = [ 0 , 0 , 0 ]
25

47
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

26 def gettop () :
27 zmax=0
28 f o r b i n O. b o d i e s :
29 i f ( i s i n s t a n c e ( b . shape , P o l y S u p e r e l l i p s o i d ) ) :
30 i f ( b . s t a t e . pos [2] >=zmax ) :
31 zmax=b . s t a t e . pos [ 2 ]
32 r e t u r n zmax
33
34
35 d e f G e n p a r t i c l e s ( a , eps , ap , mat , b o x s i z e , num_s , num_t , bottom ) :
36 global t r i a l s
37 gap=max( a , a∗ap )
38 #p r i n t ( gap )
39 num=0
40 coor =[]
41 width=( b o x s i z e [ 0 ] − 2 . ∗ gap ) / 2 .
42 l e n g t h =( b o x s i z e [ 1 ] − 2 . ∗ gap ) / 2 .
43 h e i g h t =( b o x s i z e [ 2 ] − 2 . ∗ gap )
44 #i t e r a t i o n =0
45
46 w h i l e num<num_s and t r i a l s <num_t :
47 isOK=True
48 pos = [ 0 ] ∗ 3
49 pos [ 0 ] = random . uniform (−width , width )
50 pos [ 1 ] = random . uniform ( gap−l e n g t h , l e n g t h −gap )
51 pos [ 2 ] = random . uniform ( bottom+gap , bottom+gap+h e i g h t )
52
53 f o r i in range (0 , len ( coor ) ) :
54 d i s t a n c e=sum ( [ ( ( c o o r [ i ] [ j ]− pos [ j ] ) ∗ ∗ 2 . ) f o r j i n
range ( 0 , 3 ) ] )
55 i f ( d i s t a n c e <(4.∗ gap ∗ gap ) ) :
56
57 isOK=F a l s e
58 break
59
60 i f ( isOK==True ) :
61 c o o r . append ( pos )
62 num+=1
63 t r i a l s +=1
64 #p r i n t (num)
65 return coor
66
67
68 d e f Addlayer ( a , eps , ap , mat , b o x s i z e , num_s , num_t) :
69 bottom=g e t t o p ( )
70 c o o r=G e n p a r t i c l e s ( a , eps , ap , mat , b o x s i z e , num_s , num_t , bottom )
71 f o r b in coor :
72 #xyz = np . a r r a y ( [ 1 . 0 , 0 . 5 , 0 . 8 , 1 . 5 , 1 . 2 , 0 . 4 ] ) ∗ 0 . 1
73 xyz = np . a r r a y ( [ 1 . 0 , 0 . 5 , 0 . 8 , 1 . 5 , 0 . 2 , 0 . 4 ] ) ∗ 0 . 1

48
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

74 bb=N e w P o l y S u p e r e l l i p s o i d ( [ 1 . 0 , 1 . 0 ] , xyz , mat , True , i s S p h e r e )#

75 bb . s t a t e . pos =[b [ 0 ] , b [ 1 ] , b [ 2 ] ]
76 bb . s t a t e . v e l =[0 ,0 , −1]
77 O. b o d i e s . append ( bb )
78 down ( )
79
80
81
82 p_mat = P o l y S u p e r e l l i p s o i d M a t ( l a b e l="mat1" ,Kn=1e5 , Ks=7e4 ,
f r i c t i o n A n g l e=math . atan ( 0 . 3 ) , d e n s i t y =2650 , betan =0, b e t a s =0) #
d e f i n e M a t e r i a l with d e f a u l t v a l u e s
83 wall_mat = P o l y S u p e r e l l i p s o i d M a t ( l a b e l=" wallmat " ,Kn=1e6 , Ks=7e5 ,
f r i c t i o n A n g l e =0.0 , betan =0, b e t a s =0) #d e f i n e M a t e r i a l with
default values
84 wallmat_b = P o l y S u p e r e l l i p s o i d M a t ( l a b e l=" wallmat " ,Kn=1e6 , Ks=7e5 ,
f r i c t i o n A n g l e=math . atan ( 1 ) , betan =0, b e t a s =0) #d e f i n e M a t e r i a l
with d e f a u l t v a l u e s
85
86 O. m a t e r i a l s . append (p_mat)
87 O. m a t e r i a l s . append ( wall_mat )
88 O. m a t e r i a l s . append ( wallmat_b )
89

90
91 O. b o d i e s . append ( u t i l s . w a l l ( −0.5∗ w a l l b o x s i z e [ 0 ] , a x i s =0, s e n s e =1,
m a t e r i a l = wall_mat ) )#l e f t x
92 O. b o d i e s . append ( u t i l s . w a l l ( 0 . 5 ∗ w a l l b o x s i z e [ 0 ] , a x i s =0, s e n s e =−1,
m a t e r i a l = wall_mat ) )#r i g h t x
93 O. b o d i e s . append ( u t i l s . w a l l ( −0.5∗ b o x s i z e [ 1 ] , a x i s =1, s e n s e =1,
m a t e r i a l = wall_mat ) )#f r o n t y
94 O. b o d i e s . append ( u t i l s . w a l l ( 0 . 5 ∗ b o x s i z e [ 1 ] , a x i s =1, s e n s e =−1,
m a t e r i a l = wall_mat ) )#back y
95
96 O. b o d i e s . append ( u t i l s . w a l l ( 0 , a x i s =2, s e n s e =1, m a t e r i a l =wallmat_b )
)#bottom z
97 #O. b o d i e s . append ( u t i l s . w a l l ( b o x s i z e [ 0 ] , a x i s =2, s e n s e =−1, m a t e r i a l
= wallmat_b ) )#top z
98
99 newton=N e w t o n I n t e g r a t o r ( damping = 0 . 3 , g r a v i t y = ( 0 . , 0 . , − 9 . 8 ) , l a b e l=
" newton " , i s S u p e r q u a d r i c s =2)#i s S u p e r q u a d r i c s =2 f o r Poly−
superellipsoids
100
101 O. e n g i n e s =[
102 ForceResetter () ,
103 I n s e r t i o n S o r t C o l l i d e r ( [ Bo1_PolySuperellipsoid_Aabb ( ) ,
Bo1_Wall_Aabb ( ) ] , v e r l e t D i s t = 0 . 2 ∗ 0 . 1 ) ,
104 InteractionLoop (
105 [ Ig2_Wall_PolySuperellipsoid_PolySuperellipsoidGeom () ,
Ig2_PolySuperellipsoid_PolySuperellipsoid_PolySuperellipsoidGeom
() ] ,

49
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

106 [
Ip2_PolySuperellipsoidMat_PolySuperellipsoidMat_PolySuperellipsoidPhys
() ] , # c o l l i s i o n " physics "
107 [ PolySuperellipsoidLaw () ] # c o n t a c t law −− apply f o r c e s
108 ),
109 newton
110 ]
111 O. dt=5e−5
112 #adding p a r t i c l e s
113 Addlayer ( a , eps , ap , p_mat , b o x s i z e , num_s , num_t)#Note : t h e s e
p a r t i c l e s may i n t e r s e c t with each o t h e r when we g e n e r a t e them
.
114
115 #s e t t i n g t h e D i s p l a y p r o p e r t i e s . You can go t o t h e D i s p l a y p a n e l
t o s e t them d i r e c t l y .
116 sudodem . qt . G l 1 _ P o l y S u p e r e l l i p s o i d . w i r e=F a l s e
117 sudodem . qt . G l 1 _ P o l y S u p e r e l l i p s o i d . s l i c e s =10
118 sudodem . qt . Gl1_Wall . d i v=0 #h i d e t h e w a l l s

After running the script above, the user is expected to see a packing as shown
in Fig. 3.19.

3.6 Example 5: packing of superellipses

The surface function of a superellipse in the local Cartesian coordinates can
be defined as  x 2  y 2
  +   = 1 (3.13)
rx ry
where rx and ry are referred to as the semi-major axis lengths in the direc-
tion of x, and y axies, respectively; and  ∈ (0, 2) is the shape parameters
determining the sharpness of particle edges or squareness of particle surface.
Here we simulate packing of superellipses under gravity with the following
script:
1 ##################################################
2 from sudodem import u t i l s
3 from sudodem . _ s u p e r e l l i p s e _ u t i l s import ∗
4 import numpy a s np
5 import math
6 import random
7 random . s e e d ( 1 1 2 4 5 6 4 2 )
8 ####################################################
9 #########c o n s t a n t s
10 ####################################################
11 isSphere = False
12
13 r = 0 . 5 #p a r t i c l e s i z e

50
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

Figure 3.19: Packing of poly-superellipsoids.

14 num_s=100
15 num_t=8000
16 t r i a l s =0
17 boxsize =[50.0 ,10.0]
18

19 #############################
20 mat =S u p e r e l l i p s e M a t ( l a b e l="mat1" ,Kn=1e8 , Ks=7e7 , f r i c t i o n A n g l e=
math . atan ( 0 . 2 2 5 ) , d e n s i t y =2650)
21 O. m a t e r i a l s . append ( mat )
22
23 O. b o d i e s . append ( u t i l s . w a l l ( −0.5∗ b o x s i z e [ 1 ] , a x i s =1, s e n s e =1,
m a t e r i a l = ' mat1 ' ) )#ground
24 O. b o d i e s . append ( u t i l s . w a l l ( −0.5∗ b o x s i z e [ 0 ] , a x i s =0, s e n s e =1,
m a t e r i a l = ' mat1 ' ) )#l e f t
25 O. b o d i e s . append ( u t i l s . w a l l ( 0 . 5 ∗ b o x s i z e [ 0 ] , a x i s =0, s e n s e =−1,
m a t e r i a l = ' mat1 ' ) )#r i g h t
26

27 ###
28 def gettop () :
29 ymax=0

51
 CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

30 f o r b i n O. b o d i e s :
31 i f ( i s i n s t a n c e ( b . shape , S u p e r e l l i p s e ) ) :
32 i f ( b . s t a t e . pos [1] >=ymax ) :
33 ymax=b . s t a t e . pos [ 1 ]
34 r e t u r n ymax
35
36
37 d e f G e n p a r t i c l e s ( r , mat , b o x s i z e , num_s , num_t , bottom ) :
38 global t r i a l s
39 gap=r
40 #p r i n t ( gap )
41 num=0
42 coor =[]
43 width=( b o x s i z e [ 0 ] − 2 . ∗ gap ) / 2 .
44 h e i g h t =( b o x s i z e [ 1 ] − 2 . ∗ gap )
45 #i t e r a t i o n =0
46
47 w h i l e num<num_s and t r i a l s <num_t :
48 isOK=True
49 pos = [ 0 ] ∗ 2
50 pos [ 0 ] = random . uniform (−width , width )
51 pos [ 1 ] = random . uniform ( bottom+gap , bottom+gap+h e i g h t )
52
53 f o r i in range (0 , len ( coor ) ) :
54 d i s t a n c e=sum ( [ ( ( c o o r [ i ] [ j ]− pos [ j ] ) ∗ ∗ 2 . ) f o r j i n
range ( 0 , 2 ) ] )
55 i f ( d i s t a n c e <(4.∗ gap ∗ gap ) ) :
56

57 isOK=F a l s e
58 break
59
60 i f ( isOK==True ) :
61 c o o r . append ( pos )
62 num+=1
63 t r i a l s +=1
64 #p r i n t (num)
65 return coor
66
67
68 d e f Addlayer ( r , mat , b o x s i z e , num_s , num_t) :
69 bottom=g e t t o p ( )
70 c o o r=G e n p a r t i c l e s ( r , mat , b o x s i z e , num_s , num_t , bottom )
71 f o r b in coor :
72 r r = r ∗random . uniform ( 0 . 5 , 1 . 0 )
73 bb = N e w S u p e r e l l i p s e ( 1 . 5 ∗ r , r , 1 . 0 , mat , True , i s S p h e r e )
74 bb . s t a t e . pos=b
75 bb . s t a t e . v e l =[0 , −1]
76 O. b o d i e s . append ( bb )
77 i f l e n (O. b o d i e s ) − 3 > 2 0 0 0 :

52
 A Quick Guide to SudoDEM CHAPTER 3. EXAMPLES

78 O. e n g i n e s [ − 1 ] . dead=True
79
80 O. dt = 1 e−3
81

82 newton=N e w t o n I n t e g r a t o r ( damping = 0 . 1 , g r a v i t y = ( 0 . , − 1 0 . 0 ) , l a b e l="

newton " , i s S u p e r e l l i p s e=True )
83
84 O. e n g i n e s =[
85 ForceResetter () ,
86 I n s e r t i o n S o r t C o l l i d e r ( [ Bo1_Superellipse_Aabb ( ) , Bo1_Wall_Aabb ( )
] , v e r l e t D i s t =0.1) ,
87 InteractionLoop (
88 [ Ig2_Wall_Superellipse_SuperellipseGeom ( ) ,
Ig2_Superellipse_Superellipse_SuperellipseGeom () ] ,
89 [ Ip2_SuperellipseMat_SuperellipseMat_SuperellipsePhys () ] , #
c o l l i s i o n " physics "
90 [ SuperellipseLaw () ] # c o n t a c t law −− apply f o r c e s
91 ),
92 newton ,
93 PyRunner ( command= ' Addlayer ( r , mat , b o x s i z e , num_s , num_t) ' ,
v i r t P e r i o d =0.1 , l a b e l= ' check ' , dead = F a l s e )
94 ]

After running the above script with sudodem2d, the user can get a packing
as shown in Fig. 3.20.

53
CHAPTER 3. EXAMPLES A Quick Guide to SudoDEM

Figure 3.20: Packing of superellipses under gravity.

54
 Chapter 4

Post-processing

4.1 Data
It is preferable to save a series of states during a simulation, i.e., periodically
saving data using a PyRunner:
1 def savedata () :
2 g l o b a l saving_num
3 saving_num += 1
4 #p r i n t " t e s t i n g , no data saved "
5 O. s a v e ( ' s h e a r '+s t r ( saving_num )+ ' . xml . bz2 ' )
6
7 O. e n g i n e s = O. e n g i n e s + [ PyRunner ( command = ' s a v e d a t a ( ) ' ,
i t e r P e r i o d = 1000 ) ] # i f we l o a d t h e c o n s o l i d a t i o n data ,
then we do not need t h i s l i n e b e c a u s e t h e PyRunner c a l l i n g '
s a v e d a t a ( ) ' has been d e f i n e d .
The Python module utilspost provides plentiful examples of post-processing,
including some basic functions of stress-force-fabric computation and other
pre-processing functions for visualization (e.g., writing vtk files of particles
and 3D histograms of fabric). Based on this module, users can copy and mod-
ify for a self purpose or just call the module functions in a post-processing
script, e.g.,
1 import sudodem . u t i l s p o s t a s u t i l s p o s t
2 #t h e f o l l o w i n g f u n c t i o n s w i l l output t h e a v e r a g e f r i c t i o n
m o b i l i z a t i o n i n d e x o f each s t a t e t o a s i n g l e f i l e .
3 for i in [ ' 1.0 ' , ' 1.2 ' , ' 1.4 ' , ' 1.6 ' ] :
4 u t i l s p o s t . c a l c _ a v g F r i c M o b i l I n d e x ( i , steps_num )
The Python module _superquadrics_utils also provides some auxiliary
functions, e.g.,
• outputParticles(filename) : filename, the file name for writing; output
info of each particle, i.e., shape parameters rx , ry , rz , 1 , 2 , position

55
 CHAPTER 4. POST-PROCESSING A Quick Guide to SudoDEM

x, y, z, and orientation q0,q1,q2,q3 (equivalent to qw , qx , qy , qz line by

line. respectively).

• outputParticlesbyIds(filename, ids) : filename, the file name for writ-

ing; ids, a list of ids for particles needed outputting; output info of each
particle in the list ids, i.e., shape parameters rx , ry , rz , 1 , 2 , position
x, y, z, and orientation q0,q1,q2,q3 (equivalent to qw , qx , qy , qz line by
line. respectively).

• outputWalls(filename) : filename, the file name for writing; output

positions of all walls.

• outputPOV(filename) : filename, the file name for writing; output a

POV-Ray file of particles for render using POV-Ray.

• outputVTK(filename,slices) : filename, the file name for writing; slices,

number of slices a particle needed slicing in longitude; output a vtk file
of particles for render using Paraview.

4.2 Scene Visualization

4.2.1 SudoDEM3D
For visualization, two free, open-source post-processing softwares, Paraview 1
and POV-Ray 2 , are strongly recommended.
It is preferable to reproduce a high-resolution figure of particle configu-
ration using post-processing softwares instead of a screenshot. As a demon-
stration, we load a state file named "[Link].bz2" and output a vtk file
or a pov file for post visualization.
1 O. l o a d ( " s h e a r e n d . xml . bz2 " )
2 import sudodem . _ s u p e r q u a d r i c s _ u t i l s a s s u p e r u t i l s
3 s u p e r u t i l s . outputVTK ( " s h e a r e n d . vtk " , 1 5 )

Open Paraview, then import the vtk file "[Link]" and click the button
"Apply", so that the configuration of particles is visualized as shown in Fig.
4.2. It is noteworthy that the function outputVTK processes a particle sur-
face as a combination of many (determined by the argument slices, i.e., 15
herein) polygons, which would yield a large file (around 180 Mb in the pre-
sented example). The users may find ways to reduce the vtk file for a faster
visualization in Paraview. One possible way is to show only the particles
1
[Link]
2
[Link]/

56
 A Quick Guide to SudoDEM CHAPTER 4. POST-PROCESSING

Figure 4.1: The main window of Paraview.

exposed in the view. Another possible way is to use the interface of built-in
Superquadric source in Paraview. Another alternative approach is rendering
the scene utilizing POV-Ray.
1 O. l o a d ( " s h e a r e n d . xml . bz2 " )
2 import sudodem . _ s u p e r q u a d r i c s _ u t i l s a s s u p e r u t i l s
3 s u p e r u t i l s . outputPOV ( " s h e a r e n d . pov " )

The users might revise the parameters of the camera in the POV file, i.e.,
"[Link]" here,
1 camera { // o r t h o g r a p h i c a n g l e 45
2 l o c a t i o n < 0 . 6 , 0 . 6 , 0 . 5 > ∗ 5 0 e−3
3 sky z
4 look_at <0 ,0 ,0.05 >∗10 e−3
5 r i g h t x∗ image_width / image_height
6 t r a n s l a t e < −0.04 ,0.01 , −0.1 >∗50 e−3
7 }

Light sources should be slightly adjusted as well.

1 // White background
2 background { rgb 1}
3 // Two l i g h t s with s l i g h t l y d i f f e r e n t c o l o r s
4 l i g h t _ s o u r c e { <4 ,8 ,5 >∗10 e−3 c o l o r rgb <1 ,1 ,1 >}
5 l i g h t _ s o u r c e {<12,−6>∗10e−3 c o l o r rgb <1 ,1 ,1 >}

In a terminal, run the following command to render the scene:

57
 CHAPTER 4. POST-PROCESSING A Quick Guide to SudoDEM

Figure 4.2: Configuration of particles rendered by Paraview.

1 povray [Link]

For a high-resolution render, additional arguments are appended to the com-

mand,
1 povray [Link] +A0.01 +W1600 +H1200

Image processing softwares, e.g., GIMP 3 , can be used to change image size or
make a gif animation, etc. Fig. 4.3 shows configuration of particles rendered
by POV-Ray and cropped using GIMP.
The module snapshot provides functions for visualizing poly-superellipsoids.
As an example, we load the state of final packing, and import the module
and call the outPov function as follows:
1 from sudodem import s n a p s h o t
2 s n a p s h o t . outPov ( " p a c k i n g p o l y s u p e r " , withbox=F a l s e )

we then get two files [Link] and [Link]. The first

file [Link] includes macros for the definition of poly-superellipsoid
and the particle information such as shape parameters, positions and orien-
tations, like this:
1 #macro myremain ( s i g n 1 , s i g n 2 , s i g n 3 )
2 clipped_by { p l a n e{− s i g n 1 ∗x , 0 } }
3 clipped_by { p l a n e{− s i g n 2 ∗y , 0 } }
4 clipped_by { p l a n e{− s i g n 3 ∗ z , 0 } }
3
[Link]

58
 A Quick Guide to SudoDEM CHAPTER 4. POST-PROCESSING

Figure 4.3: Configuration of particles rendered by POV-Ray.

5 #end
6 #macro randomColor ( r , g , b )
7 t e x t u r e { pigment { c o l o r rgb<r , g , b> t r a n s m i t para_trans } f i n i s h {
phong para_phong }}
8 #end
9 #macro o c t a n t S u p e r ( ep1 , ep2 , rx , ry , rz , s i g n 1 , s i g n 2 , s i g n 3 )
10 s u p e r e l l i p s o i d { <ep1 , ep2>
11 #i f ( s i n g l e C o l o r = f a l s e )
12 randomColor ( rand ( Random_r ) , rand (Random_g) , rand (Random_b) )
13 #end
14 myremain ( s i g n 1 , s i g n 2 , s i g n 3 ) s c a l e <rx , ry , rz >}
15 #end
16 #macro p o l y S u p e r E l l i p s o i d ( ep1 , ep2 , a1 , a2 , b1 , b2 , c1 , c2 , t1 , t2 , t3 , r1 ,
r2 , r3 , r , g , b )
17 union {
18 o c t a n t S u p e r ( ep1 , ep2 , a2 , b2 , c2 , −1 , −1 , −1)
19 o c t a n t S u p e r ( ep1 , ep2 , a1 , b2 , c2 ,1 , −1 , −1)
20 o c t a n t S u p e r ( ep1 , ep2 , a2 , b1 , c2 , −1 ,1 , −1)
21 o c t a n t S u p e r ( ep1 , ep2 , a1 , b1 , c2 , 1 , 1 , − 1 )
22 o c t a n t S u p e r ( ep1 , ep2 , a2 , b2 , c1 , −1 , −1 ,1)
23 o c t a n t S u p e r ( ep1 , ep2 , a1 , b2 , c1 , 1 , − 1 , 1 )
24 o c t a n t S u p e r ( ep1 , ep2 , a2 , b1 , c1 , − 1 , 1 , 1 )
25 o c t a n t S u p e r ( ep1 , ep2 , a1 , b1 , c1 , 1 , 1 , 1 )
26 r o t a t e <r1 , r2 , r3>
27 t r a n s l a t e <t1 , t2 , t3>
28 #i f ( s i n g l e C o l o r = t r u e )
29 randomColor ( r , g , b )
30 #end
31 }
32 #end
33 p o l y S u p e r E l l i p s o i d ( 1 . 0 0 0 0 0 0 0 e +00 ,1.0000000 e +00 ,1.0000000 e −01 ,

59
 CHAPTER 4. POST-PROCESSING A Quick Guide to SudoDEM

5 . 0 0 0 0 0 0 0 e −02 ,8.0000000 e −02 ,1.5000000 e −01 ,2.0000000 e −02 ,

4 . 0 0 0 0 0 0 0 e −02 , −3.2218276 e −01 ,1.8390308 e −01 ,3.1674861 e −01 ,
−167.727379 ,16.417442 , −104.719219 ,0.783099 ,0.394383 ,0.840188)
34 ...
35 ...

The users are not suggested to edit this file though. But the second file
[Link] includes some editable parameters as introduced above,
and looks like this:
1 // u s i n g t h e command : povray ∗ ∗ . pov +A0 . 0 1 +W1600 +H1200
2 #i n c l u d e " c o l o r s . i n c "
3 #d e c l a r e s i n g l e C o l o r = t r u e ;
4 #d e c l a r e Random_r = s e e d ( 1 4 3 2 ) ;
5 #d e c l a r e Random_g = s e e d ( 7 2 4 2 ) ;
6 #d e c l a r e Random_b = s e e d ( 9 9 1 2 ) ;
7 #d e c l a r e para_trans = 0 ;
8 #d e c l a r e para_phong = 1 . 0 ;
9 camera {
10 l o c a t i o n < 1.39244459003 , 0.950562440499 , 0.605226784362 >
11 sky z
12 r i g h t −x∗ image_width / image_height
13 look_at < 0 . 6 1 2 8 7 2 3 8 2 4 5 1 , 0 . 3 4 9 2 5 5 1 6 0 6 3 7 , 0 . 4 3 0 0 2 1 9 6 6 0 5 3 >
14 }
15 #i n c l u d e " p a c k i n g e l l i p s e s . i n c "
16 l i g h t _ s o u r c e { <4 ,8 ,5 >∗10 c o l o r rgb <1 ,1 ,1 >}
17 l i g h t _ s o u r c e {<12,−6>∗10 c o l o r rgb <1 ,1 ,1 >}
18 l i g h t _ s o u r c e { <0, 2 , 10> White }
19 background { rgb 1}
20 p l a n e { z , −5 pigment { c h e c k e r Green White }} // you may comment
out t h i s l i n e t o remove t h e p l a n e

Then render it in a fresh terminal by typing

1 povray [Link] +A0.01 +W1600 +H1200

you will get a nice figure with higher quality as Fig. 4.4.

4.2.2 SudoDEM2D
Module _superellipse_utils provides a function drawSVG (see Sec. 5.2.2) to
dump the particle profiles and contact forces chains to a SVG file.
1 from sudodem import _ s u p e r e l l i p s e _ u t i l s
2 O. l o a d ( " your data f i l e " )
3 drawSVG( " t e s t . svg " , d r a w _ f o r c e c h a i n=True ,
4 force_line_color =(1.0 ,0 ,0) ,
5 solo_color =(93.0/255 ,152.0/255 ,208.0/255) )

60
A Quick Guide to SudoDEM CHAPTER 4. POST-PROCESSING

Figure 4.4: Packing of poly-superellipsoids in Example 4 rendered by POV-

Ray.

The output file [Link] is editable by using a text editor and/or Inkscape
(recommended). Fig. 4.5 shows an exemplified snapshot of a packing of
elliptic particles with periodic boundary conditions. Refer to the keywords
of drawSVG for configuring the output figure.

61
CHAPTER 4. POST-PROCESSING A Quick Guide to SudoDEM

Figure 4.5: Particles and contact force chains dumped by drawSVG.

62
Chapter 5

Python Class Reference

5.1 SudoDEM3D
5.1.1 Basic classes
class State(inherits Serializable)
State of a body (spatial configuration, internal variables).

• angMom(=Vector3r::Zero())
Current angular momentum

• angVel(=Vector3r::Zero())
Current angular velocity

• blockedDOFs
Degress of freedom where linear/angular velocity will be always con-
stant (equal to zero, or to an user-defined value), regardless of applied
force/torque. String that may contain ‘xyzXYZ’ (translations and ro-
tations).

• dict() → dict
Return dictionary of attributes.

• inertia(=Vector3r::Zero())
Inertia of associated body, in local coordinate system.

• isDamped(=true)
Damping in Newtonintegrator can be deactivated for individual parti-
cles by setting this variable to FALSE. E.g. damping is inappropriate
for particles in free flight under gravity but it might still be applicable
to other particles in the same simulation.

63
CHAPTER 5. PYTHON CLASS REFERENCE
A Quick Guide to SudoDEM

• mass(=0)
Mass of this body
• ori
Current orientation.
• pos
Current position.
• refOri(=Quaternionr::Identity())
Reference orientation
• refPos(=Vector3r::Zero())
Reference position
• se3(=Se3r(Vector3r::Zero(), Quaternionr::Identity()))
Position and orientation as one object.

5.1.2 Module _superquadrics_utils

• NewSuperquadrics2(rx , ry , rz , 1 , 2 , mat, rotate, isSphere)
Generate a Superquadric with isSphere option.
• rx , ry , rz : float, semi-major axis lengths in x, y, and z axies.
• 1 , 2 : float, shape parameters.
• mat: Material, a material attached to the particle.
• rotate: bool, with a random orientation or not.
• isSphere: bool, particle is spherical or not.
• NewSuperquadrics_rot2(rx , ry , rz , 1 , 2 , mat, qw , qx , qy , qz , isSphere)
Generate a Super-ellipsoid with specified quaternion components: qw ,
qx , qy , qz with isSphere option.
• rx , ry , rz : float, semi-major axis lengths in x, y, and z axies.
• 1 , 2 : float, shape parameters.
• mat: Material, a material attached to the particle.
• isSphere: bool, particle is spherical or not.
• qw , qx , qy , qz : float, components of a quaternion (qw , qx , qy , qz )
representing the orientation of a particle.
• outputWalls(filename)
Output positions of all walls in the scene.

64
A Quick Guide to SudoDEM
CHAPTER 5. PYTHON CLASS REFERENCE

• filename: string, the name (with path) of the output file.

• getSurfArea(id, w, h)
Get the surface area of a super-ellipsoid by a given particle id with
resolution w and h (e.g., w=10,h=10). Larger w and h will yield more
accurate results.

• id: int, the id of the given particle.

• w, h: int, w×h slices the particle surface will be discretized into.

• outputParticles(filename)
Output information including id, shape parameters, position and ori-
entation of all particles into a file.

• filename: string, the name (with path) of the output file.

• outputParticlesbyIds(filename, ids)
Output particles by a list of ids, referring to the last function. [Only
for superellipsoids]

• filename: string, the name (with path) of the output file.

• ids: list of int, a list of particle ids, e.g., [10, 24, 22].

• NewPolySuperellipsoid(eps, rxyz, mat, rotate, isSphere, inertiaScale =

1.0)
Generate a PolySuperellipsoid.

• eps: Vector2, shape parameters, i.e., [1 , 2 ].

• rxyz: Vector6, semi-major axis lengths in x, y, and z axies, i.e.,[rx− ,
rx+ , ry− , ry+ , rz− , rz+ ].
• mat: Material, a material attached to the particle.
• rotate: bool, with a random orientation or not.
• isSphere: bool, particle is spherical or not.
• inertialScale: float, scale the moment of inertia only. Do not set
it if you are not sure the side effect.

• NewPolySuperellipsoid_rot(eps, rxyz, mat, qw , qx , qy , qz , isSphere, in-

ertiaScale = 1.0)
Generate a PolySuperellipsoid with specified quaternion components:
qw , qx , qy and qz .

• eps: Vector2, shape parameters, i.e., [1 , 2 ].

65
CHAPTER 5. PYTHON CLASS REFERENCE
A Quick Guide to SudoDEM

• rxyz: Vector6, semi-major axis lengths in x, y, and z axies, i.e.,[rx− ,

rx+ , ry− , ry+ , rz− , rz+ ].
• mat: Material, a material attached to the particle.
• isSphere: bool, particle is spherical or not.
• qw , qx , qy and qz : float, components of a quaternion (qw , qx , qy ,
qz ) representing the orientation of a particle.
• inertialScale: float, scale the moment of inertia only. Do not set
it if you are not sure the side effect.

• outputPOV(filename)
Output Superellipsoids (including spheres) into a POV file for post-
processing using POV-ray. [Only for superellipsoids and spheres]

• filename: string, the name (with path) of the output file.

• PolySuperellipsoidPOV(filename, ids = [], scale = 1.0)

Output PolySuperellipsoids into a POV file for post-processing using
POV-ray. This function only output the particles’ information to a
POV file like *.inc, and the user may need to add some other info such
as camera etc. See the module snapshot.

• filename: string, the name (with path) of the output file.

• ids: list of int, a list of particle ids. A NULL list will output all
particles by default.
• scale: float, scale the length unit.

• outputVTK(filename, slices)
Output Superellpisoids into a VTK file for post-processing using Par-
aview. [Only for superellipsoids]

• filename: string, the name (with path) of the output file.

• slices: int, into how many slices the particle surface will be dis-
cretized.

5.1.3 Module _gjkparticle_utils

• GJKSphere(radius,margin,mat)

• radius, float, radius of the sphere.

• margin, float, a small margin added to the surface for assisting
contact detection.

66
A Quick Guide to SudoDEM
CHAPTER 5. PYTHON CLASS REFERENCE

• mat, Material, a specified material attached to the particle.

• GJKPolyhedron(vertices,extent,margin,mat, rotate)

• vertices, a list of vertices, and each vertex is a list, e.g., [x, y, z]

for coordinates.
• extent, a list of three components to scale the x, y, z dimensions;
valid only if vertices has no elements (i.e., []); a randomly shaped
polyhedron will be generated based on Voronoi tessellation.
• margin, float, a small margin added to the surface for assisting
contact detection.
• mat, Material, a specified material attached to the particle.
• rotate, bool, with random orientation or not.

• GJKCone(radius, height, margin, mat, rotate)

• radius, float, base radius of the cone.

• height, float, height of the cone.
• margin, float, a small margin added to the surface for assisting
contact detection.
• mat, Material, a specified material attached to the particle.
• rotate, bool, with random orientation or not.

• GJKCylinder(radius, height, margin, mat, rotate)

• radius, float, radius of the cylinder.

• height, float, height of the cylinder.
• margin, float, a small margin added to the surface for assisting
contact detection.
• mat, Material, a specified material attached to the particle.
• rotate, bool, with random orientation or not.

• GJKCuboid(extent, margin,mat, rotate)

• extent, a list of three floats, extent of the cuboid in the x, y, z

dimensions.
• margin, float, a small margin added to the surface for assisting
contact detection.
• mat, Material, a specified material attached to the particle.

67
CHAPTER 5. PYTHON CLASS REFERENCE
A Quick Guide to SudoDEM

• rotate, bool, with random orientation or not.

Note: a margin is used to expand a particle for achieving a better compu-

tational efficiency, which however should be large enough for ensuring that
collision occurs in this small region but small enough for shape fidelity.

5.1.4 Module snapshot

• outputBoxPov(filename, x, y, z, r=0.001, wallMask=[1,0,1,0,0,1])
Output *.POV of a cubic box (x=[x_min, x_max], y=[y_min, y_max],
z=[z_min, z_max]) for post-processing in Pov-ray.

• filename: string, the name (with path) of the output file.

• x, y, z: Vector2, ranges of the cubic box along the three (x, y, z)
directions.
• r: float, thickness of the box edge.
• wallMask, list of 6 int, corresponding to the six walls at (x-, x+,
y-, y+, z-, z+). 1 to show and 0 to hide the wall.

• outPov(fname,transmit=0,phong=1.0,singleColor=True,ids=[],scale=1.0,
withbox = True)

• filename, string, the name (with path) of the output file.

• transmit, float, parameter for transparency of particles, ranging
between 0 and 1.
• phong, float, parameter for highlighting the surface, ranging be-
tween 0 and 1.
• singleColor, bool, whether to use a single color for all particles. If
not, particles’ colors would be the values defined in Shape.
• ids, list of int, a list of particle ids that will be output. If no item
are given, then all particles will be output by default.
• scale, float, scale the unit.

The user is referred to POV-Ray’s Doc for more details on the two
parameters (transmit and phong) and others that will be written to
the output file by this function.

68
A Quick Guide to SudoDEM
CHAPTER 5. PYTHON CLASS REFERENCE

5.2 SudoDEM2D
5.2.1 Basic classes
Class State
State of a body (spatial configuration, internal variables).

• angMom(=Vector2r::Zero())
Current angular momentum

• angVel(=Vector2r::Zero())
Current angular velocity

• blockedDOFs
Degress of freedom where linear/angular velocity will be always con-
stant (equal to zero, or to an user-defined value), regardless of applied
force/torque. String that may contain ‘xyZ’ (translations and rota-
tions).

• dict() → dict
Return dictionary of attributes.

• inertia(=Vector2r::Zero())
Inertia of associated body, in local coordinate system.

• isDamped(=true)
Damping in Newtonintegrator can be deactivated for individual parti-
cles by setting this variable to FALSE. E.g. damping is inappropriate
for particles in free flight under gravity but it might still be applicable
to other particles in the same simulation.

• mass(=0)
Mass of this body

• ori
Current orientation. The orientation and rotation of a particle is rep-
resented by a Rotation2d object (see Rotation2D in the Eigen library),
which has the following members and functions.
Class Rotation2d

– angle: float, rotation represented by an angle in rad.

– toRotationMatrix(): Matrix2r, returns an equivalent 2x2 rotation
matrix.

69
CHAPTER 5. PYTHON CLASS REFERENCE
A Quick Guide to SudoDEM

– Identity(): Roation2d, returns an identity rotation.

– smallestAngle(): float, returns the rotation angle in [−π, π]
– inverse(): Roation2d, returns the inverse rotation
– smallestPositiveAngle(): float, returns the rotation angle in [0, 2π]

• pos
Current position.

• refOri(=Rotation2D::Identity())
Reference orientation

• refPos(=Vector2r::Zero())
Reference position

• se2(=Se2r(Vector2r::Zero(), Rotation2d::Identity()))
Position and orientation as one object.

• vel(=Vector2r::Zero()) Current linear velocity.

5.2.2 Module _superellipse_utils

• NewSuperellipse(rx, ry, epsilon, material, rotate, isSphere, z_dim =1.0)
Create a Superellipse

– rx, ry: float, semi-axis length of the particle

– epsilon: float, shape parameter of a superellipse
– material: Material instance
– rotate: bool, rotate the particle with random orientation?
– isSphere: bool, is the superellipse a disk? A disk will speed up
the computation.
– z_dim: float, the virtual length at the z dimension, and set to 1.0
by default.

• NewSuperellipse_rot(x, y, epsilon, material, miniAngle, isSphere, z_dim

= 1.0)
Create a superellipse with certain orientation

– rx, ry: float, semi-axis length of the particle

– epsilon: float, shape parameter of a superellipse

70
A Quick Guide to SudoDEM
CHAPTER 5. PYTHON CLASS REFERENCE

– material: Material instance

– miniAngle: orientation of the particle specified by the angle be-
tween its rx-axis and the global x axis
– isSphere: bool, is the superellipse a disk? A disk will speed up
the computation.
– z_dim: float, the virtual length at the z dimension, and set to 1.0
by default.
• outputParticles(filename)
output particle info(rx,ry,eps,x,y,rotation angle) to a text file
• drawSVG(filename, width_cm = 10.0, slice = 20, line_width = 0.1,
fill_opacity = 0.5, draw_filled = true, draw_lines = true, color_po =
false, po_color = Vector3r(1.0,0.0,0.0), solo_color = Vector3r(0,0,0),
force_line_width = 0.001, force_fill_opacity = 0, force_line_color =
Vector3r(1.0,1.0,1.0), draw_forcechain = false)
Output particle profiles and contact force chains into a SVG file.
– filename: string, the name of the output file
– width_cm: float, width (in centimeters) in the SVG header
– slice: int, to how many slices a Superellipse will be discretized
– line_width: float, the line width of the edge of a Superellipse
– fill_opacity: float, the opacity of the fill color in a Superellipse
– draw_filled: bool, wheter to fill a Superellipse with a color
– draw_lines: bool, whether to draw the out profile of a Superellipse
– color_po: bool, whether to use a fill color to identify the orienta-
tion of a Superellipse
– po_color: Vector3r, the fill color to identify the orientation of a
Superellipse
– solo_color: Vector3r, a solor color to fill a Superellipse. If its norm
is zero, then use the color defined by po_color or shape→color
– force_line_width: float, the line width of the force chain with the
average normal contact force
– force_fill_opacity: float, the opacity of force chains
– force_line_color: Vector3r, the color of the force chain
– draw_forcechain: bool, whether to draw the contact force chain
that will be superposed with the particles

71
CHAPTER 5. PYTHON CLASS REFERENCE
A Quick Guide to SudoDEM

5.2.3 Module _utils

• unbalancedForce(useMaxForce=false)
Compute the ratio of mean (or maximum, if *useMaxForce*) summary
force on bodies and mean force magnitude on interactions. For perfectly
static equilibrium, summary force on all bodies is zero (since forces
from interactions cancel out and induce no acceleration of particles);
this ratio will tend to zero as simulation stabilizes, though zero is never
reached because of finite precision computation. Sufficiently small value
can be e.g. 1e-2 or smaller, depending on how much equilibrium it
should be.

• getParticleVolume2D()
Compute the total volume (area for 2D) of particles in the scene.

• getMeanCN()
Get the mean coordination number of a packing

• changeParticleSize2D(alpha)
expand a particle by a given coefficient

• getVoidRatio2D(cellArea=1)
Compute 2D void ratio. Keyword cellArea is effective only for aperioidc
cell

• getStress2D(z_dim=1)
Compute overall stress of periodic cell

• getFabricTensorCN2D()
Fabric tensor of contact normal

• getFabricTensorPO2D()
Fabric tensor of particle orientation along the rx axis

• getStressAndTangent2D(z_dim=1,symmetry=true)
Compute overall stress of periodic cell using the same equation as func-
tion getStress. In addition, the tangent operator is calculated using the
equation: Sijkl = V c (kn ni lj nk ll + kt ti lj tk ll ) float volume: same as in
1
P
function getStress bool symmetry: make the tensors symmetric.
return: macroscopic stress tensor and tangent operator as py::tuple

• getStressTangentThermal2D(z_dim=1, symmetry =true)

Compute overall stress of periodic cell using the same equation as func-
tion getStress. In addition, the tangent operator is calculated using

72
A Quick Guide to SudoDEM
CHAPTER 5. PYTHON CLASS REFERENCE

the equation: Sijkl = V1 c (kn ni lj nk ll + kt ti lj tk ll ) float volume: same

P
as in function getStress bool symmetry: make the tensors symmet-
ric. Finally, the thermal conductivity tensor is calculated based on the
formular in PFC. macroscopic stress tensor and tangent operator as
py::tuple

5.2.4 Module utils

• disk(center, radius, z_dim=1, dynamic=None, fixed=False, wire=False,
color=None, highlight=False, material=-1,mask=1)
Create a disk with given parameters; mass and inertia computed auto-
matically.

– center: Vector2, center

– radius: float, radius
– dynamic: float, deprecated, see "fixed"
– fixed: float, generate the body with all DOFs blocked?
– material: specify ‘[Link]‘; different types are accepted:
∗ int: [Link][material] will be used; as a special case,
if material==-1 and there is no shared materials defined,
[Link]() will be assigned to [Link][0]
∗ string: label of an existing material that will be used
∗ ‘Material‘ instance: this instance will be used
∗ callable: will be called without arguments; returned Material
value will be used (Material factory object, if you like)
– int mask: ‘[Link]‘ for the body
– wire: display as wire disk?
– highlight: highlight this body in the viewer?
– Vector2-or-None: body’s color, as normalized RGB; random color
will be assigned if “None“.

• wall(position, axis, sense=0, color=None, material=-1, mask=1)

Return ready-made wall body.

– position: float-or-Vector3 , center of the wall. If float, it is the

position along given axis, the other 2 components being zero
– axis: {0,1} , orientation of the wall normal (0,1) for x,y

73
CHAPTER 5. PYTHON CLASS REFERENCE
A Quick Guide to SudoDEM

– sense: {-1,0,1} , sense in which to interact (0: both, -1: negative,

+1: positive; see ‘Wall‘)

See ‘[Link]‘’s documentation for meaning of other parameters.

• fwall(vertex1, vertex2, dynamic=None, fixed=True, color=None, high-

light=False, noBound=False, material=-1, mask=1, chain=-1)
Create fwall with given parameters.

– vertex1: Vector2, coordinates of vertex1 in the global coordinate

system.
– vertex2: Vector2, coordinates of vertex2 in the global coordinate
system.
– noBound: bool, set ‘[Link]‘
– color: Vector3-or-None , color of the facet; random color will be
assigned if ‘None‘.

See ‘[Link]‘’s documentation for meaning of other parameters.

74
Acknowledgments

The work was partially supported by the Hong Kong Scholars Program
(2018), the National Natural Science Foundation of China (by Project No.
51679207, 51909095), Research Grants Council of Hong Kong (by GRF Project
No. 16205418, TBRS Project No. T22-603/15N and CRF Project No.
C6012-15G). The following open-source tools but not limited to are acknowl-
edged: Ubuntu, YADE, Python, Boost, Eigen, Voro++, Paraview, Pov-Ray,
GIMP, InkScape, Veusz and LaTeX.

75