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

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

## Initialization of the grid: Grid object
grid = sim.grid

## Dimension of the grid
## Arguments: (xmin, xmax, ymin, ymax)
grid.setRange(0, 2.6, 0, 2.6)

## Element numbers
## Arguments: (numberx, numbery)
grid.setNumEle(26, 26)

## Print necessary information of grid
grid.printInfo()

## Element size
elemsize = 0.1

## Dampings
pdamping = 0.1   ## particle damping
gdamping = 0.    ## node damping
pic = 0.1        ## PIC fraction for mixed format 
                 ## between FLIP and PIC
                 ## 0 for pure FLIP, 1 for pure PIC

## Target stress for the consolidation stage
q = -1e5

## Material 1: Mohr-Coulomb material model 
mat1 = task.Mc()
mat1.rho = 1                  ## density (kg/m3)
mat1.infinistate = False      ## True:  finite strain theory adopted
                              ## False: infinitesimal starin theory adopted
                            
## Material 2: Mohr-Coulomb material model
phi = 20         ## friction angle (deg)
c = 0.1          ## cohesive (Pa)
psi = 10         ## dilation angle (deg)

mat2 = task.Mc()
mat2.rho = 2650               ## density (kg/m3)
mat2.E = 1e5                  ## Young's modulus (Pa)
mat2.nu = 0.2                 ## Poisson's ratio
mat2.c = c                    ## cohesive (Pa)
mat2.phi = phi                ## friction angle (deg)
mat2.psi = psi                ## dilation angle (deg)
mat2.tension = 0.2            ## tension cut off (Pa)
mat2.infinistate = False      ## True:  finite strain theory adopted
                              ## False: infinitesimal starin theory adopted

## Material array as member of Analysis object
sim.mats = [mat1, mat2]

## Ids for each material
sim.mats[-1].setIndex(len(sim.mats)-1)
sim.mats[-2].setIndex(len(sim.mats)-2)

## Shape objects:
## Derived classes of shape are: Rectangle, Disk, Ring, Polygon and DemPolygon

## Rectangle
## Arguments: (xmin, xmax, ymin, ymax, boundary_shape)
left = task.Rectangle(0.8, 0.85, 0.3, 2.3, task.MPTTRACSHAPE)

## MPTTRACSHAPE means that tranction boundary condition 
## will be applied to material points inside the left object
## Arguments: (faceid, direction, signed_magnitude)

## Face id:
## --3--
## |   |
## 4   2
## --1--

## Direction options: XDIR, YDIR, NORMALDIR, TAGENTDIR
## For NORMALDIR, negative q means compression and positive q for tension

left.setTracBc(4, task.NORMALDIR, q)

## Predefined consolidation stress
## Arguments: (consolidation_stress_x, consolidation_stress_y)
left.setConsolidation(q, q)

## Bind the material points to material objects
## Arguments: 
## (material_type, material_id, material_points_number_per_element)
left.setMatProperty(task.MCMAT, 1, 2)

## Dampings for material points
## Arguments: (particle_damping, node_damping, PIC_fraction)
left.setDamping(pdamping, gdamping, pic)

right = task.Rectangle(1.75, 1.8, 0.3, 2.3, task.MPTTRACSHAPE)
right.setTracBc(2, task.NORMALDIR, q)
right.setConsolidation(q, q)
right.setMatProperty(task.MCMAT, 1, 2)
right.setDamping(pdamping, gdamping, pic)

## Boundary condition is not always necessary
## NONEBC denotes no boundary BC provided
soil = task.Rectangle(0.85, 1.75, 0.3, 2.3, task.NONEBC)
soil.setMatProperty(task.MCMAT, 1, 2)
soil.setConsolidation(q, q)
soil.setDamping(pdamping, gdamping, pic)

## RIGIDSHAPE claims particular material points: rigid material points
## Rigid material points only serve the boundary condition,
## they ARE NOT involved in any mechanical calculations
## Arguments: (xmin, xmax, ymin, ymax, RIGIDSHAPE)
bot = task.Rectangle(0.7, 1.9, 0.2, 0.3, task.RIGIDSHAPE)

## SetBc controls if the velocity in x or y direction are fixed
## Arguments: (isxfixed, velx_value, isyfixed, vely_value)
bot.setBc(True, 0, True, 0)
bot.setMatProperty(task.RIGIDMAT, 0, 2)

top = task.Rectangle(0.7, 1.9, 2.3, 2.4, task.RIGIDSHAPE)
top.setBc(True, 0, True, -0.02)
top.setMatProperty(task.RIGIDMAT, 0, 2)

## Summary: left and right shapes provide the stress servo
##          and top shape is the shear loading plate
##          bot for the fixed velocity boundary condition
##          soil is the sample undertaking the shear loading

## Shapes array as a member of Analysis object
sim.shapes = [left, right, top, bot, soil]

## 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(15., 1e-3)

## MPM interpolation method,
## options: LINEAR, GIMP, LINEARCPDI, BSPLINECPDI
sim.setMethod(task.GIMP)

## Simulation resutl archive
## When to start the archive (sec)
sim.writer.setStartTime(0)     
 
## Time to stop the archive (sec)                                
sim.writer.setEndTime(17)      

## Frequency for the file archive (sec)                                 
sim.writer.setIntervalTime(0.1)       

## Where to store the archived datas
## Arguments: (directory_name, file_name)                          
sim.writer.setDir("./Results/", 'biaxial')                          
   
## Simulation task, options: MANAL, TANAL, PHASEANAL
## MANAL: mechanical simulation
## TANAL: thermal simulation
## PHASEANAL: phase transition simulation
sim.simTypeMasks = task.MANAL                                    

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

  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) 

  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

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
length = 1.0
eleNum = 20
elemsize = length / eleNum
pic = 2
xrange = [0, length+2.0*elemsize]
yrange = xrange

## Dimension of the grid
## Arguments: (xmin, xmax, ymin, ymax)
sim.grid.setRange(xrange[0], xrange[1], yrange[0], yrange[1])

## Element numbers
## Arguments: (numberx, numbery)
sim.grid.setNumEle(eleNum+2, eleNum+2)
sim.grid.printInfo()


mc=task.Mc()
mc.rho = 1
mc.E = 5e5

## Thermal properties: specific heat and 
## conductivity coefficient tensor, K
## Arguments: (specific_heat, conduction tensor K)

## The use can initialize the K by four arguments [Kxx, Kxy, Kyx, Kyy]
## If only two arguments are given for the tensor
## they are for Kxx and Kyy, [Kxx, 0., 0., Kyy]
mc.setThermalProp(10.0, [1.0,1.0])

## Material array as member of Analysis object
sim.mats=[mc]

## Ids for each material
sim.mats[-1].setIndex(len(sim.mats)-1)

## Shapes:
left = task.Rectangle(0*elemsize, 1.*elemsize, 1*elemsize, (eleNum+1)*elemsize, task.RIGIDSHAPE)

## Thermal boundary condition, temperature
left.setThermalBc(100.0)
left.setMatProperty(task.RIGIDMAT,0, pic)

right = task.Rectangle(xrange[1] - 1.*elemsize, xrange[1]-0*elemsize, 1*elemsize, (eleNum+1)*elemsize, task.RIGIDSHAPE)
right.setThermalBc(100.0)
right.setMatProperty(task.RIGIDMAT, 0, pic)

bot = task.Rectangle(0*elemsize, xrange[1] -0*elemsize, 0*elemsize, 1.*elemsize, task.RIGIDSHAPE)
bot.setThermalBc(100.0)
bot.setMatProperty(task.RIGIDMAT, 0, pic)

top = task.Rectangle(0*elemsize, xrange[1] -0*elemsize, yrange[1]-1*elemsize, yrange[1]-0.*elemsize, task.RIGIDSHAPE)
top.setThermalBc(100.0)
top.setMatProperty(task.RIGIDMAT, 0, pic)

soil = task.Rectangle(1*elemsize, xrange[1] - 1.*elemsize, 1.*elemsize, (eleNum+1)*elemsize, task.NONEBC)
soil.setMatProperty(task.MCMAT, 0, pic)
soil.setDamping(0.2, 0, 0.1)

## Summary: four sides are imposed by fixing temperature boundary condition
## as 100 Celsius, and soil has a initial temperature as 0 Celsius

## Shapes array as a member of Analysis object
sim.shapes = [left, right, bot, top, soil]

## 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, 1e-4)

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

## 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.TANAL

## Where to store the archived datas
## (directory_name, file_name)   
sim.writer.setDir("./Results", "thermal2d")

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