1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
|
import sys
## You may need to specify the MPM package path
## if it is not installed in PYTHONPATH
sys.path.append("/path/to/your/mpm/package")
## GUI settings
app.core.setAnalysisType(AnalysisType = "MPM")
## Import MPM python package
import mpm as task
## Initialization of the core: Analysis object
sim=app.core.sim
## Geometry of mesh
elemsize = 0.25
width = 10
height = 5
## Initialization of the grid: Grid object
grid = sim.grid
## Dimension of the grid
## Arguments: (xmin, xmax, ymin, ymax)
grid.setRange(0, width, 0, height)
## Element numbers
## Arguments: (numberx, numbery)
grid.setNumEle(40, 20)
## Material 1: Mohr-Coulomb material model
mat1 = task.Mc()
mat1.rho = 1
mat1.infinistate = False
## Material 2: Mohr-Coulomb material model
mat2 = task.Mc()
mat2.rho = 1300.
mat2.E = 5e5
mat2.nu = 0.4
mat2.c = 0.
mat2.phi = 22
mat2.psi = 0.01
mat2.tension = 0.01
mat2.infinistate = False
## Contact stiffness for the frictional base
mat2.Kn = 5e8 ## normal contact stiffness
mat2.Ks = 2.5e7 ## shear contact stiffness
mat2.friction = 0.5 ## friction coefficient
## Material array as member of Analysis object
sim.mats = [mat1, mat2]
## Assign the id for each material
sim.mats[-1].setIndex(len(sim.mats)-1)
sim.mats[-2].setIndex(len(sim.mats)-2)
## Dampings
pdamping = 1.0
gdamping = 1.0
pic = 1.0
## Timestep
dt = 1e-4
## Shape objects:
## GRAVITYSHAPE defines a gravity BC, all the material points
## inside the shape will bear a body force = [0, -9.8]
soil = task.Rectangle(2.*elemsize, 2.*elemsize+2, 2.*elemsize, 2.*elemsize + 4, task.GRAVITYSHAPE)
## SetBC for the GRAVITYSHAPE is:
## Arguments: (isxsetgravity, gravity_x, isysetgravity, gravity_y)
soil.setBc(False, 0, True, -9.8)
soil.setMatProperty(task.MCMAT, 1, 2)
soil.setDamping(pdamping, gdamping, pic)
left = task.Rectangle(elemsize, 2.*elemsize, 2.*elemsize, 2.*elemsize + 5, task.RIGIDSHAPE)
left.setBc(True, 0, False, 0)
left.setMatProperty(task.RIGIDMAT, 0, 1)
right = task.Rectangle(2.*elemsize+2, 3.*elemsize+2, 2.*elemsize, 2.*elemsize + 5, task.RIGIDSHAPE)
right.setBc(True, 0, False, 0)
right.setMatProperty(task.RIGIDMAT, 0, 1)
## DemPolygon DOES NOT hold any material points,
## but it will impose the contact boundary conditions.
## Force between material points and DemPolygon stands on their
## penetration and the relative veloicty,
## which is similar to DEM contact law
## Create a DemPoly by inputing all their vertex
## in anti-clockwise or clock-wise order
## Make sure the vertexes follow the same direction
## Any intersection between two lines formed by three adjacenting
## material points is not allowed
pts = [0, elemsize, width, elemsize, width, 2.*elemsize, 0., 2.*elemsize]
## Constructor of DemPolygon:
## Arguments:
## (timestep, search_range, update_iteration, normal_contact_stiffness,
## shear_contact_stiffness, friction_coefficient)
foot = task.DemPolygon(dt, 10, 100, 5e6, 5e6, 0.5);
## Set the vertex array
foot.setPts(pts)
## UsePolygonContact determines if material points is regarded
## as sphere or polygon
## True: polygon False: sphere
foot.usePolygonContact(False)
## Set the DemPolygon velocity boundary condition
foot.setBc(True, 0, True, 0)
foot.setMatProperty(task.RIGIDMAT, 0, 1)
## Summary: left and right shape to provide the fixed velocity
## or displacement constraint, the bot works like
## a friction base and only apply a vertical displacement
## constraint. And soil is the samples that will suffer
## from a collapse
## Shapes array as a member of Analysis object
sim.shapes = [soil, left, foot, right]
## Analysis settings
## CPU cores (at least 1)
## SetNumThreads must be called explicitly,
## otherwise undeifned behaviour will appear
sim.setNumThreads(1)
## Set the total time elapsed and timestep
## Arguments: (tottal_time, timestep)
sim.setTimeDt(1, dt)
## Simulation resutl archive
## When to start the archive (sec)
sim.writer.setStartTime(0)
## Time to stop the archive (sec)
sim.writer.setEndTime(10)
## Frequency for the file archive (sec)
sim.writer.setIntervalTime(0.1)
## Where to store the archived datas
## (directory_name, file_name)
sim.writer.setDir("./Results/", 'granular_collapse')
## MPM interpolation method,
## options: LINEAR, GIMP, LINEARCPDI, BSPLINECPDI
sim.setMethod(task.GIMP)
## Simulation task, options: MANAL, TANAL, PHASEANAL
## MANAL: mechanical simulation
## TANAL: thermal simulation
## PHASEANAL: phase transition simulation
sim.simTypeMasks = task.MANAL
def gravityHookFunc(steps = 100):
## Get the current iteration steps by sim.iter()
partg = -9.8*sim.iter()/float(steps)
## Try to print information to app, use the app.print
app.print('Current g:', partg)
## Set the accumulated increasing gravity
soil.setBc(False, 0, True, partg)
## The PyHook allows the user to run the customized task
## after each MPM cycle
pyhook1 = task.PyHook()
## State the hooked function by string
pyhook1.command = 'gravityHookFunc()'
## Set the hook task interval and total_runnings
## Arguments: (interval, total_runnings)
## The elapsed iteration = interval * total_runnings
pyhook1.reset(iterPeriod = 1, totalRuns = 100)
## Dead=True indicates the hook stopped, False means activated
pyhook1.dead = False
def removeRightWall():
## Note that the rigidshape can be removed at any time
## by setting the bool flag to false
right.setBc(False, 0, False, 0)
## Set all the damping to 0 to capture a dynamic simulation
soil.setDamping(0, 0, 0)
## Before the calculation,
## run PreAnalysis to initialize the analysis core
sim.PreAnalysis()
## Run specific steps
## If GUI is activated,
## app.run(1)
## else:
## sim.run(1)
|