Example 1: packing of super-ellipsoids

The surface function of a superellipsoid in the local Cartesian coordinates can be defined as $$\label{eqsurface} \Big(\big|\frac{x}{r_x}\big|^{\frac{2}{\epsilon_1}} + \big|\frac{y}{r_y}\big|^{\frac{2}{\epsilon_1}}\Big)^{\frac{\epsilon_1}{\epsilon_2}} + \big|\frac{z}{r_z}\big|^{\frac{2}{\epsilon_2}} = 1$$ where $r_x, r_y$ and $r_z$ are referred to as the semi-major axis lengths in the direction of x, y, and z axies, respectively; and $\epsilon_i (i=1,2)$ are the shape parameters determining the sharpness of particle edges or squareness of particle surface. Varying $\epsilon_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{reference-type=“ref” reference=“modsuperquadrics”}):

• NewSuperquadrics2($r_x$, $r_y$, $r_z$, $\epsilon_1$, $\epsilon_2$, mat, rotate, isSphere)

• NewSuperquadrics_rot2($r_x$, $r_y$, $r_z$, $\epsilon_1$, $\epsilon_2$, mat, $q_w$, $q_x$, $q_y$, $q_z$, isSphere)

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

###############################################
# granular packing of super-ellipsoids
# We position particles at a lattice grid, 
# then particles free fall into a cubic box.
###############################################
# import some modules
from sudodem import _superquadrics_utils

from sudodem._superquadrics_utils import *
import math
import random as rand
import numpy as np

########define some parameters#################
isSphere=False

num_x = 5
num_y = 5
num_z = 20
R = 0.1

#######define some auxiliary functions#########
               
#define a latice grid
#R: distance between two neighboring nodes
#num_x: number of nodes along x axis
#num_y: number of nodes along y axis
#num_z: number of nodes along z axis
#return a list of the postions of all nodes
def GridInitial(R,num_x=10,num_y=10,num_z=20):
   pos = list()
   for i in range(num_x):
      for j in range(num_y):
         for k in range(num_z):
            x = i*R*2.0
            y = j*R*2.0
            z = k*R*2.0
            pos.append([x,y,z])       
   return pos

#generate a sample
def GenSample(r,pos):
   for p in pos:
		epsilon1 = rand.uniform(0.7,1.6)
		epsilon2 = rand.uniform(0.7,1.6)
		a = r
		b = r*rand.uniform(0.4,0.9)#aspect ratio
		c = r*rand.uniform(0.4,0.9)
		
		body = NewSuperquadrics2(a,b,c,epsilon1,epsilon2,p_mat,True,isSphere)
		body.state.pos=p
		O.bodies.append(body)

#########setup a simulation####################
# material for particles
p_mat = SuperquadricsMat(label="mat1",Kn=1e5,Ks=7e4,frictionAngle=math.atan(0.3),density=2650,betan=0,betas=0)
# betan and betas are coefficients of viscous damping at contact, no viscous damping with 0 by default.
# material for the side walls
wall_mat = SuperquadricsMat(label="wallmat",Kn=1e6,Ks=7e5,frictionAngle=0.0,betan=0,betas=0)
# material for the bottom wall
wallmat_b = SuperquadricsMat(label="wallmat",Kn=1e6,Ks=7e5,frictionAngle=math.atan(1),betan=0,betas=0)
# add materials to O.materials
O.materials.append(p_mat)
O.materials.append(wall_mat)
O.materials.append(wallmat_b)

# create the box and add it to O.bodies
O.bodies.append(utils.wall(-R,axis=0,sense=1, material = wall_mat))#left wall along x axis
O.bodies.append(utils.wall(2.0*R*num_x-R,axis=0,sense=-1, material = wall_mat))#right wall along x axis
O.bodies.append(utils.wall(-R,axis=1,sense=1, material = wall_mat))#front wall along y axis
O.bodies.append(utils.wall(2.0*R*num_y-R,axis=1,sense=-1, material = wall_mat))#back wall along y axis
O.bodies.append(utils.wall(-R,axis=2,sense=1, material =wallmat_b))#bottom wall along z axis

# create a lattice grid
pos = GridInitial(R,num_x=num_x,num_y=num_y,num_z=num_z) # get positions of all nodes
# create particles at each nodes
GenSample(R,pos)

# create engines
# define a Newton engine
newton=NewtonIntegrator(damping = 0.1,gravity=(0.,0.,-9.8),label="newton",isSuperquadrics=1) # isSuperquadrics: 1 for superquadrics

O.engines=[
   ForceResetter(), InsertionSortCollider([Bo1_Superquadrics_Aabb(),Bo1_Wall_Aabb()],verletDist=0.2*0.1),
   InteractionLoop(
      [Ig2_Wall_Superquadrics_SuperquadricsGeom(),    Ig2_Superquadrics_Superquadrics_SuperquadricsGeom()], 
      [Ip2_SuperquadricsMat_SuperquadricsMat_SuperquadricsPhys()], # collision "physics"
      [SuperquadricsLaw()]   # contact law
   ),
   newton
   
]
O.dt=5e-5

With a GUI controller, we can click ‘show 3D’ button at the panel ‘Simulation’ 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{reference-type=“ref” reference=“figguidisplay”}.

Figs. 3.4{reference-type=“ref” reference=“figexp11”} - 3.6{reference-type=“ref” reference=“figexp13”} show simulations of packing of superellipsoids at different states for the presented example.

Several properties and methods for a shape Superquadrics are listed below:

• Properties:

  • $r_x$, $r_y$, $r_z$: float, semi-major axis lengths in x, y, and z axies.

  • $\epsilon_1$, $\epsilon_2$: float, shape parameters.

  • isSphere: bool, particle is spherical or not.

• Methods:

  • getVolume(): float, return particle’s volume.

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

volume = O.bodies[100].shape.getVolume()

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

O.bodies[100].dict()

The output is like follows:

{'bound': <Aabb instance at 0x56072ec83da0>,
 'chain': -1,
 'clumpId': -1,
 'flags': 3,
 'groupMask': 1,
 'id': 100,
 'iterBorn': 0,
 'material': <SuperquadricsMat instance at 0x56072e6ed5f0>,
 'shape': <Superquadrics instance at 0x56072ec83b60>,
 'state': <State instance at 0x56072ec839e0>,
 '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:

O.bodies[100].state.dict()

The output is as follows:

{'angMom': Vector3(0,0,0),
 'angVel': Vector3(0,0,0),
 'blockedDOFs': 0,
 'densityScaling': 1.0,
 'inertia': Vector3(0.0045389863901528615,0.007287422614611933,
                           0.0067053718092146865),
 'isDamped': True,
 'mass': 3.225715855242557,
 'refOri': Quaternion((1,0,0),0),
 'refPos': Vector3(0,0,0),
 'se3': (Vector3(0,0.8,3),
  Quaternion((0.3970010520699211,0.5339678220536823,
                     -0.7465042060609055),2.1754719537556495)),
 'vel': Vector3(0,0,0)}

For an object Shape, we list all properties which may vary for different child classes.

O.bodies[100].shape.dict()

The output is as follows:

{'color': Vector3(0.6407663308273844,0.07426042997942327,
                            0.05872324344642611),
 'eps': Vector2(1.3975848801318045,1.5376309088540607),
 'eps1': 1.3975848801318045,
 'eps2': 1.5376309088540607,
 'highlight': False,
 'isSphere': False,
 'rx': 0.1,
 'rxyz': Vector3(0.1,0.06469580193313924,0.07572423080016706),
 'ry': 0.06469580193313924,
 'rz': 0.07572423080016706,
 'wire': False}

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

for b in O.bodies: # loop all bodies in the simulation
    # we have to check if the present body is a superellipsoid or others, e.g., a wall
    if isinstance(b.shape, Superquadrics): # if it is a superellipsoid, then to do
        pid = b.id # get the id of particle b
        pos = b.state.pos # get the position of particle b
        print pid, pos